0

今のところ、それを行うために独自の関数を定義して文書化しようとしましたが、コードのテストで問題が発生しており、それが正しいかどうかは実際にはわかりません. BioPython、re、またはその他でいくつかの解決策を見つけましたが、これを歩留まりで機能させたいと本当に思っています。

#generator for GenBank to FASTA
def parse_GB_to_FASTA (lines):
    #set Default label
    curr_label = None
    #set Default sequence
    curr_seq = ""
    for line in lines:
        #if the line starts with ACCESSION this should be saved as the beginning of the label
        if line.startswith('ACCESSION'):
            #if the label has already been changed
            if curr_label is not None:
                #output the label and sequence
                yield curr_label, curr_seq
                ''' if the label starts with ACCESSION, immediately replace the current label with
                the next ACCESSION number and continue with the next check'''
            #strip the first column and leave the number
            curr_label = '>' + line.strip()[12:]
        #check for the organism column
        elif line.startswith ('  ORGANISM'):
            #add the organism name to the label line
            curr_label = curr_label + " " + line.strip()[12:]
        #check if the region of the sequence starts
        elif line.startswith ('ORIGIN'):
            #until the end of the sequence is reached
            while line.startswith ('//') is False:
                #get a line without spaces and numbers
                curr_seq += line.upper().strip()[12:].translate(None, '1234567890 ')
    #if no more lines, then give the last label and sequence            
    yield curr_label, curr_seq
4

1 に答える 1

0

私は非常に大きな GenBank ファイルを扱うことが多く、(何年も前に) BioPython パーサーが脆すぎて、(当時) 数十万件のレコードを異常なレコードでクラッシュすることなく処理できないことがわかりました。

開いているファイルから次のレコード全体を返す純粋な python(2) 関数を作成し、1k チャンクを読み取り、次のレコードを取得できるようにファイル ポインタを残します。これを、この関数を使用する単純な反復子と、fasta バージョンを取得する fasta(self) メソッドを持つ GenBank Record クラスに結び付けました。

YMMV ですが、次のレコードを取得する関数はここにあり、使用するイテレータ スキームでプラグインできるはずです。fasta への変換に関する限り、上記の ACCESSION と ORIGIN の取得と同様のロジックを使用できます。または、次を使用してセクションのテキスト (ORIGIN など) を取得できます。

sectionTitle='ORIGIN'    
searchRslt=re.search(r'^(%s.+?)^\S'%sectionTitle,
                     gbrText,re.MULTILINE | re.DOTALL) 
sectionText=searchRslt.groups()[0]

ORGANISM のようなサブセクションでは、左側に 5 つのスペースのパッドが必要です。

主な問題に対する私の解決策は次のとおりです。

def getNextRecordFromOpenFile(fHandle):
    """Look in file for the next GenBank record
    return text of the record
    """
    cSize =1024
    recFound = False
    recChunks = []
    try:
        fHandle.seek(-1,1)
    except IOError:
        pass
    sPos = fHandle.tell()

    gbr=None
    while True:
        cPos=fHandle.tell()
        c=fHandle.read(cSize)
        if c=='':
            return None
        if not recFound:

            locusPos=c.find('\nLOCUS')
            if sPos==0 and c.startswith('LOCUS'):
                locusPos=0
            elif locusPos == -1:
                continue
            if locusPos>0:
                locusPos+=1
            c=c[locusPos:]
            recFound=True
        else:
            locusPos=0

        if (len(recChunks)>0 and
            ((c.startswith('//\n') and recChunks[-1].endswith('\n'))
             or (c.startswith('\n') and recChunks[-1].endswith('\n//'))
             or (c.startswith('/\n') and recChunks[-1].endswith('\n/'))
             )):
            eorPos=0
        else:
            eorPos=c.find('\n//\n',locusPos)

        if eorPos == -1:
            recChunks.append(c)
        else:
            recChunks.append(c[:(eorPos+4)])
            gbrText=''.join(recChunks)
            fHandle.seek(cPos-locusPos+eorPos)
            return gbrText
于 2014-11-05T01:16:27.357 に答える