3
$\begingroup$

I am using prokka to annotate a bacterial genome:

prokka ecoli.fa

Prokka is outputting a tab-separated file (called PROKKA_12142017.tsv) with differing numbers of elements in each row:

locus_tag    ftype    gene    EC_number    product
BOFHCHGE_00001    CDS    hypothetical protein
BOFHCHGE_00002    CDS    rdgB_1    3.6.1.9    dITP/XTP pyrophosphatase
BOFHCHGE_00003    CDS    hypothetical protein
BOFHCHGE_00004    CDS    hypothetical protein
BOFHCHGE_00005    CDS    hypothetical protein

Does anybody know how I can get round this to import into a dataframe in R?

Is there a way to ask prokka to add something in the gene column?

Improve this question
$\endgroup$
1
  • 1
    $\begingroup$ Please make sure to send a bug report to the prokka authors, since this should be fixed on their end. $\endgroup$
    – Devon Ryan
    Dec 14, 2017 at 15:14

5 Answers 5

Reset to default
5
$\begingroup$

Author here - this was a bug which was fixed in August:

See the issue in Github.

You will need to use the latest Github HEAD version until the next release.

Improve this answer
$\endgroup$
6
$\begingroup$

You don't really need to get around anything, your example can be loaded correctly with read.delim(). The rows with missing values are filled in with blanks. You may, however, prefer the readr package, which handles this a bit more elegantly (it'll tell you which values were missing and not do the annoying string->factor conversion).

Improve this answer
$\endgroup$
5
  • $\begingroup$ Thanks, useful but it puts 'hypothetical protein' under the gene column when it actually belongs in the product column. The gene column should be empty for hypothetical proteins. $\endgroup$
    – Biomagician
    Dec 14, 2017 at 15:01
  • $\begingroup$ Oh, I thought that was intentional. In that case prokka has a bug and you'll need to test for "" as the product and just assign accordingly. $\endgroup$
    – Devon Ryan
    Dec 14, 2017 at 15:02
  • $\begingroup$ And I get 'Sec-independent protein translocase protein TatA' for EC number, so that's definitely wrong. I never see an empty element before a non-empty element. $\endgroup$
    – Biomagician
    Dec 14, 2017 at 16:05
  • 3
    $\begingroup$ You can use the stringsAsFactors = FALSE to prevent converting strings to factors in read.delim/read.csv... $\endgroup$
    – llrs
    Dec 14, 2017 at 21:09
  • 1
    $\begingroup$ @Llopis Quite true, but the readr package has that as the default :) $\endgroup$
    – Devon Ryan
    Dec 14, 2017 at 21:12
3
$\begingroup$

As Devon said, read.delim can deal with this perfectly well as long as the file is properly formatted (so missing fields still end with the field separator, like \t\t). For example, with this input:

Field1  Field2  Field3
Row1field1  Row1field2  Row1field3
Row2field1      Row2field3

Or, to show the fields more clearly:

Field1\tField2\tField3
Row1field1\tRow1field2\tRow1field3
Row2field1\t\tRow2field3

R can import:

> df<-read.delim("file",header=TRUE, sep="\t")
> df
      Field1     Field2     Field3
1 Row1field1 Row1field2 Row1field3
2 Row2field1            Row2field3

Alternatively, if you really want to add something to the missing columns, you can use:

$ awk -F"\t" '($3==""){$3="none" }1;' file
locus_tag   ftype   gene    EC_number   product
BOFHCHGE_00001  CDS hypothetical    protein
BOFHCHGE_00002  CDS rdgB_1  3.6.1.9 dITP/XTP    pyrophosphatase
BOFHCHGE_00003  CDS hypothetical    protein
BOFHCHGE_00004  CDS hypothetical    protein
BOFHCHGE_00005  CDS hypothetical    protein

Or, to add none to any empty field, not just the 3rd one:

