使用biopython的SeqIO解析genbank文件格式

Question

我试图解析一个蛋白质genbank文件格式，下面是一个示例文件（example.protein.gpff）

LOCUS       NP_001346895             208 aa            linear   PRI 20-JAN-2018
DEFINITION  intercellular adhesion molecule 2 precursor [Cercocebus atys].
ACCESSION   NP_001346895
VERSION     NP_001346895.1
DBSOURCE    REFSEQ: accession NM_001359966.1
KEYWORDS    RefSeq.
SOURCE      Cercocebus atys (sooty mangabey)
  ORGANISM  Cercocebus atys
            Eukaryota; Metazoa; Chordata; Craniata; Vertebrata; Euteleostomi;
            Mammalia; Eutheria; Euarchontoglires; Primates; Haplorrhini;
            Catarrhini; Cercopithecidae; Cercopithecinae; Cercocebus.
REFERENCE   1  (residues 1 to 208)
  AUTHORS   Palesch D, Bosinger SE, Tharp GK, Vanderford TH, Paiardini M,
            Chahroudi A, Johnson ZP, Kirchhoff F, Hahn BH, Norgren RB, Patel
            NB, Sodora DL, Dawoud RA, Stewart CB, Seepo SM, Harris RA, Liu Y,
            Raveendran M, Han Y, English A, Thomas GWC, Hahn MW, Pipes L, Mason
            CE, Muzny DM, Gibbs RA, Sauter D, Worley K, Rogers J and Silvestri
            G.
  TITLE     Sooty mangabey genome sequence provides insight into AIDS
            resistance in a natural SIV host
  JOURNAL   Nature 553 (7686), 77-81 (2018)
   PUBMED   29300007
COMMENT     VALIDATED REFSEQ: This record has undergone validation or
            preliminary review. The reference sequence was derived from
            KY308194.1.

            ##Evidence-Data-START##
            Transcript exon combination :: KY308194.1 [ECO:0000332]
            RNAseq introns              :: single sample supports all introns
                                           SAMN02045730, SAMN03085078
                                           [ECO:0000348]
            ##Evidence-Data-END##
FEATURES             Location/Qualifiers
     source          1..208
                     /organism="Cercocebus atys"
                     /db_xref="taxon:9531"
     Protein         1..208
                     /product="intercellular adhesion molecule 2 precursor"
                     /calculated_mol_wt=21138
     sig_peptide     1..19
                     /inference="COORDINATES: ab initio prediction:SignalP:4.0"
                     /calculated_mol_wt=1999
     Region          24..109
                     /region_name="ICAM_N"
                     /note="Intercellular adhesion molecule (ICAM), N-terminal
                     domain; pfam03921"
                     /db_xref="CDD:252248"
     Region          112..>167
                     /region_name="Ig"
                     /note="Immunoglobulin domain; cl11960"
                     /db_xref="CDD:325142"
     CDS             1..208
                     /gene="ICAM2"
                     /coded_by="NM_001359966.1:1..627"
                     /db_xref="GeneID:105590766"
ORIGIN
        1 mssfgfgtlt malfalvccs gsdekafevh mrleklivkp kesfevncst tcnqpevggl
       61 etslnkilll eqtqwkhyli snishdtvlw chftcsgkqk smssnvsvyq pprqvfltlq
      121 ptwvavgksf tiecrvpave pldsltlsll rgsetlhsqt frkaapalpv lrelgmkfiq
      181 lcprrglagt mppsrpwcpa athwsqgc
//
LOCUS       NP_001280013             406 aa            linear   MAM 22-JAN-2018
DEFINITION  26S proteasome regulatory subunit 8 [Dasypus novemcinctus].
ACCESSION   NP_001280013 XP_004456848
VERSION     NP_001280013.1
DBSOURCE    REFSEQ: accession NM_001293084.1
KEYWORDS    RefSeq.
SOURCE      Dasypus novemcinctus (nine-banded armadillo)
  ORGANISM  Dasypus novemcinctus
            Eukaryota; Metazoa; Chordata; Craniata; Vertebrata; Euteleostomi;
            Mammalia; Eutheria; Xenarthra; Cingulata; Dasypodidae; Dasypus.
