DDBJ_MSS形式からGFF3を作る

Biopythonを使ってGFF３ファイルで遊ぶ

GFF3形式のファイルはゲノム上のアノテーション情報を扱う場面などでよく利用されますが、parent 関係があるなど、扱いはちょっと大変です。しかしBiopythonのGFF parsing機能をつかうと、便利に使えるようになりましたので。その扱い方法を記述しておきます。

Biopythonは配列データをSeqRecord形式で扱います。
SeqRecordはSeqとIDの情報を必須として持ち、さらにSeqFeature項目を持つことが出来ます。
Seqは塩基配列情報、IDは配列名です。

例えばSeqRecordデータを１から作る場合、以下の様になります

'''python: make_SeqRecord.py
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord

seq = Seq("GATCGATCGATCGATCGATC")
id= "ID_1"
rec = SeqRecord(seq, id) #Serecored()関数を使ってSeqRecordオブジェクトであるrecを作成
'''

実際はSeqRecordはfasta からつくることが多いと思います。
fastaから作る場合は　SeqIO()関数を使います。
multiple fastaの場合は１配列づつ読み込まれていきます。

'''python: imp_fasta.py
from Bio import SeqIO
for record in SeqIO.parse("example.fasta", "fasta"): #example.fastaファイルを開く。"fasta"で形式を指定
print(record.id)
'''

アノテーション情報はSeqFeature形式で扱われます。
SeqFeatureは　FeatureLocation項目で領域やstrandを指定し、情報タイプ（gene やmRNA等）をtype項目で指定します。他のNote, ID, product などは辞書形式のqualifiers項目で指定します。
具体的には、以下の様な感じです。

'''python: make_feature.py
from BCBio import GFF
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.SeqFeature import SeqFeature, FeatureLocation

qualifiers = {"source": "prediction", "score": 10.0, "product": "hypothetical protein","ID": "gene1"} #辞書形式でqualifersを作成
feature = SeqFeature(FeatureLocation(0, 20), type="gene", strand=1, qualifiers=qualifiers)
'''

これも実際は既に出来ているgenbankデータや、gffファイルから作成することになります。

Ciconia.v5.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-


import pprint
import sys
from Bio import SeqIO
from BCBio import GFF
from BCBio.GFF import GFFExaminer
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.SeqFeature import SeqFeature, FeatureLocation


in_fasta_file  = sys.argv[1]
in_gff_file  = sys.argv[2]
out_file = sys.argv[3]

PreSubfeatureID = "new_contig"


##出力用データを初期化
out_seq_record = SeqRecord(Seq(""),"")


