BiopythonでGenbankファイル操作：CDSをすべてFASTAに書き出す

Biopython：Genbank file parsing cheat sheetというタイトルで前回書きましたが、一度山に登って見渡した時に広がる風景をあとから書いたようなもので、これから登ってみようという場合には、自分の経験から照らし合わせても、あまり役に立たなかったかもしれません。

BiopythonによるGenbankファイルのparsingでは、情報がどのような構造に格納されるのかを理解することがカギとなります。どんな収納棚（object）があって、その下にどのような構造の引き出し（attributes）が用意されているか。引き出しはどうやって呼び出すか、またその中身は順序（index）で指定できるのか（list）、名前（key）で呼べるのか（dict）。いま自分が見ているものが何なのかを知るには、type()関数を使うことです。()の中に知りたい変数を入れてどんどん調べましょう。

さて、このポストでは、FEATURESセクションでCDSとしてアノテートされている部分について、片っ端から該当領域を抜き出してFASTAに書き出す、をやってみましょう。

早速コードに飛び込みたいところですが、準備運動しましょう。そのほうがごはんが進み、もとい、理解がはかどります。

まずは、FEATURESセクションの構造をとりあげて、各部分の名前を理解します。試しに比較的小さめの、ヒトのミトコンドリアゲノム（NC_012920.1）を例にして、FEATUREセクションにズームインしてみましょう。

FEATURES             Location/Qualifiers
     source          1..16569
                     /organism="Homo sapiens"
                     /organelle="mitochondrion"
                     /mol_type="genomic DNA"
                     /isolation_source="caucasian"
                     /db_xref="taxon:9606"
                     /tissue_type="placenta"
                     /country="United Kingdom: Great Britain"
                     /note="this is the rCRS"
     D-loop          complement(join(16024..16569,1..576))

              ##### 途中省略 #####

     gene            3307..4262
                     /gene="ND1"
                     /gene_synonym="MTND1"
                     /nomenclature="Official Symbol: MT-ND1 | Name:
                     mitochondrially encoded NADH dehydrogenase 1 | Provided
                     by: HGNC:HGNC:7455"
                     /db_xref="GeneID:4535"
                     /db_xref="HGNC:HGNC:7455"
                     /db_xref="MIM:516000"
     CDS             3307..4262
                     /gene="ND1"
                     /gene_synonym="MTND1"
                     /note="NADH dehydrogenase, subunit 1 (complex I); TAA stop
                     codon is completed by the addition of 3' A residues to the
                     mRNA"
                     /codon_start=1
                     /transl_except=(pos:4261..4262,aa:TERM)
                     /transl_table=2
                     /product="NADH dehydrogenase subunit 1"
                     /protein_id="YP_003024026.1"
                     /db_xref="GeneID:4535"
                     /db_xref="HGNC:HGNC:7455"
                     /db_xref="MIM:516000"
                     /translation="MPMANLLLLIVPILIAMAFLMLTERKILGYMQLRKGPNVVGPYG
                     LLQPFADAMKLFTKEPLKPATSTITLYITAPTLALTIALLLWTPLPMPNPLVNLNLGL
                     LFILATSSLAVYSILWSGWASNSNYALIGALRAVAQTISYEVTLAIILLSTLLMSGSF
                     NLSTLITTQEHLWLLLPSWPLAMMWFISTLAETNRTPFDLAEGESELVSGFNIEYAAG
                     PFALFFMAEYTNIIMMNTLTTTIFLGTTYDALSPELYTTYFVTKTLLLTSLFLWIRTA
                     YPRFRYDQLMHLLWKNFLPLTLALLMWYVSMPITISSIPPQT"

              ##### 以下省略 #####

まず、FEATURES以下、左側にsource、D-loop、gene、CDSが見えます。これはfeatureの名前（typeといいます）をあらわしていて、それぞれが個別のアノテーションブロックを形成しています。例えばCDSの行を右に見ていくと、3307..4262があります。アノテーションは配列の部分に対して与えられるものなので、どこなの？についての答えがまず書かれています。位置情報以下、スラッシュ/以降、=の前までが、qualifierで、=の後ろがvalueです。qualifierにはgene、codon_start、translationなどいろいろなものが見えますね。一つのアノテーションブロックにつき、複数のqualifierがあって、それぞれにvalueがぶら下がっているという形です。わかったような、わからないような。。。

ひとつのアノテーションブロックには、下のような情報がぶらさがっています。

• 名前 ==> CDS, D-loopなどがひとつ
• location ==> 3307..4262などがひとつ
• /qualifier=valueがいくつか

ゆっくりズームアウトしてから、コードに入っていきます。
このようなアノテーションブロックがいくつも集まって、FEATURESセクションを形成しています。ひとつのGenbank エントリー（record）には、ひとつだけFEATURESセクションがあります。Genbankファイル自体には、複数のrecordを含めることもできることを覚えておいてください。

さて、ラジオ体操第一が終わりました。

必要なパーツを部分的に書いていきましょう。
（最後にまとめたものをのせます）

まず、そのままでは身軽（すぎ）なpythonさんに、今回の作業で必要な道具箱を持ってくるように伝えます。

import sys
import os
from Bio import SeqIO
from Bio.SeqRecord import SeqRecord

処理するファイル（path）は、コマンドラインから渡せるほうが便利です。最初のargumentとして拾うようにします。

gb_file = sys.argv[1]

