BiopythonのSeqIOを用いてfasta形式のファイルから繰り返し遺伝子の配列を取り出す

biopythonのSeqIOを用いてfastaファイルから部分的な配列を取り出したいです

gtfファイルから遺伝子の領域の範囲(例えば、500から600のような数字）を取り出して、それをもとにfastaファイルからdescriptionを指定して配列部分の該当箇所(今回は500から600)の配列を得たいです。

その際、配列の種類が2つあるので条件によって処理の一部を分けようと思い、関数内でSeqIO.parseを使いました。
ところが、forループの1回目は上手くいくにもかかわらず2回目以降は関数内の「for record SeqIO.parse(fasta_in, 'fasta'):」のループを素通りしています。

SeqIOはもしかしてforループの中で繰り返して使ったりはできないのでしょうか？

発生している問題・エラーメッセージ

UnboundLocalError: local variable 'Sequence' referenced before assignment

該当のソースコード

Python(Biopython）


from sys import *
from Bio.Seq import Seq
from Bio import SeqIO


#--------------------
# 1.GTFファイルからid, starand, frame, CDSの始まりと終わりを取ってきて2次元配列にする
# 2.>id strand frameの次の行からCDSの配列をfastaファイルから取ってきてprintする
#--------------------


#--------------------
# make functions
#--------------------

def connect_strand(small, fasta_in): # small = large[i], fasta_in = fasta_file
    A = ' '.join(small)
    Desc = '>' + ' ' + A
    for record in SeqIO.parse(fasta_in, 'fasta'):
        print(record.seq)
        if small[0] in record.description:
            # print('AAAAA')
            BtoE = small[3:]
            Sequence = ''
            for j in range(0, len(BtoE), 2):    # j番目が始まり、j+1番目が終わり
                B = int(BtoE[j])
                E = int(BtoE[j+1])
                Sequence += record.seq[B-1:E]                        
            
        else:
            pass

    return Desc, Sequence


#-------------------
# get information of CDS place
#-------------------
with open('GTF_file.txt', 'r') as gtf_file, open('FASTA_file.txt', 'r') as fasta_file:
    counter = 0
    large = []
    
    for line in gtf_file:
        tmp = line.strip().split()
        if ((tmp[2] == "start_codon") or (tmp[2] == "stop_codon")) and counter ==0:
            small = []
            small.append(tmp[0])   # chromosome
            small.append(tmp[6])   # strand
            small.append(tmp[7])   # frame
            counter = 1

        elif tmp[2] == "CDS" and counter ==1:
            small.append(tmp[3])    # begin
            small.append(tmp[4])    # end

        elif ((tmp[2] == "start_codon") or (tmp[2] == "stop_codon")) and counter ==1:
            counter = 0
            large.append(small)

        else:
            pass

    # large: [[chromosome(fasta desc), strand(+ or -), frame(0, 1 or 2), begin, end,..., begin, end],..., [chromosome, ..., begin, end]]


    for i in range(len(large)):     # ここから該当箇所
        if large[i][1] == '+':   # forward strand
            D, S = connect_strand(large[i], fasta_file)
            print(D)
            print(S)
        elif large[i][1] == '-':    # reverse strand
            D, S = connect_strand(large[i], fasta_file)
            S = S.reverse_complement()
            print(D)
            print(S)
        else:
            pass

試したこと

「for i in range(len(large)): # ここから該当箇所」以降の部分を関数を使用せずに同じ内容を繰り返して書いたり、予約語のチェックもしたりしました。

補足情報（FW/ツールのバージョンなど）

python3.7

行動規範の内容に同意します

回答1件

ベストアンサー

ところが、forループの1回目は上手くいくにもかかわらず2回目以降は関数内の「for record SeqIO.parse(fasta_in, 'fasta'):」のループを素通りしています。
SeqIOはもしかしてforループの中で繰り返して使ったりはできないのでしょうか？

ライブラリの実装によりますが、SeqIO.parseにて渡されたファイルを全読込しているとファイル位置が末尾に達します。
その結果、次回の呼び出し以降はデータなしになりえます。
よってconnect_strand関数の最初にfasta_in.seek(0)にてファイル位置を先頭に戻すことで解消される可能性があります。

投稿2021/09/03 06:58