CoverM 的相对丰度为什么变了？用 100 条读段实跑覆盖深度、覆盖广度和未比对分母

同一份宏基因组数据，过滤前后平均覆盖深度相同，相对丰度却可能明显变化。遇到这种情况，先检查参与计算的读段集合、参考基因组长度和过滤条件，往往能够解释差异。尤其是提前删除未比对记录，或者把覆盖广度门槛设得更严格，都可能改变最终丰度表的含义。

本文安装 CoverM 0.8.0，构造一份可以手工核算的百条读段比对文件，实际运行三组覆盖统计，并另外检查读段计数与广度过滤的组合。实验直接使用设计好的比对记录，借助 samtools 转成排序后的 BAM，重点验证 CoverM 在这些输入条件下怎样计算和报告结果。

先分清深度、广度和相对丰度

平均覆盖深度描述参考序列每个位置平均被多少条读段覆盖。覆盖广度描述参考序列中有多少比例的位置至少得到一次覆盖。两者可以同时很高，也可能一个高、另一个低。

例如，一千条读段集中覆盖同一个短区域，会让这个区域堆积出很高的深度，但基因组其他位置仍然没有证据。若将所有位置一起平均，得到的平均深度可以大于一，覆盖广度却仍然很低。因此，判断某个基因组在样本中是否得到充分支持时，两个指标值得一起检查。

CoverM 可以按片段或基因组汇总多种覆盖指标。基因组模式的相对丰度还会考虑不同参考长度和实际得到支持的读段比例。使用这些指标时，需要同时记录计算方法和过滤参数。CoverM 官方项目说明

本次所有条件都把片段两端排除长度设为零。这样，参考的全部位置都参与平均深度计算，方便手工复算。正式项目是否排除片段两端，应结合比对方式、片段长度分布和末端覆盖情况决定，并保持不同样本之间的统计口径一致。

一百条读段如何构成这个实验

参考集合包含两个教学基因组。甲的长度为一千碱基，乙的长度为两千碱基。每条读段长度都是一百碱基，所有已比对读段均为无错配的完整匹配。

甲获得十条读段，分别覆盖连续的十个一百碱基区间，因此每个位置恰好覆盖一次。乙获得五十条读段，全部堆在同一个一百碱基区间。其余四十条读段以未比对记录的形式保存在文件中。这里使用可控输入，是为了把统计规则与比对算法的差异分开观察。

由输入设计可直接算出，甲的平均深度为一，覆盖广度为百分之百。乙的累计覆盖碱基数为五千，除以两千碱基得到平均深度二点五，但实际覆盖位置只有一百，覆盖广度为百分之五。

附带脚本自动创建参考序列、比对记录和基因组定义表，固定随机种子后可以重复生成相同输入。比对记录中的参考名称与定义表完全一致，所有读段标识唯一，并且没有次级或补充比对记录，便于集中观察本次关注的两个参数变化。

第一组保留全部记录，关闭覆盖广度过滤

基因组定义表使用两列，第一列为基因组名称，第二列为属于该基因组的参考片段名称。实际项目中，一个基因组可以对应多条片段，需要把全部对应关系列出来。整理定义表时，应检查片段是否遗漏，以及同一个片段是否意外重复归入多个基因组。

printf 'Genome_A\ta\nGenome_B\tb\n' > genomes.tsv

本次运行的核心命令如下。输入已经排序，命令直接读取 BAM 并统计平均深度、覆盖广度和相对丰度。参数值及程序版本可以在附带日志和结果记录中核对。本次 CoverM 发布版本

coverm genome \
  --bam-files all_reads.bam \
  --genome-definition genomes.tsv \
  --methods mean covered_fraction relative_abundance \
  --min-covered-fraction 0 \
  --contig-end-exclusion 0 \
  --threads 2 \
  --output-file all_no_breadth_filter.tsv

实际输出中，甲的平均深度为一，乙为二点五，广度分别为一和零点零五，与手算一致。甲的相对丰度约为百分之十七点一四，乙约为百分之四十二点八六，未比对行约为百分之四十。

理解这组丰度时，可以先按平均深度分配两者的相对份额。甲占一除以三点五，乙占二点五除以三点五，再乘以六十条已比对读段占全部一百条的比例。这就得到工具本次报告的两个百分比。

mapped_fraction = 60 / 100
a_percent = 1 / (1 + 2.5) * mapped_fraction * 100
b_percent = 2.5 / (1 + 2.5) * mapped_fraction * 100
print(a_percent, b_percent)

甲只有十条已比对读段，却得到约百分之十七点一四的丰度，原因就在于参考长度参与了覆盖计算。将基因组相对丰度与原始读段占比直接等同，会忽略这一步归一化。

第二组加入百分之十的覆盖广度门槛

第二次运行保持同一份 BAM，唯一改动是将最小覆盖广度设为百分之十。在本次版本的命令中使用 --min-covered-fraction 10，日志确认采用百分之十门槛。

乙只有百分之五的参考位置获得覆盖，因此被置为零。输出中的平均深度、覆盖广度和相对丰度都变成零。甲仍然保留，平均深度为一、广度为一，但相对丰度变成百分之十，未比对行变成百分之九十。

第二组的未比对行尤其值得认真解释。输入文件仍有六十条带有比对位置的记录，读段标志也保持原样。表格里报告的百分之九十，是结合基因组层面的过滤条件后得到的统计结果，包含了没有通过这一门槛的参考所对应的读段影响。