Genbankファイルを開いて、中にあるrecordをひとつづつ取り上げていきます。ひとつしかないよ、という方も、がまんしてください。このほうが汎用性があるんです。for ループを使います。recordをひとつづつ拾いながら、このループの直下で以下のことをしておきましょう。

• 配列全体をとる
• Accessionをとる

for record in SeqIO.parse(gb_file, 'genbank'):
    # get entire sequence
    parental_seq = record.seq

    # get accession
    acc = record.id

ここからです。やりたいことは、CDS該当配列を抜きだす、です。FEATURESセクションにあるアノテーションブロックたちは、record.featuresで取り出せます。この中に複数のアノテーションブロック情報（feature）が入っています。二段目のforループを用意します。

また、featureの名前はfeature.typeで、locationはfeature.location、qualifier=valueはfeature.qualifiers（辞書として返ります）でアクセスできます。

名前がCDSのものにつき、locationをとり、これを使って部分配列を取得します。

for feature in record.features:
    # pick only 'CDS' annotation
    if feature.type == 'CDS':
        # get sequence slice with location
        seq_slice = feature.location.extract(parental_seq)

最後の行がナイスですね。これ一行で全体の配列から対応する部分配列を抜き出せます。

さあ、ここまで来たら向こうの島まで泳ごう、というところまできました。あとは何が必要でしょうか。結果はfastaファイルとして得たいので、考えることは以下の通り。

• 出力ファイル名
  • 入力Genbankファイルに複数のrecordが入っていても個別に出力できるよう、accessionを冠した出力ファイルを作りたい。また、見つかったCDSの個数もファイル名に含めると便利かもしれない。
• 抜き出した該当部分のfasta id line
  • ここはご自由に、という感じですが、なんのCDSなのよ、がわからないと困るので、また、fasta id lineがそれぞれuniqueになっていないと後々トラブルのもとになるので、qualifierからgene=、product=に対応する部分を含め、位置とstrand情報も含めたい。

方針が固まったところで、最終コードです。

import sys
from Bio import SeqIO
from Bio.SeqRecord import SeqRecord

# command line args
gb_file = sys.argv[1] # input gb file path


def makeSeqRecord(seq_obj, id_for_fasta):
    record = SeqRecord(seq_obj, id=id_for_fasta, description="")
    return record


for record in SeqIO.parse(gb_file, 'genbank'):
    # basket to collect CDS slices
    CDS_records = []
    # get entire seq
    parental_seq = record.seq # Seq obj
    # get accession
    acc = record.id

    for feature in record.features:
        # pick only CDS
        if feature.type == 'CDS':
            # get seq_slice with location obj
            seq_slice = feature.location.extract(parental_seq)

            # strand info can be accessed by
            strand_int = feature.location.strand
            if strand_int == 1:
                strand = 'top'
            elif strand_int == -1:
                strand = 'btm'
            else:
                strand = 'cannotbe'
            # get start and end
            # As ponnhide-san pointed out (thank you, ponnhide-san!), 
            # this is better in case the start and end
            # is before and after the position 1 (circular genome).
            start = feature.location.parts[0].start.position + 1
            end   = feature.location.parts[-1].end.position
            #start = feature.location.start.position + 1
            #end = feature.location.end.position

            # get gene annotation
            if 'gene' in feature.qualifiers:
                gene = feature.qualifiers['gene'][0] # value is in a list
            else:
                gene = 'NA'

            # get product annotation
            if 'product' in feature.qualifiers:
                product = feature.qualifiers['product'][0]
            else:
                product = 'NA'

            # construct fasta id line
            CDS_id_line = acc + '_' + gene + ' ' + product
            CDS_id_line += ' (' + str(start) + '..' + str(end) + ', strand=' + strand + ')'

            # construct SeqRecord object and put it in a busket
            CDS_records.append(makeSeqRecord(seq_slice, CDS_id_line))

    # output file name construction
    check = 1
    CDS_count = len(CDS_records)
    if CDS_count == 1:
        out_fasta = acc + '_all-CDS_1.fa'
    elif CDS_count >=2:
        out_fasta = acc + '_all-CDS_' + str(CDS_count) + '.mfa'
    else: # por si acaso...
        out_fasta = acc + '_NO-CDS-FOUND.txt'
        check = 0

    # write CDS_records into a file
    with open(out_fasta, mode='w') as ofa:
        if check == 1:
            SeqIO.write(CDS_records, ofa, 'fasta')
        else:
            ofa.write('NO CDS ANNOTATION FOUND.')

出力ファイル名が、おせっかいすぎかもしれません。.faでも.mfaでもどっちでもええやん、という場合にはそうしてください。Fastaについては、私は一つなのかマルチなのか、区別しています。まったくCDSが見つからなかった場合にも、だんまりを決め込むのではなく、.txtファイルを出すようにしています。この類の処理では、あとからみて何が起こっていたのか記録を残すことが大切なことが多いので、見てすぐわかるようにしてみました。処理のログを別ファイルにだすのもいいかもしれません。

まとめ

• recordをひとつひとつとりあげる

• featureたち（日本語複数形）をrecord.featuresでつかむ

• ひとつひとつのfeatureにつき

  • feature.type
  • feature.location
  • feature.qualifiers
    がもれなくついてくるので、ほしいものにアクセスして処理に使う。

• seq_slice = feature.location.extract(parental_seq)が便利。

このあと、Biopython：Genbank file parsing cheat sheetをもう一度みていただくと、もう少し味がすると思います。

Genbank file parsingはやっかいです。たまにしかやらないとすぐに忘れます。そのときのためのヒントになれば幸いです。