So I'm trying to actively learn how to use Biopython, and python for that matter. I made a script that extracts the Genbank accession number and the country from a .gb (genbank) file and puts the accession number and country in a text file. Here is my script:
from Bio import SeqIO
gb_file = "denv1.gb"
for gb_record in SeqIO.parse(open(gb_file,"r"), "genbank") :
printgb_record.name),
for feat in gb_record.features:
if feat.type == 'source':
source = gb_record.features[0]
for qualifiers in source.qualifiers:
if qualifiers == 'country':
country = source.qualifiers['country']
print(country[0])
Now the output text file, every accession number with a country is listed perfectly, but when a country is not present for a particular sequence in the genbank file, it does this:
Instead I would like it to be blank where the name of the country would go. Hopefully this makes sense and isn't too confusing. Can anyone help me out? As I said I'm trying to learn, so pointers or links to help would be greatly appreciated!
1 answer
Ok here's what I've ended up with:
from Bio import SeqIO
gb_file = "test.gbk"
# parse genbank file
for gb_record in SeqIO.parse(gb_file, "genbank") :
# now do something with the record
for feat in gb_record.features:
if feat.type == 'source':
source = gb_record.features[0]
for qualifiers in source.qualifiers:
if qualifiers == 'country':
country = str(source.qualifiers['country'][0])
elif country == '':
country = 'none'
printgb_record.name + ' ' + country)
This works on the test data. The file as you provided it returns:
KT279761 Haiti
EU482591 USA: Puerto Rico
HQ332181 Venezuela
FJ810415 Venezuela: Aragua
If I manually change one of the gbk entries so that it doesn't contain any country data:
##Assembly-Data-END##
FEATURES Location/Qualifiers
source 1..10735
/organism="Dengue virus 1"
/mol_type="viral cRNA"
/strain="Haiti/1207/2014"
/serotype="1"
/host="Homo sapiens"
/db_xref="taxon:11053"
/country=""
/collection_date="27-Nov-2014"
5'UTR 1..94
The same code will return:
KT279761 none
EU482591 USA: Puerto Rico
HQ332181 Venezuela
FJ810415 Venezuela: Aragua
However, without seeing one of your genbanks which doesnt contain the country information I'm not sure if this will work 100% of the time. If the /country= field is missing altogether I think it may still break. For your own practice you may want to alter the code to try and handle that (another simple elif statement would probably do it).
Log in to answer this question.
-
Getting cds files from ncbi using accession numbers
written by lily •
-
How can I print and write the strain /isolate/voucher number of a SeqRecord objec in biopython?
written by diogobio.dh •
-
Unable to extract country name from genbank file
written by Akbar •
-
Getting proteins from genomic interval of a genbank file
written by fhsantanna
-
Retrieve qualifiers plasmid and pathovar
written by felipelira3 •
-
Extracting Gene nucleotide sequences from a Genbank files using Biopython
written by tpaisie
-
Filtering sequences from a Fasta file using Biopython
written by tpaisie
-
If GI number in list1 is present in genbank file, make new list and append location information
written by pawlowac
-
How To Extract Title And /Isolation_Source From A List Of Genbank Accession Numbers
written by samlambrechts299
-
Refseq Genbank To Fasta Format Failing With Contig Fields
written by Käpt'N Silico •
Could you format the code using in the biostar post editor ?
There! Thanks for letting me know!
Couple of other things to point out. with
SeqIO.parse()you don't also need to callopen(). My suggestion would be to throw in a check for existence of location... have a look at this thread from the other day: C: How do can I use Biopython and SeqIO to parse out multiple genes from several NCOn line 5, you're making sure that feature exists before you do something with it. My guess is that your code is failing to detect that a source isnt found in the next record, and populating it with the last thing source evaluated to from a previous record.
So I need to add something to my code to tell it when no country is present in the file it should print nothing? Are you showing that in the above code?
I haven't tested your code but that would be my guess. I could be wrong :p
The code above is similar to what you're doing. It extracts info from Genbanks.
I believe you need a statement something like:
Can you upload your genbank that you're trying this on so I can replicate what you've got?
So the entire file is to large, but there is a shortened version of the file: http://gobin.io/OBcO
ok so i added this to my script:
and it does leave a the accession numbers country column blank but it also adds an extra line inbetween each accession number and country:
I suspect this will be to do with carry over of a newline character in the case where the source can actually be found from the file, but not if it's blank.