用 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 食谱和其他帖子中提供的示例，我已经能够拼凑出以下代码，这些代码将生成读取长度的描述性统计信息，并为我提供 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")

行，因为您已经有了 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]

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

Parsing fasta file with biopython to count number sequence reads belonging to each ID

python

parsing

fasta

biopython