Parsing a genbank file and outputting specific feature information to a csv using BioPython

Question

So I am trying to parse through a genbank file, extract particular feature information and output that information to a csv file. The example genbank file looks like this:

SBxxxxxx.LargeContigs.gbk

LOCUS       scaffold_31            38809 bp    DNA              UNK 01-JAN-1980
DEFINITION  scaffold_31.
ACCESSION   scaffold_31
VERSION     scaffold_31
KEYWORDS    .
SOURCE      .
      ORGANISM  .
COMMENT     ##antiSMASH-Data-START##
            Version      :: 6.1.1
            Run date     :: 2022-09-21 11:09:55
            ##antiSMASH-Data-END##
FEATURES             Location/Qualifiers
            protocluster    26198..38809
                            /aStool="rule-based-clusters"
                            /category="terpene"
                            /contig_edge="True"
                            /core_location="[36197:37079](-)"
                            /cutoff="20000"
                            /detection_rule="(Terpene_synth or Terpene_synth_C or
                            phytoene_synt or Lycopene_cycl or terpene_cyclase or NapT7
                            or fung_ggpps or fung_ggpps2 or trichodiene_synth or TRI5)"
                            /neighbourhood="10000"
                            /product="terpene"
                            /protocluster_number="1"
                            /tool="antismash"


  

Now for the output file, I want to create a csv with 3 columns. One column will have the Scaffold information (ie. scaffold_31), the second column will have the category value in the protocluster feature (ie. /category = "terpene") and the third column will have the product value in the protocluster feature (ie. /product="terpene")

This is what I have so far for code. I know I can sort through the feature.qualifiers in the protocluster feature to get the category and product. I am not sure how to extract the scaffold information. I am completely new to parsing through gene bank files so have little knowledge in this domain. Thanks in advance for any assitance!

import Bio
from Bio import SeqIO
import os

input = "/Path to SBxxxxxx.LargeContigs.gbk"
output = open("output.csv", "w")

if not os.path.exists(output):
     for record in SeqIO.parse(input, "genbank")
          for feature in record.features:
              if feature.type == "protocluster" and "category" and "product" in feature.qualifiers:
                  outfile = feature.qualifiers["category"][0] + "," + feature.qualifiers["product"][0] + "\n"
                  output.write(outfile)

```

Answer 1

Depending on which field you want to pull the "scaffold_31" text from, you have a few options:

GenBank Text	BioPython object
LOCUS	record.name
DEFINITION	record.description
ACCESSION	record.accessions[0]
VERSION	record.version

Python's built in dir() function is handy for figuring out this kind of thing.

If you're working with a draft flat file (like BankIt gives you just before submitting) note that some of those are placeholders that get updated with the actual accession info when it's finalized. So your "scaffold_31" text will only show up I think in the DEFINITION line in the end if I remember right.

Answer 2

I've used SARS-CoV-2 (Genbank: PA544053), because there was no Genbank entry given in the OPs question. If my example is representative (might not be) I think its about the object attributes. To understand the object I listed its attributes,

print(record[0].__dict__.keys())

Output

dict_keys(['_seq', 'id', 'name', 'description', 'dbxrefs', 'annotations', '_per_letter_annotations', 'features'])

The docs and @jesse's very kind response says there's a 'accession' attribute (Biopython docs below). @Jesse did mention dir() which was cool. I couldn't find record[0].accession or perhaps record[0].accessions and the OP might have had the same problem.

In my example there is an 'annotations' attribute and beneath that was 'accession' accessed via,

print(record[0].annotations['accessions'])

Output

['PA544053']

There it is ... Here's the full code including the CSV package, I'm using efetch so it'll just copy and paste and run.

from Bio import Entrez
from Bio import SeqIO
import csv
from pathlib import Path

fileout = Path('/mypath/mydir', 'test.csv')
handle = Entrez.efetch(db="nucleotide", id="PA544053", rettype="gb",retmode="text") #, 
record = [rec for rec in SeqIO.parse(handle,'genbank')]
fieldnames = ['Accession', 'db_xref', 'organism']

# print(record[0].__dict__.keys())
print(record[0].annotations['accessions'])

with open (fileout, 'w') as f:
    csvwriter = csv.writer(f, delimiter="\t", quotechar='"')
    csvwriter.writerow(fieldnames)    
    for feature in record[0].features:
        if feature.type == 'organism' or "mol_type" and "db_xref" in feature.qualifiers:
            csvwriter.writerow([''.join(record[0].annotations['accessions']), ''.join(feature.qualifiers['db_xref']), ''.join(feature.qualifiers['organism'])])

Output:

Accession      db_xref       organism
PA544053       taxon:2697049 Severe acute respiratory syndrome coronavirus 2

---

Biopython docs A likely reason for the question is the missing attribute is described in the official docs. Using a GenBank object (not SeqIO) there is certainly an accession attribute, https://biopython.org/docs/1.75/api/Bio.GenBank.html

The official Biopython page ...

from Bio import GenBank
with open("GenBank/NC_000932.gb") as handle: 
    record = GenBank.read(handle)    
    print(record.accession)

I think the basis of the question is to associate the accession number with the biochemical/genetic info. Note, I don't know the difference between SeqIO and GenBank objects.