如何使用 biopython 提取和修剪 fasta 序列

Question

大家好，我是 python 新手，正在努力使用 biopython 完成一项小任务。我有两个文件 - 一个包含 id 列表和相关编号。例如

id.txt

tr_F6LMO6_F6LMO6_9LE  25
tr_F6ISE0_F6ISE0_9LE  17
tr_F6HSF4_F6HSF4_9LE  27
tr_F6PLK9_F6PLK9_9LE  19
tr_F6HOT8_F6HOT8_9LE  29

第二个文件包含一个大的 fasta 序列。例如下面

fasta_db.fasta

    >tr|F6LMO6|F6LMO6_9LEHG Transporter
    MLAPETRRKRLFSLIFLCTILTTRDLLSVGIFQPSHNARYGGMGGTNLAIGGSPMDIGTN
    PANLGLSSKKELEFGVSLPYIRSVYTDKLQDPDPNLAYTNSQNYNVLAPLPYIAIRIPIT
    EKLTYGGGVYVPGGGNGNVSELNRATPNGQTFQNWSGLNISGPIGDSRRIKESYSSTFYV

   >tr|F6ISE0|F6ISE0_9LEHG peptidase domain protein OMat str.  
    MPILKVAFVSFVLLVFSLPSFAEEKTDFDGVRKAVVQIKVYSQAINPYSPWTTDGVRASS
    GTGFLIGKKRILTNAHVVSNAKFIQVQRYNQTEWYRVKILFIAHDCDLAILEAEDGQFYK

   >tr|F6HSF4|F6HSF4_9LEHG hypothetical protein,  
    MNLRSYIREIQVGLLCILVFLMSLYLLYFESKSRGASVKEILGNVSFRYKTAQRKFPDRM
    LWEDLEQGMSVFDKDSVRTDEASEAVVHLNSGTQIELDPQSMVVLQLKENREILHLGEGS


   >tr|F6PLK9|F6PLK9_9LEHG Uncharacterized protein mano str. 
   MRKITGSYSKISLLTLLFLIGFTVLQSETNSFSLSSFTLRDLRLQKSESGNNFIELSPRD
   RKQGGELFFDFEEDEASNLQDKTGGYRVLSSSYLVDSAQAHTGKRSARFAGKRSGIKISG

我想将第一个文件中的 id 与第二个文件匹配，并在删除长度（从 1 到 25，在 eq 中）后在新文件中打印那些匹配的 seq。

例如输出[25（与id关联的值，第一个文件），当id匹配时，aa从开始删除]。

fasta_pruned.fasta

>tr|F6LMO6|F6LMO6_9LEHG Transporter     
LLSVGIFQPSHNARYGGMGGTNLAIGGSPMDIGTNPANLGLSSKKELEFGVSL
PYIRSVYTDKLQDPDPNLAYTNSQNYNVLAPLPYIAIRIPITEKLTYGGGVYV
PGGGNGNVSELNRATPNGQTFQNWSGLNISGPIGDSRRIKESYSSTFYV

Biopython 食谱在我刚接触 python 编程时远远超出我的想象。感谢您提供的任何帮助。

我试过了，搞砸了。就是这样。

from Bio import SeqIO
from Bio import Seq

f1 = open('fasta_pruned.fasta','w')


lengthdict = dict() 
with open("seqid_len.txt") as seqlengths:
    for line in seqlengths:
        split_IDlength  = line.strip().split(' ')
        lengthdict[split_IDlength[0]] = split_IDlength[1]


with open("species.fasta","rU") as spe:
    for record in SeqIO.parse(spe,"fasta"):
        if record[0] == '>' :
            split_header = line.split('|')
            accession_ID = split_header[1]
            if accession_ID in lengthdict:
                f1.write(str(seq_record.id) + "\n")
                f1.write(str(seq_record_seq[split_IDlength[1]-1:]))



f1.close()

Answer 1

您的代码几乎包含所有内容，除了一些妨碍它提供所需输出的小东西：

• 您的文件id.txt 在id 和位置之间有两个空格。你取第二个元素，在这种情况下它是空的。

• 当文件被读取时，它被解释为一个字符串，但你希望位置是一个整数

  lengthdict[split_IDlength[0]] = int(split_IDlength[-1])
  

• 您的 id 非常相似但不相同，唯一相同的部分是 6 个字符的标识符，可用于映射两个文件（在您认为它有效之前请仔细检查）。拥有相同的键使映射更容易。

  f1 = open('fasta_pruned.fasta', 'w')
  
  fasta = dict()
  with open("species.fasta","rU") as spe:
      for record in SeqIO.parse(spe, "fasta"):
          fasta[record.id.split('|')[1]] = record
  
  lengthdict = dict() 
  with open("seqid_len.txt") as seqlengths:
      for line in seqlengths:
          split_IDlength  = line.strip().split(' ')
          lengthdict[split_IDlength[0].split('_')[1]] = int(split_IDlength[1])
  
  for k, v in lengthdict.items():
      if fasta.get(k) is None:
          continue
      print('>' + k)
      print(fasta[k].seq[v:])
      f1.write('>{}\n'.format(k))
      f1.write(str(fasta[k].seq[v:]) + '\n')
  
  f1.close()
  

输出：

>F6LMO6
LLSVGIFQPSHNARYGGMGGTNLAIGGSPMDIGTNPANLGLSSKKELEFGVSLPYIRSVYTDKLQDPDPNLAYTNSQNYNVLAPLPYIAIRIPITEKLTYGGGVYVPGGGNGNVSELNRATPNGQTFQNWSGLNISGPIGDSRRIKESYSSTFYV
>F6ISE0
LPSFAEEKTDFDGVRKAVVQIKVYSQAINPYSPWTTDGVRASSGTGFLIGKKRILTNAHVVSNAKFIQVQRYNQTEWYRVKILFIAHDCDLAILEAEDGQFYK
>F6HSF4
YFESKSRGASVKEILGNVSFRYKTAQRKFPDRMLWEDLEQGMSVFDKDSVRTDEASEAVVHLNSGTQIELDPQSMVVLQLKENREILHLGEGS
>F6PLK9
IGFTVLQSETNSFSLSSFTLRDLRLQKSESGNNFIELSPRDRKQGGELFFDFEEDEASNLQDKTGGYRVLSSSYLVDSAQAHTGKRSARFAGKRSGIKISG
>F6HOT8