Pythonのfastaヘッダーの対応するGI番号からNCBIのアクセッション番号を取得

Question

GenBankでGI番号を段階的に削除していて、次の形式でヘッダーを編集した場所にいくつかのfastaファイルが保存されています。

>SomeText_ginumber 

次のように私はNCBIおよび出力から各GIのためのヘッダーを持つファイルを、対応するアクセッション番号を取得することができることを、理想的にはPythonで、私もこれで始めるには考えてきませんが、方法があります：

>SomeText_accessionnumber 

別のexaファイル形式のmple：

>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) 

Answer 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)