Samtools mpileup中读段名称对结果的影响机制解析

在基因组数据分析流程中，samtools mpileup是一个常用的工具，用于生成每个基因组位置的覆盖信息。近期有用户报告了一个有趣的现象：当使用不同工具（bcl2fastq和BCL Convert）对同一批测序数据进行解复用后，虽然生成的BAM文件内容几乎完全相同（仅读段名称不同），但samtools mpileup的输出结果却出现了差异。本文将深入解析这一现象背后的技术原理。

现象描述

用户对Illumina NovaSeq 6000测序数据分别使用bcl2fastq和BCL Convert进行解复用，随后使用相同版本的bwa和samblaster处理，生成两个BAM文件。这两个文件经samtools view逐条比对显示，除读段名称外所有记录完全一致。samtools stats和samtools flagstats的输出也完全相同。

然而，当使用samtools mpileup时，某些基因组位置的覆盖信息出现了差异。具体表现为：

1. 相同位置的覆盖碱基符号分布不同
2. 有时使用正向链读段，有时使用反向链读段
3. 在13146个不同位置中，有10个位置出现了这种现象

技术原理解析

这种现象源于samtools mpileup处理重叠读段对的特殊机制。当一对读段在基因组上存在重叠区域时，htslib会默认移除其中一个读段的覆盖信息，以避免重复计数。这一行为由以下因素决定：

1. 历史处理方式：早期版本单纯基于质量值决定保留哪个读段，但这可能导致严重的链特异性偏差（strand bias），特别是在扩增子测序等应用中。

2. 改进后的随机化策略：为避免系统性的链偏好性，当前版本采用了一种半随机化的选择方式。具体实现是通过对读段名称进行哈希计算（一种简单的位扰乱算法），从而决定保留哪个读段。这种设计确保了：

   • 相同输入会产生相同输出（可重复性）
   • 不会引入链特异性偏差
   • 处理结果与读段名称相关

3. 影响范围：这种现象仅影响存在重叠的读段对。在大多数测序协议中，读段对重叠的情况相对少见，因此对整体结果影响有限。

解决方案与建议

用户可以通过以下方式控制这一行为：

1. 忽略重叠移除：使用--ignore-overlaps-removal参数可以完全禁用重叠读段的移除功能。这样处理会：

   • 确保结果不受读段名称影响
   • 可能导致某些区域的覆盖度被重复计算
   • 适用于变异检测等对覆盖度精确性要求不高的场景

2. 评估应用场景：是否需要移除重叠读段取决于：

   • 数据类型：扩增子测序通常需要移除以避免系统偏差
   • 分析目的：表达谱分析通常需要移除，而变异检测可能不需要
   • 覆盖度水平：低覆盖度数据中1x和2x的差异更为关键

最佳实践建议

1. 保持流程一致性：在整个项目中统一使用同一种解复用工具，避免混用bcl2fastq和BCL Convert。

2. 明确分析需求：根据具体应用决定是否使用--ignore-overlaps-removal参数。

3. 结果验证：对于关键位点，建议检查重叠读段的处理方式是否影响结论。

4. 文档记录：在方法部分明确记录使用的参数和处理策略，确保结果可重复。

这一现象虽然看起来令人意外，但实际上体现了samtools团队为平衡技术精确性和生物学合理性所做的设计考量。理解这一机制有助于研究人员更好地解释分析结果，并做出适当的技术选择。