2016-08-31 8 views
0

私は巨大なfastaファイルを持っていますが、シーケンスの開始と終了の塩基対座標を知っていれば、その一部だけを抽出する必要があります。また、それは1行あたり60 bpの長さのfasta形式でなければなりません。これは私の試みです。もしそれが改善されたと思われるならお知らせください。bp座標に基づいてfastaシーケンスの一部を抽出します

from Bio import SeqIO 

inFile = open('full_chr.fa','r') 
fw=open("part.fa",'w') 
line_width = 60 
for record in SeqIO.parse(inFile,'fasta'): 
    fw.write(">" + record.id + "\n") 
    fww = (str(record.seq[600130000:602000000]) + '\n') 
    for i in xrange(0,len(fww),line_width): 
     fw.write(str(fww[i:i+line_width]) + '\n') 
fw.close() 
+0

FASTAフォーマットがライン –

+0

はいあたり60bpの長さに制限されていません。これは、80本のbpのラインの例ですそれ以外の場合は開かない。 – user3224522

+0

私はbedtools getfastaを推奨しています。http://bedtools.readthedocs.io/en/latest/content/tools/getfasta.html –

答えて

2

それはように簡単です:SeqIO.writeはすでに60文字長の折り返しを使用しています

from Bio import SeqIO 


record = SeqIO.read("Chromosome.fas", "fasta") 

with open("output.fas", "w") as out: 
    SeqIO.write(record[100:500], out, "fasta") 

。行折り返しを操作する場合は、FastaWriterオブジェクトを使用します。 、60塩基対の長さを使用することをお勧めしてgeditでそれを開くために、私にとって

from Bio import SeqIO 
from Bio.SeqIO.FastaIO import FastaWriter 


record = SeqIO.read("Chromosome.fas", "fasta") 

with open("output.fas", "w") as out: 
    writer = FastaWriter(out, wrap=80) 
    writer.write_header() 
    writer.write_record(record[100:500]) 
    writer.write_footer() 
+0

上記のスクリプトはあなたと同じ出力を出力しますが、あなたの出力はよりシンプルです。 – user3224522

関連する問題