使用biopython的SeqIO模块发现独特的克隆

Question

我正在研究 DNA 的下一代测序 (NGS) 分析。我正在使用 SeqIO Biopython 模块来解析 Fasta 格式的 DNA 库。我只想过滤唯一的克隆（唯一的记录）。为此，我使用以下 python 代码。

seen=[]
unique_clones=[]
records=list(SeqIO.parse('DNA_library', 'fasta'))
for record in records:
  if str(record.seq) not in seen:
    seen.append(str(record.seq))
    unique_clones.append(record)
SeqIO.write(unique_clones, 'unique_clones.fasta', 'fasta')

一个 DNA 文库大约有 100 万条记录，而我有 100 多个文库需要分析。使用此代码过滤 1 个库需要 2 个多小时。这段代码提取唯一克隆的速度似乎非常慢。还有其他方法可以过滤唯一记录吗？

---

对于没有任何生物信息学工作经验的Python编码员

fasta格式的克隆有两个参数，ID（>id-number）和记录（ATCG），如下所示：

>id-No-1
ATCGGGCTAAATTCGACTGCAGT

>id-No-2
ATCGGGCTAAATTCGACTGCAGT

我只是想要根据记录过滤唯一克隆并想要打印唯一克隆（id 和记录）。

如果我的问题或解释不够清楚，请告诉我。

Answer 1

我没有你的文件，所以我无法测试你将获得的实际性能增益，但这里有一些对我来说很慢的事情：

• 行 records=list(SeqIO.parse('DNA_library', 'fasta')) 将记录转换为记录列表，这可能听起来无伤大雅，但如果您有数百万条记录，则会变得昂贵。根据 docs，SeqIO.parse(...) 返回一个迭代器，因此您可以直接迭代它。
• 跟踪查看的记录时，使用set而不是list。使用 in 执行成员资格检查时，列表必须迭代每个元素，而集合则在恒定时间内执行操作（更多信息此处）。

通过这些更改，您的代码将变为：

seen_records = set()
records_to_keep = []

for record in SeqIO.parse('DNA_library', 'fasta'):
  seq = str(record.seq)
  if seq not in seen_records:
    seen_records.add(seq)
    records_to_keep.append(record)

SeqIO.write(records_to_keep, 'unique_clones.fasta', 'fasta')

Answer 2

不需要像上面提到的那样使用 set，因为该代码中的 cpython 机器会检查这里是否存在 if seq not in saw_records: 以及这里的 seen_records。添加（seq）（sys.stdin）。所以我会尝试这个：

from Bio.SeqRecord import SeqRecord

seen_seqs = list()
IDs = list()

for record in SeqIO.parse('DNA_library', 'fasta'):
    if record.seq not in seen_seqs:
        seen_seqs.append(record.seq)
        IDs.append(record.id)

SeqIO.write((SeqRecord(seq, id=ID, name="", description="") for seq, ID in zip(seen_seqs, IDs)),
            'unique_clones.fasta',
            'fasta')

如果这对你来说仍然太慢，请发表评论，因为据我所知，这对 Cython 来说是有用的。当您在 python 中处理大数据时，建议在其 C/C++ 包和包中处理它们。不要使用 python 的循环和类似的东西。