I have around 1000 gbff (genbank) files and I want to extract from each:
- Accession number e.g.
LOCUS NZ_CP011636 - tax ID e.g.
taxon:571 - all the tRNA anticodons (zero or multiple per file) e.g.
/anticodon=pos:complement(1141190..1141192),aa:Met,seq:cat)(I actually just want theaa:Met,seq:catpart)
(note that the latter sometimes spills across two lines).
and compile output into a table containing rows for each file.
Example of file: https://www.ncbi.nlm.nih.gov/nuccore/NZ_CP011636
grep -e "LOCUS" -e "taxon:" -e "anticodon=" Klebsiella_oxytoca.gbff > output.txt
This sort of works but grep doesn't work when #3 spills onto two lines and also I need to store all the output from multiple files. Thanks for any help!
1 answer
This is probably best done using BioPython for a more 'readable' code but if you wish, you can do this using Entrez Direct on the command line using this command
## <filename.txt> should have nucleotide accessions; one per line
epost -db nuccore -input <filename.txt> \
| efetch -db nuccore -format gb -mode xml -style withparts \
| xtract -pattern GBSeq -tab '\n' -element GBSeq_accession-version \
-group GBFeature -if GBFeature_key -equals 'source' \
-block GBQualifier -if GBQualifier_name -equals 'db_xref' \
-tab '\n' -element GBQualifier_value \
-group GBFeature -if GBFeature_key -equals 'tRNA' \
-block GBQualifier -if GBQualifier_name -equals 'anticodon' \
-tab '\n' -element GBQualifier_value
## sample output
NZ_CP011636.1
taxon:571
(pos:complement(60036..60038),aa:Glu,seq:ttc)
(pos:complement(118496..118498),aa:Thr,seq:ggt)
To keep the output size manageable I piped it tohead. Remove that and you will see all 86 features
I have updated the command above to clean it up a little bit, remove the final head command and add an epost step to read nucleotide accessions from a file. Replace <filename.txt> with your file. Good luck!
This sounds like a formatting issue. See what happens if you paste the entire command as a single line without the backslashes.
Log in to answer this question.