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

Question

如果以前有人问过这个问题，我很抱歉，但是我已经搜索了几天，在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文件的大小超过1GB。我已经找到了解决这个问题的方法，通过基于重复ID删除重复项或使用bash，但遗憾的是我不能在我的计算机上这样做。 这项任务是研究项目，不是家庭作业或任务

提前感谢您的帮助

Answer 1

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

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 ]

Answer 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] )

Answer 3

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

使用Biopython，但与fasta文件一起使用，文件头为：

“>；标题类型请参见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的建议，使用列表的一组istead，速度更快：

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的，但是速度较慢，因为如果需要的话，序列计数可以发布它