0

Genbank で GI 番号を段階的に廃止し、次の形式でヘッダーを編集した場所にいくつかの fasta ファイルが保存されているという警告が表示され続けます。

>SomeText_ginumber

どこから始めればよいかわかりませんが、理想的にはPythonを使用して、NCBIから各giに対応するアクセッション番号を取得し、次のようにヘッダー付きのファイルを出力できる方法はありますか?

>SomeText_accessionnumber

ファイル形式の別の例を次に示します。

>Desulfovibrio_fructosivorans_ferredoxin_492837709
MTKSVAPTVTIDGKVVPIEGERNLLELIRKVRIDLPTFCYHSELSVYGACRLCLVEVKNRGIMGAC
>Oxalobacteraceae_bacterium_AB_14_hypothetical_protein_522195384
MIKLTVNGIPVEVDEGATYLDAANKAGVHIPTLCYHPRFRSHAVCRMCLVHVAGSSRPQAACIGKA

編集/更新:

from Bio import Entrez
from time import sleep
import sys
import re
Entrez.email = ''

seqs = open(sys.argv[1],"r")

for line in seqs:
    if line.startswith('>'):
        gi = re.findall("\d{5,}",line)
        matches = len(gi)
        #print(matches)
        if matches < 2:
            if gi:
                handle = Entrez.efetch(db="nucleotide", id=gi, retmode="xml")
                records = Entrez.read(handle)
                print(line[0:line.rfind('_') + 1] + records[0]['GBSeq_primary-accession'])
                sleep(1)
        elif matches >= 2:
             print("Error - More than one ginumber in header!")
             break
        else:
            seq=line.rstrip()
            print(seq)
    else:
        seq1=line.rstrip()
        print(seq1)
4

1 に答える 1

1

BioPythonを使用してみてください。

次のスニペットで開始できます。最初にヘッダー (ヘッダーのアンダースコアの後の部分) から GI を取得し、GenBank からデータを取得し、古いヘッダーを印刷しますが、アクセッション番号と残りの入力シーケンスを印刷します:)

これは 2 つの例で機能しますが、データが増えると失敗する可能性があります (GI の欠落など)。また、ヘッダーと同様に、アクセッション番号にはアンダースコアが含まれているため、後で解析が複雑になります。おそらく、アンダースコアを別のものに置き換えるか、別のセパレーターを追加してください。

from Bio import Entrez
from time import sleep
Entrez.email = 'your@email'

seqs = """>Desulfovibrio_fructosivorans_ferredoxin_492837709
MTKSVAPTVTIDGKVVPIEGERNLLELIRKVRIDLPTFCYHSELSVYGACRLCLVEVKNRGIMGAC
>Oxalobacteraceae_bacterium_AB_14_hypothetical_protein_522195384
MIKLTVNGIPVEVDEGATYLDAANKAGVHIPTLCYHPRFRSHAVCRMCLVHVAGSSRPQAACIGKA"""


for line in seqs.splitlines():
    if line.startswith('>'):
        gi = line[line.rfind('_') + 1:]
        handle = Entrez.efetch(db="nucleotide", id=gi, retmode="xml")
        records = Entrez.read(handle)
        print(line[0:line.rfind('_') + 1] + records[0]['GBSeq_primary-accession'])
        sleep(1)
    else:
        print(line)
于 2016-09-14T10:19:05.063 に答える