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

Question

我是 python 编程新手，尝试解析 fasta 文件并计算属于文件中每个 ID 的读取次数(此处使用玩具示例，但计划用于宏基因组 seq 文件，每个文件有 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

您可以只拥有一个字典 id -> 出现次数并在迭代记录时更新它。另外，我认为你可以删除:

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

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

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]