파이썬에서 fasta 헤더의 해당 GI 번호에서 NCBI의 Accession 번호 받기

Question

GenBank에서 GI 번호를 단계적으로 제거하고 다음 형식으로 헤더를 편집 한 위치에 여러 개의 fasta 파일이 저장된다는 경고가 계속 표시됩니다.

>SomeText_ginumber 

나는 다음과 같이 심지어 헤더 파일을이로 시작하지만 방법이 이상적으로 파이썬 함께, 내가 NCBI에서 각 GI에 해당하는 가입 번호를 얻을 수 있음을 출력 할 생각도 없어 :

>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에서 데이터를 가져 와서 이전 헤더를 인쇄하지만 액세스 번호와 나머지 입력 시퀀스를 완료하면 완료 :)

This 귀하의 두 가지 예를 위해 일하지만 더 많은 데이터 (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)