biopythonを用いたシーケンスによるfastaファイルからの遺伝子抽出

Question

私は解析中のfastaファイルを持っています。このファイルは、異なる細菌株由来の同じ遺伝子に属するいくつかの配列からなる。私は何をしたいのですか？スクリプトは、シーケンスが参照seqと等しいかどうかをチェックすることです。その情報で私は新しいファイルを作りたいと思っていますが、それぞれのシーケンスは1つだけです。

def checking_sequences(): 
gene_records=list(SeqIO.parse('/files/gene_A.fasta', 'fasta')) 
ref_id=gene_records[-1].id 
ref_seq=gene_records[-1].seq 
#print gene_records[-1].description 
output_handle=open('//files/' + 'corrected_gene_1', 'a') 
print len(gene_records) 
count=0 
dif_count=0 
reference_list=[] 

for gene_record in gene_records: 
    #count+=1 
    if len(gene_record.seq) == len(ref_seq): 
    #print len(gene_records.seq) 
    #print len(ref_seq) 
     print 'Found all lengths are equal'      
    else: 
     print 'Found %s sequence with different lengths' % (gene_records.description) 

    ###checking sequence equality 
    if gene_record.seq==ref_seq: 
     count+=1 
     gene_record.id=gene_record.id +'_0' 
     reference_list.append(gene_record) 
     ref_count=reference_list.count(gene_record.seq) 
     print 'There are %i sequences are equal to the reference sequence' %(count)  
    else:  
     dif_count+=1 
     reference_list.append(gene_record.seq) 
     seq_count=reference_list.count(gene_record.seq) 
     gene_record.id=gene_record.id +'_'+ str(dif_count) 
     print 'Found %i different that ref_seq' % (seq_count) 
     print 'xxxxxxxxxxx' 




     #print seq_count 
     #print len(reference_list) 
    SeqIO.write(gene_record, output_handle, 'fasta') 


output_handle.close() 

checking_sequences（）いくつかの明確化のために ：

original file       desire output 
    >gene_1 strainIDx      >gen1_strainIDx 
    seqA         seqA 
    >gene_1 strainIDy      >gene_1 strainIDy 
    seqB          seqB 
    >gene_1 strainIDz 
    seqA 

私はちょうど私がそれぞれの1つの配列を持ちたいIDについては気にしません。私は 'ブレーク'を使用しようとしましたが、私が望む出力は得られません。ヘルプ

Answer 1

コメントを認めるようになります。私はそれがないか持たないかわからなかったので、私はハッシュのことを聞いたことがありません。

参考：

組み込み関数hash(object)
は整数としてハッシュ値を返します 。
辞書検索中に辞書キーをすばやく比較するために使用されます。

SO QA Built in python hash() function

---

質問：私は新しいファイルが、それぞれの1つの配列のみを生成します。

には、例えば、ルックアップテーブルを使用します。

lookup = Set() 

with open('/files/' + 'corrected_gene_1', "w") as handle: 
    for record in SeqIO.parse('/files/gene_A.fasta', "fasta"): 
     seq_hash = hash(str(record.seq)) 

     if not seq_hash in lookup: 
      # Not in lookup, save 
      lookup.add(seq_hash) 
      SeqIO.write(records, handle, "fasta") 
     else: 
      # already saved - Skip 
      pass