使用 Python 删除 FASTA 中的重复序列答案

【问题标题】：Remove duplicated sequences in FASTA with Python使用 Python 删除 FASTA 中的重复序列
【发布时间】：2021-03-03 18:12:21
【问题描述】：

如果之前有人问过这个问题，我深表歉意，但我已经搜索了好几天，在 Python 中找不到解决方案。

我有一个大的 fasta 文件，包含标题和序列。

>cavPor3_rmsk_tRNA-Leu-TTA(m) range=chrM:2643-2717 5'pad=0 3'pad=0 strand=+ repeatMasking=none
GTTAAGGTGGCAGAGCCGGTAATTGCATAAAATTTAAGACTTTACTCTCA
GAGGTTCAACTCCTCTCCTTAACAC

>cavPor3_rmsk_tRNA-Gln-CAA_ range=chrM:3745-3815 5'pad=0 3'pad=0 strand=- repeatMasking=none
AGAGGGTCATAAAGGTTATGGGGTTGGCTTGAAACCAGCTTTAGGGGGTT
CAATTCCTTCCTCTCT

>cavPor3_rmsk_tRNA-Ser-TCA(m) range=chrM:6875-6940 5'pad=0 3'pad=0 strand=- repeatMasking=none
AGAGGGTCATAAAGGTTATGGGGTTGGCTTGAAACCAGCTTTAGGGGGTT
CAATTCCTTCCTCTCT

这是文件外观的一小部分。我只想保留第一个条目（标题和序列），如果你可以看到最后两个条目的序列是相同的。

输出如下所示：

>cavPor3_rmsk_tRNA-Leu-TTA(m) range=chrM:2643-2717 5'pad=0 3'pad=0 strand=+ repeatMasking=none
GTTAAGGTGGCAGAGCCGGTAATTGCATAAAATTTAAGACTTTACTCTCA
GAGGTTCAACTCCTCTCCTTAACAC

>cavPor3_rmsk_tRNA-Gln-CAA_ range=chrM:3745-3815 5'pad=0 3'pad=0 strand=- repeatMasking=none
AGAGGGTCATAAAGGTTATGGGGTTGGCTTGAAACCAGCTTTAGGGGGTT
CAATTCCTTCCTCTCT

问题在于 FASTA 文件的大小超过 1 GB。我找到了通过基于重复 ID 删除重复项或使用 bash 来解决此问题的方法，但遗憾的是我无法在我的计算机上执行此操作。此任务是针对研究项目，而不是家庭作业或任务。

提前感谢您的帮助！

【问题讨论】：

不应该FASTA文件如：> SEQUENCE_1 MTEITAAMVKELRESTGAGMMDCKNALSETNGDFDKAVQLLREKGLGKAAKKADRLAAEG LVSVKVSDDFTIAAMRPSYLSYEDLDMTFVENEYKALVAELEKENEERRRLKDPNKPEHK IPQFASRKQLSDAILKEAEEKIKEELKAQGKPEKIWDNIIPGKMNSFIADNSQLDSKLTL MGQFYVMDDKKTVEQVIAEKEKEFGGKIKIVEFICFEVGEGLEKKTEDFAAEVAAQL> SEQUENCE_2 SATVSEINSETDFVAKNDQFIALTKDTTAHIQSNSLQSVEELHSSTINGVKFEEYLKSQI ATIGENLVVRRFATLKAGANGVVNGYIHTNGRVGVVIAAACDSAEVASKSRDLLRQICMH 跨度>
bioinformatics.stackexchange.com/questions/13708/…
所以你应该读取第一个序列并附加到新文件，然后读取第二个序列检查新文件，如果不存在，将其附加到它，依此类推，直到第一个文件结束
biostars.org/p/250123 问题：使用 python 处理大型 fasta 文件的最佳方法。
biostars.org/p/710 : 问题：在 Python 中解析 Fasta 文件的正确方法

标签： python duplicates biopython fasta

【解决方案1】：

这是从这里复制的：Remove Redundant Sequences from FASTA file in Python

使用 Biopython，但适用于标题为以下内容的 fasta 文件：

'> header' 类型见FAsta Format Wiki

from Bio import SeqIO
import time

start = time.time() 

seen = []
records = []

for record in SeqIO.parse("INPUT-FILE", "fasta"):  
    if str(record.seq) not in seen:
        seen.append(str(record.seq))
        records.append(record)


#writing to a fasta file
SeqIO.write(records, "OUTPUT-FILE", "fasta")
end = time.time()

print(f"Run time is {(end- start)/60}")

按照 MattMDo 的建议，使用集合而不是列表更快：

seen = set()
records = []

for record in SeqIO.parse("b4r2.fasta", "fasta"):  
    if record.seq not in seen:
        seen.add(record.seq)
        records.append(record)

我有一个较长的使用 argparser 但速度较慢，因为如果需要计数序列可以发布它

【讨论】：

OP 的 FASTA 条目格式正确，只是由于 > 字符表示 Markdown 中的块引用而没有显示。我已经重新格式化了原来的帖子。
关于你的代码，对于大文件，如果seen 是一个集合而不是一个列表，它会更加更有效。
@MattDMo 它将如何应对 GB 大小的文件问题？？
这取决于机器有多少内存，但任何最近的机器都应该有至少 2-4GB 的 RAM，所以这里应该没有问题。不过，一般来说，如果数据结构的大小超过可用内存，您将获得MemoryError。
@MattDMo biostars.org/p/710：在 Python 中解析 Fasta 文件的正确方法。根据这个我的想法是可行的（不管它是否慢得要命）：如果内存不足，只需创建一次''fasta_sequences = SeqIO.parse(open(input_file),'fasta')''并附加序列1到循环中的新文件，每次重新创建输出文件的解析器，以检查输入文件的 n 序列是否要附加到输出文件。我错了吗？我在哪里可以得到一个 2GB 的 fasta 文件来试用它？大多数序列存储库即使对于blast比对结果也有大小限制

【解决方案2】：

如果行都具有相同的格式，因此在序列之前有 6 个空格分隔的字段，那么这很容易。您必须将所有唯一值存储在内存中。

memory = set()
for line in open('x.txt'):
    if len(line) < 5:
        continue
    parts = line.split(' ', 6)
    if parts[-1] not in memory:
        print( line.strip() )
        memory.add( parts[-1] )

【讨论】：

【解决方案3】：

如果您想在删除重复项的同时保留两个标题，您可以使用：

input1 = open("fasta.fasta")
dict_fasta = {record.id:record.seq for record in SeqIO.parse(input1,"fasta")}

tmp = {}
for key, value in dict_fasta.items():
  if value in tmp:
    tmp[value].append(key)
  else:
    tmp[value] = [ key ]

【讨论】：