因此，如果想知道原始比对程序究竟有多少读段获得位置，应查看输入比对文件的标志统计。如果想知道在当前过滤条件下，还有多少读段未计入合格参考，则要结合 CoverM 的输出口径。报告这两个数字时，把它们各自对应的步骤写清楚，读者才容易理解差异。

这组对照也展示了广度门槛的用途。五十条读段集中在一个短区域，给出了相当高的平均深度，但支持范围有限。门槛能够帮助识别这种局部信号，至于它来自保守区域、重复序列还是其他原因，还需要检查真实比对位置和参考内容。

第三组提前删掉未比对记录

第三组先用 samtools 删除四十条未比对记录，得到只含六十条已比对记录的文件，再使用第一组的统计参数。甲和乙的平均深度及覆盖广度全部相同，丰度却分别升到约百分之二十八点五七和百分之七十一点四三。

samtools view -b -F 4 \
  -o mapped_only.bam all_reads.bam

从输入文件的角度看，分母已经只剩下已比对读段。对于 CoverM 能看到的这些记录，全部都有比对位置，所以未比对行报告为零，两个基因组的相对丰度加起来为百分之百。

如果不同批次中，有些文件保留未比对记录，有些文件已经做过类似过滤，再把丰度表直接合并，得到的差异可能同时包含文件处理方式的影响。开始多样本比较之前，可以先核对每份 BAM 的总记录数、已比对记录数和原始测序读段统计，明确丰度所使用的分母。

也应注意，本例采用单端、单一主比对记录。真实双端数据以及包含次级、补充比对的文件，需要依据实际过滤规则核对读段和记录的对应关系。把输入限定条件写入分析记录，可以帮助后续复查计数差异来自哪里。

一个实跑中遇到的参数组合问题

本次还尝试把读段计数与平均深度放在同一条命令里输出。当最小覆盖广度大于零时，本次版本明确报错，提示计数估计器不计算覆盖广度，因此不能与这个门槛组合使用。

处理方式是把原始读段计数另行计算，在计数命令中将最小覆盖广度设为零。实际得到甲十条、乙五十条，与输入设计完全一致。这样可以同时保留原始计数和经过广度条件筛选的丰度表，后续再按基因组名称连接。

coverm genome \
  --bam-files all_reads.bam \
  --genome-definition genomes.tsv \
  --methods count \
  --min-covered-fraction 0 \
  --output-file read_counts.tsv

将两张表合并时，建议给字段注明计算条件。某个基因组原始计数大于零、筛选后的平均深度却为零，在本次实验里完全可以解释。若只保留最后一张表，就会失去追踪局部支持信号的线索。

还可以用一个独立的逐位点统计来核对汇总值。本次补跑了下面的命令，要求把零覆盖位置也输出。甲得到一千个位点，深度总和为一千，所有位置均有覆盖。乙得到两千个位点，深度总和为五千，其中只有一百个位点得到覆盖，最高位点深度为五十。由这些逐位点数值重新求平均，可以再次得到前面的深度和广度。

samtools depth -aa all_reads.bam > per_base_depth.tsv

逐位点结果还能帮助判断覆盖是否集中。两个参考即使拥有相同的平均深度，也可能分别对应均匀覆盖与局部高峰。实际排查时，可以按窗口汇总深度，查看信号是否集中在少数保守基因、重复区域或片段末端，并与相关区域的序列注释对应起来。

展示结果时还应区分小数与百分数。覆盖广度列中的零点零五对应百分之五，相对丰度列已经以百分数报告。程序输出的未比对比例有时会显示为三十九点九九九九九六，这是浮点数表示造成的微小差异。表格展示时可以合理四舍五入，进行自动核对时则使用适当容差，同时保留原始输出以便追溯。

迁移到真实数据时，怎样逐项核对

首先，检查基因组定义和参考比对之间的名称对应关系。建索引后又重命名片段，或者在合并 FASTA 时遗漏了部分片段，都会影响统计对象。保存参考文件校验值和定义表，可以让丰度结果追溯到实际使用的参考集合。

其次，检查参考之间的近缘关系。多个高度相似的基因组可能竞争同一批读段，结果依赖比对和过滤策略。可以先分析参考集合的冗余程度，再检查多重比对、比对质量及支持位置的分布。一次汇总丰度计算无法独立解决所有近缘参考之间的归属问题。

再次，保持样本之间的参数和输入处理方式一致。除了本次测试的广度门槛，还应记录读段匹配条件、片段末端处理、是否保留未比对记录及实际使用的统计方法。比较结果时，先比较这些条件，再解释样本的生物学差异。

CoverM 的便利之处在于同一套参考可以输出多种可配合解释的统计量。相应地，结果也需要连同条件一起阅读。对于出现突变的样本，可以从丰度退回平均深度，再退回覆盖广度和具体位置，逐步定位变化。

本次脚本会生成完整输入，保存三组丰度统计、独立读段计数以及不兼容参数组合的错误信息。读者可以先手算这两个参考的结果，再将同样的检查顺序用于自己的多样本矩阵。下载实跑脚本与本次结果

当样本数量增加、需要统一参考和批量管理比对结果时，可了解云生信的生物信息学分析计算服务。