时间戳

首页 > 解决方案 > 根据其 id 用间隙替换 FASTA 序列

python - 根据其 id 用间隙替换 FASTA 序列

问题描述

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

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

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

例如如果 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和对应的sequence来看，我想改变sequence来引入gap。

标签： pythonbioinformaticsbiopythonfasta

解决方案

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

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

推荐阅读