python - 用 biopython 解析 fasta 文件以计算属于每个 ID 的数字序列读取

Question

我是 python 编程的新手，正在尝试解析一个 fasta 文件并计算属于文件中每个 ID 的读取次数(在此处使用玩具示例，但计划用于每个具有 1-100,000 次读取的宏基因组序列文件主题)。我想获得的输出将是一个类似于以下内容的文本文件:

Total reads: 8
Mean read length: 232.5
Median: 234.5
Mode: 250
Max: 250
Min: 209

Sample   Count
001-00   1
002-00   4
003-00   3
Etc.

通过 biopython cookbook 和其他帖子中提供的示例，我已经能够拼凑出以下代码，这些代码将生成读取长度的描述性统计信息，并为我提供单次读取的 SampleID 和读取长度，但我似乎无法理解如何最好地计算每个 ID 在文件中出现的次数并按上述格式设置输出格式(用制表符分隔样本和制表符字段)。到目前为止我放在一起的代码是:

from Bio import SeqIO
import statistics

records = list(SeqIO.parse("test_fasta.fasta", "fasta"))
print("Total reads: %i" % len(records))

sizes = [len(rec) for rec in SeqIO.parse("test_fasta.fasta", "fasta")]

print("Mean read length:", statistics.mean(sizes))
print("Median:", statistics.median(sizes))
print("Mode:", statistics.mode(sizes))
print("Max:", max(sizes))
print("Min:", min(sizes))
print()

generator = SeqIO.parse("test_fasta.fasta", "fasta")
print("Sample", "\t", "Read Length")
for seqrecord in generator:
    idkeep, rest = seqrecord.id.split(';',1)
    print(idkeep, "\t", len(seqrecord))

但这当然给出了读取长度而不是计数。任何关于如何解决这个问题的想法都将不胜感激(对于我迄今为止所写内容中明显低效的任何想法也是如此)。

Answer 1

你可以只拥有一个 dict id -> number of occurrences 并在遍历记录时更新它。另外，我认为你可以删除:

generator = SeqIO.parse("test_fasta.fasta","fasta")

行，因为您已经有了 records 列表。您还可以更改行:

sizes = [len(rec) for rec in SeqIO.parse("test_fasta.fasta", "fasta")]

到:

sizes = [len(rec) for rec in records]

所以你只解析一次 fasta 文件。

我会添加到您的代码中:

occurrences_dict = {}
for seqrecord in records:
    idkeep, rest = seqrecord.id.split(';',1)
    occurrences_dict[idkeep] = occurrences_dict.get(idkeep, 0) + 1

然后您可以打印每个 ID 出现的次数:

for k in occurrences_dict:
    print k, occurrences_dict[k]