awk -F"\t" '{for(i=1;i<=NF;i++){if($i==""){$i="none"; }}}1;' file > newfile
Improve this answer
$\endgroup$
10
  • $\begingroup$ I fear that the .tsv file is mal-formatted then. df<-read.delim("file",header=TRUE, sep="\t") will misimport and EC column will end up having product names, gene having hypothetical protein, etc. $\endgroup$
    – Biomagician
    Dec 14, 2017 at 15:34
  • 1
    $\begingroup$ @charlesdarwin if it is indeed malformatted, yes. However, I would be very surprised if that were the case. Using empty fields as described above is standard practice. You can check easily enough though. FInd one of the lines with a missing field (using 12 here as an example), and then run awk 'NR==12' | od -c and make sure you see the \t\t where the empty field is. If you don't see that, this becomes a text parsing problem. I suggest you post a new question with some example lines and we'll see what we can do. $\endgroup$
    – terdon ♦
    Dec 14, 2017 at 15:39
  • $\begingroup$ my 6th line has a missing EC number. My file is file.txt and my data frame is called df, what is the command I should run in the terminal with awk? $\endgroup$
    – Biomagician
    Dec 14, 2017 at 15:46
  • 1
    $\begingroup$ @charlesdarwin (sorry, I just saw this for some reason) it depends. If there is a \t\t for the empty field, you can do awk -F"\t" '($4==""){$4="none" }1;' file.txt > newfile.txt. If that doesn't work, you might need to ask a new question, or come into the site chatroom and ping me and I'll see if I can help. $\endgroup$
    – terdon ♦
    Dec 15, 2017 at 10:45
  • $\begingroup$ newfile.txt is a mess. The element 1,1 is now BOFHCHGE_00001 CDS hypothetical protein none. I've written to the developer. I don't think it should be that difficult. I am using the GFF file instead. It is well formatted and contains all the information, too. Thanks for your help. $\endgroup$
    – Biomagician
    Dec 15, 2017 at 12:07
1
$\begingroup$

Since the missing values don't have field separators, the easiest option is to rearrange the columns so that the ones with missing values are at the end. They will then be the ones that are filled with missing values when you import them into R. You can use use fill=TRUE with read.table() or use readr, which will complain about missing columns, but fill them in with missing values automatically:

awk -F'\t' -v OFS='\t' '{print($1,$2,$5,$3,$4)}' file.txt |perl -pe 's/\t+/\t/g' >newfile.txt

With readr:

read_tsv("newfile.txt")

# A tibble: 5 x 5
       locus_tag ftype                  product   gene EC_number
           <chr> <chr>                    <chr>  <chr>     <chr>
1 BOFHCHGE_00001   CDS     hypothetical protein   <NA>      <NA>
2 BOFHCHGE_00002   CDS dITP/XTP pyrophosphatase rdgB_1   3.6.1.9
3 BOFHCHGE_00003   CDS     hypothetical protein   <NA>      <NA>
4 BOFHCHGE_00004   CDS     hypothetical protein   <NA>      <NA>
5 BOFHCHGE_00005   CDS     hypothetical protein   <NA>      <NA>

The perl regex is necessary because awk adds tabs after the missing values, which defeats the point of reordering the columns

Improve this answer
$\endgroup$
2
  • $\begingroup$ I did as you said but the data frame is still wrong. On line 6, gene and product are interchanged. $\endgroup$
    – Biomagician
    Dec 17, 2017 at 18:12
  • $\begingroup$ Since you only provided 5 lines in the example, I can't comment on the results for the 6th line $\endgroup$
    – heathobrien
    Dec 18, 2017 at 19:49
1
$\begingroup$

all. I know this question was asked 2 years ago, but I wrote a script in R, with a little bash, to parse the tsv file in the old version of prokka. In the new version of prokka this problem is fixed.

First use bash to get the rows where ftype is CDS.

grep "CDS" all_prokka_annotations.tsv > all_prokka_annotations_CDS.tsv

Next we can parse the TSV with the CDS lines in R.

prokka_annotations = readLines("all_prokka_annotations_CDS.tsv")
prokka_annotations_split = strsplit(prokka_annotations,split="\t")
# first get the protein product name. 
prokka_annotations_product = sapply(prokka_annotations_split, function(x) x[length(x)])
# now the I have all protein product names remove them
prokka_annotations_split = sapply(prokka_annotations_split, function(x) x[-length(x)])
#next get locus tag and ftype
locus_tag_cds = sapply(prokka_annotations_split, function(x) x[c(1,2)])
locus_tag_cds = t(locus_tag_cds)
# now remove locus tag and ftype
prokka_annotations_split = sapply(prokka_annotations_split, function(x) x[-c(1,2)])
# next get EC number
ec_number = lapply(prokka_annotations_split, function(x) x[grep(".",x,fixed=T)])
ec_number = as.character(ec_number)
ec_number[ec_number=="character(0)"] <- NA
## now get the gene names if available. to do this remove the EC numbers. whats left are gene names
prokka_gene_names = lapply(prokka_annotations_split, function(x) x[!grepl(".",x,fixed=T)])
prokka_gene_names = as.character(prokka_gene_names)
prokka_gene_names[prokka_gene_names=="character(0)"] <- NA
prokka_parsed_tsv_correct = cbind(locus_tag_cds,product=prokka_annotations_product,gene=prokka_gene_names,EC=ec_number)

```
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.