in_handle = open(in_gff_file)
for now_GFF_contig in GFF.parse(in_handle):  ##GFFファイルを開く
    for now_fasta in SeqIO.parse(in_fasta_file, "fasta"):   ##fastaファイルを開く
        if now_GFF_contig.id == now_fasta.id:               ##fastaとqiitaGFFの両ファイルのシーケンスIDが一致したら
            out_fasta_contig = now_fasta #出力用のSeqrecord を作成
            out_fasta_contig.features = [] #contigのfeaturesを初期化
            print("==============================Contig " + now_GFF_contig.id + " started===========================================================")
            for now_GFF_feature in now_GFF_contig.features: #GFFのサブフューチャーを一つ筒読み込む
                if now_GFF_feature.type == "CDS": #CDStyepだった時に
                    print ("found CDS type feature from " + str(now_GFF_feature.qualifiers['Note']))
                    if PreSubfeatureID == "new_contig": #新コンティグだった場合
                        #TopfetureとしてmRNA項目を作成する
                        Count_ID = 1
                        mRNA_location_START = now_GFF_feature.location.start
                        mRNA_location_END = now_GFF_feature.location.end
                        mRNA_location_STRAND = now_GFF_feature.location.strand
                        mRNA_temp_PRODUCT = now_GFF_feature.qualifiers['product']
                        mRNA_temp_ID = now_GFF_feature.qualifiers['Note']
                        mRNA_ID = mRNA_temp_ID[0]
                        mRNA_Note = now_GFF_feature.id
                        print("!!!!!!!!!new mRNA name is ===  " + mRNA_ID)
                        sub_qualifiers = {"source": "prediction", "Note": mRNA_Note,"ID":mRNA_ID, "product":mRNA_temp_PRODUCT}
                        top_mRNA_feature = SeqFeature(FeatureLocation(mRNA_location_START, mRNA_location_END, strand = mRNA_location_STRAND),  id = mRNA_ID, type="mRNA", strand = mRNA_location_STRAND, qualifiers=sub_qualifiers)
                        print(top_mRNA_feature.id)
                        #subfeatureとしてCDS項目を作成しmRNA_featureの下にいれる
                        #CDSデータを元にsubfeatureを作り出す
                        out_GFF_featur = now_GFF_feature
                        CDS_location_Note = now_GFF_feature.qualifiers['Note']
                        out_GFF_featur.qualifiers['ID'] = CDS_location_Note[0] + "_" + str(Count_ID)
                        top_mRNA_feature.sub_features = []
                        top_mRNA_feature.sub_features.append (out_GFF_featur)
                        #PreSubfeatureIDにCDSのNote項目をいれる（同じNote=transcriptのときまとめてmRNAにするため）
                        PreSubfeatureID = now_GFF_feature.qualifiers['Note']
                    elif PreSubfeatureID == now_GFF_feature.qualifiers['Note']: #同じNote=transcriptのとき
                        Count_ID += 1
                        out_GFF_featur = now_GFF_feature
                        CDS_location_Note = now_GFF_feature.qualifiers['Note']
                        out_GFF_featur.qualifiers['ID'] = CDS_location_Note[0] + "_" + str(Count_ID)
                        top_mRNA_feature.sub_features.append (out_GFF_featur) #CDSをsubfeatureとして追加する
                        mRNA_location_END = now_GFF_feature.location.end #mRNAの終了地点を更新する
                        mRNA_temp_PRODUCT = now_GFF_feature.qualifiers['product']
                        mRNA_temp_ID = now_GFF_feature.qualifiers['Note']
                        CDS_location_PRODUCT = now_GFF_feature.qualifiers['product']
                        mRNA_ID = mRNA_temp_ID[0]
                        #mRNA_Note = now_GFF_feature.id
                    elif PreSubfeatureID != now_GFF_feature.qualifiers['Note']: #違うnote=transcriptにあたった場合
                        #今までのmRNAを出力する
                        sub_qualifiers = {"source": "prediction", "Note": mRNA_Note, "ID":mRNA_ID, "product":mRNA_temp_PRODUCT}
                        top_end_mRNA_feature = SeqFeature(FeatureLocation(mRNA_location_START, mRNA_location_END, strand = mRNA_location_STRAND),  id = mRNA_ID, type="mRNA", strand = mRNA_location_STRAND, qualifiers=sub_qualifiers)
                        print(top_end_mRNA_feature)
                        top_end_mRNA_feature.sub_features = top_mRNA_feature.sub_features
                        out_fasta_contig.features.append(top_end_mRNA_feature) #作成したmRNA（＋CDS）をSeqRecoredのfeaturesとして追加する
                        #新しいTopfetureとしてmRNA項目を作成する
                        Count_ID = 1
                        mRNA_location_START = now_GFF_feature.location.start
                        mRNA_location_END = now_GFF_feature.location.end
                        mRNA_location_STRAND = now_GFF_feature.location.strand
                        mRNA_temp_ID = now_GFF_feature.qualifiers['Note']
                        mRNA_ID = mRNA_temp_ID[0]
                        mRNA_temp_PRODUCT = now_GFF_feature.qualifiers['product']
                        mRNA_Note = now_GFF_feature.id
                        out_GFF_featur = now_GFF_feature
                        CDS_location_Note = now_GFF_feature.qualifiers['Note']
                        out_GFF_featur.qualifiers['ID'] = CDS_location_Note[0] + "_" + str(Count_ID)
                        top_mRNA_feature = SeqFeature()
                        top_mRNA_feature.sub_features=[out_GFF_featur] #CDSをsubfeatureとして追加する
                        #PreSubfeatureIDにCDSのNote項目をいれる（同じNote=transcriptのときまとめてmRNAにするため）
                        PreSubfeatureID = now_GFF_feature.qualifiers['Note']
            #最後のmRNAを出力する
            sub_qualifiers = {"source": "prediction", "Note": mRNA_Note,"ID":mRNA_ID, "product":mRNA_temp_PRODUCT}
            top_end_mRNA_feature = SeqFeature(FeatureLocation(mRNA_location_START, mRNA_location_END, strand = mRNA_location_STRAND),  id = mRNA_ID, type="mRNA", strand = mRNA_location_STRAND, qualifiers=sub_qualifiers)
            print(top_end_mRNA_feature)
            top_end_mRNA_feature.sub_features = top_mRNA_feature.sub_features
            out_fasta_contig.features.append(top_end_mRNA_feature) #作成したmRNA（＋CDS）をSeqRecoredのfeaturesとして追加する
            print (out_fasta_contig)
            with open(out_file, "a") as out_handle:
                GFF.write([out_fasta_contig], out_handle)
    print("+++++++++++++Contig " + now_GFF_contig.id + " finished +++++++++++++++++") #1contigが終了
    PreSubfeatureID = "new_contig" #新しいcontigに入った！
in_handle.close