COMMENT     VALIDATED REFSEQ: This record has undergone validation or
            preliminary review. The reference sequence was derived from
            AAGV03083693.1.
            On May 28, 2014 this sequence version replaced XP_004456848.1.

            Sequence Note: The RefSeq transcript and protein were derived from
            genomic sequence to make the sequence consistent with the reference
            genome assembly. The genomic coordinates used for the transcript
            record were based on alignments.

            ##Evidence-Data-START##
            RNAseq introns :: mixed/partial sample support SAMN00668203,
                              SAMN00744121 [ECO:0000350]
            ##Evidence-Data-END##
FEATURES             Location/Qualifiers
     source          1..406
                     /organism="Dasypus novemcinctus"
                     /db_xref="taxon:9361"
     Protein         1..406
                     /product="26S proteasome regulatory subunit 8"
                     /calculated_mol_wt=45495
     Region          4..404
                     /region_name="RPT1"
                     /note="ATP-dependent 26S proteasome regulatory subunit
                     [Posttranslational modification, protein turnover,
                     chaperones]; COG1222"
                     /db_xref="CDD:224143"
     CDS             1..406
                     /gene="PSMC5"
                     /coded_by="NM_001293084.1:1..1221"
                     /db_xref="GeneID:101445299"
ORIGIN
        1 maldgpeqme leegkagsgl rqyylskiee lqlivndksq nlrrlqaqrn elnakvrllr
       61 eelqllqeqg syvgevvram dkkkvlvkvh pegkfvvdvd knidindvtp ncrvalrnds
      121 ytlhkilpnk vdplvslmmv ekvpdstyem iggldkqike ikevielpvk hpelfealgi
      181 aqpkgvllyg ppgtgktlla ravahhtdct firvsgselv qkfigegarm vrelfvmare
      241 hapsiifmde idsigssrle ggsggdsevq rtmlellnql dgfeatknik vimatnridi
      301 ldsallrpgr idrkiefppp neearldilk ihsrkmnltr ginlrkiael mpgasgaevk
      361 gvcteagmya lrerrvhvtq edfemavakv mqkdseknms ikklwk
//

该格式具有重复的record（由//分隔），其中每个record是一种蛋白质。每个record都有几个部分，其中有一个带有几个固定字段的FEATURES部分，例如source、CDS和Region，这些值引用特定于该record的信息

我感兴趣的是使用biopython的SeqIO将这个文件解析成一个dataframe，其中列出了每个记录ID的gene、db_xref和coded_by的值，以及CDS字段中的organism和^}值，以及Region字段中的{}值。除了Regions字段可能在record的FEATURES部分出现多次外，CDS和source字段在record的FEATURES部分只出现一次

到目前为止，我的失败尝试如下所示：

from Bio import SeqIO
filename = "example.protein.gpff"
for record in SeqIO.parse(filename, "genbank"):
  for feature in record.features:
    if feature.type == "CDS":
      symbol = feature.qualifiers.get("gene", ["???"])[0]
      gene_id = feature.qualifiers.get("db_xref", ["???"])[0]
      gene_id = re.sub('GeneID:', '', gene_id)
      transcript_id = feature.qualifiers.get("coded_by", ["???"])[0]
      transcript_id = re.sub(':.*', '', transcript_id)
      if feature.type == "source":
        species_name = feature.qualifiers.get("organism", ["???"])[0]
        species_id = feature.qualifiers.get("db_xref", ["???"])[0]
        species_id = re.sub('taxon:', '', species_id)
      if feature.type == "Region":
        cdd_id = feature.qualifiers.get("db_xref", ["???"])[0]
        cdd_id = re.sub('CDD:', '', cdd_id)
      print("%s,%s,%s,%s,%s,%s,%s" % (record.id, cdd_id, transcript_id, symbol, gene_id, species_name, species_id))

我想要得到的结果dataframe（对于上面的example.protein.gpff）是：

    record_id  CDS_coded_by   CDS_db_xref CDS_gene              source_organism source_db_xref Region_db_xref
1 NP_001346895  NM_001359966.1:1..627 GeneID:105590766  ICAM2      Cercocebus atys taxon:9531 CDD:252248
2 NP_001346895  NM_001359966.1:1..627 GeneID:105590766  ICAM2      Cercocebus atys taxon:9531 CDD:325142
3 NP_001280013 NM_001293084.1:1..1221 GeneID:101445299  PSMC5 Dasypus novemcinctus taxon:9361 CDD:224143

Answer 1

似乎处理此文件格式的最简单方法是将其转换为JSON格式（例如，使用Bio），然后使用各种JSON解析器读取它（如R中的rjson包，它将JSON文件解析为record中的list）