如何使用生物技术提取和修剪fasta序列

Question

大家好，我是刚开始用生物技术做一项小任务的蟒蛇，我有两个文件--一个包含I列表和相关的number.eg

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 sequences.eg

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] =int(拆分的IDlength-1)
• 您的it非常相似，但不完全相同，唯一相同的部分是6个字符标识符，它可以用于映射两个文件(在假设该文件正常工作之前进行双重检查)。拥有相同的键使映射更加容易。 f1 =open(‘fasta _pruned.Quicka，'w') fasta= dict()，open("species.fasta"，"rU")作为spe:用于SeqIO.parse(spe，"fasta")：fasta[record.id.split(‘AC.26 ') 1]=记录长度t= dict()，其中open("seqid_len.txt")表示长度: split_IDlength = line.strip().split(’‘) lengthdict[split_IDlength.split('_')1] = int(split_IDlength1) for k，V在lengthdict.items()中:如果fasta.get(k)为None:继续fasta.get(‘>’+ k)打印(fastak.seqv：)f1.write(‘>{}\n’..format(K)) f1.write(str(fastak.seqv:) + '\n') f1.lose()

输出：

>F6ISE0 >F6HSF4 >F6PLK9 >F6HOT8