根据带有间隙的id替换FASTA序列

Question

我在文件中有2000个fasta序列，如下所示：

 >T1        
AQSFDRATVFEARRAGYQRESRVALGKSTGVLEWHVFHVWAPRETTILVETLSQIECSGKGIADRRQENPLTI------ATAI--TSQLLELVDIVIMDFKAITQFFL     
>T484      
AQSFDRATVFEARRAGYQREARVALGKSTGKLEWQVFHVWAPRETTILVETLSQLECSGKGIADRRQENPLKI------ATAI--TSQLLELVDIVIMDFKAITQFFL     
>T582      
AQSFDRATVFEKRRAGYQREARVALGKSTGKLEWQVFHVWAPRETTILVETLSQLECSGKGIADRRQENPLKI------ATAI--TSQLLELVDIVIMDFKAITQFFL     
>T1424     
AQSFDRATVFEKRRAGYQLEARVALGKSTGKLEWQVFHVWAPRETTILVETLSQLECSGKGIANRRQENPLKI------ATAI--TSQLLELVDILIMDFKAITQFFL     
>T1552     
AQSFDRATIFEKRRAGYQIEARVALGKSTGKLEWQVYHAWAPRETTILVETLSQLENAGKGVANRRHENPLKI------ATAI--TSQLLELVDILMMDFKAITQFFL     
.
.
.

我希望从N个中随机抽取f个序列（N = 2000个序列）。

例如，如果f = 2，我从2000个序列中随机抽取2个序列。

f=2 

l=[] 

for i in range(f):
    x=randint(1, N)
    l.append(x)

在我的清单l中，例如，我将有[291，566]。然后，我要绘制第291和566序列：

> T1424
AQSFDRATVFEKRRAGYQLEARVALGKSTGKLEWQVFHVWAPRETTILVETLSQLECSGKGIANRRQENPLKI ------ ATAI - TSQLLELVDILIMDFKAITQFFL
> T1552
AQSFDRATIFEKRRAGYQIEARVALGKSTGKLEWQVYHAWAPRETTILVETLSQLENAGKGVANRRHENPLKI ------ ATAI - TSQLLELVDILMMDFKAITQFFL

我想用长度为55的间隙（“-”）替换这些序列：

> T1424
-----------------------------------------------------------------------------------------------------------------------------
> T1552
-----------------------------------------------------------------------------------------------------------------------------

我尝试使用此代码：

from Bio import SeqIO 
from random import *

records = list(SeqIO.parse("fichier1.txt", "fasta")) 
#print(records[0].id)  # first record
#print(records[0].seq)

N=len(records) #2000 

f=2 

l=[] 

for i in range(f):
    x=randint(1, N)
    l.append(x)

d={} 
for i2 in l:
    print(records[i2].id,records[i2].seq)
    d[str(records[i2].id)]=str(records[i2].seq)
    with open("fichier1.txt") as fichier, open("newfile.txt", "w") as newfile:
        texte = fichier.read()

        new_text += texte.replace(str(records[i2].seq), "---------------------------------------")
        newfile.write(new_text)

print(d)

没有用，因为有时文件中的序列可能相同，但标识符却不同。

根据ID和相应的序列，我想更改序列以引入间隔。

Answer 1

我可以解决您的编码问题（尽管我仍然不相信这是处理您的科学问题的最佳方法，即缺少某些核糖体蛋白）。无论如何，这是一个解决方案，而不是依靠序列的唯一性，我只是跟踪序列的索引：

import random

from Bio import SeqIO
from Bio.Seq import Seq

num_to_sample = 2

# WARNING: puts all sequences in memory
lst = list(SeqIO.parse("fichier1.txt", "fasta"))

sample = set(random.sample(range(len(lst)), num_to_sample))

for idx, record in enumerate(lst):
    if idx in sample:
        record.seq = Seq(len(record.seq) * "-")
    print(record.format("fasta"), end='')

示例输出：

>T1
AQSFDRATVFEARRAGYQRESRVALGKSTGVLEWHVFHVWAPRETTILVETLSQIECSGK
GIADRRQENPLTI------ATAI--TSQLLELVDIVIMDFKAITQFFL
>T484
AQSFDRATVFEARRAGYQREARVALGKSTGKLEWQVFHVWAPRETTILVETLSQLECSGK
GIADRRQENPLKI------ATAI--TSQLLELVDIVIMDFKAITQFFL
>T582
------------------------------------------------------------
------------------------------------------------
>T1424
------------------------------------------------------------
------------------------------------------------
>T1552
AQSFDRATIFEKRRAGYQIEARVALGKSTGKLEWQVYHAWAPRETTILVETLSQLENAGK
GVANRRHENPLKI------ATAI--TSQLLELVDILMMDFKAITQFFL