甜瓜重测序项目高重复率异常排查与定量归因报告

Posted on | Post modified | In

对 Picard 报告的高重复率进行归因。

问题背景：指标冲突与业务影响

异常现象

在甜瓜重测序样本（参考基因组：DHL92_v3.6.1）的标准分析流程中，我们观察到数据质控指标存在显著的不一致：

  • Fastp 统计（基于序列相似性）：样本重复率处于 2%~3% 的正常基线水平
  • Picard MarkDuplicates 统计（基于5’端比对坐标）：样本重复率畸高，普遍达到 12%~20% 以上
  • 换v4基因组后 picard统计的重复率降至8%左右

业务影响

该指标冲突直接影响对测序文库物理质量的判定，并可能导致大量有效数据在变异检测前被错误过滤。如不明确原因，可能：

  1. 误判为实验事故，触发不必要的数据重测
  2. 在变异检测中过度过滤，损失真实生物学信号
  3. 影响客户对数据质量的信任度

排查目标

对 Picard 报告的高重复率进行归因。

排查思路：从“指标冲突”到“空间分布”

物理本质分析：序列唯一性与坐标碰撞的对立

我们观察到的核心矛盾是两个工具对“重复”定义的偏差：

  • Fastp：对比的是 Read 序列本身。如果重复率低，说明文库中绝大多数分子的序列是独特的（Unique）。
  • Picard：对比的是 Read 比对到参考基因组后的 5’端坐标。

当时的思考点： 如果序列本身是独特的，但比对后的坐标却撞在了一起，这说明不同的物理分子被强行比对到了同一个基因组位置。

诊断启发式：全局随机误差 vs. 局部系统偏好

既然确定了是“坐标碰撞”而非“序列重复”，接下来的逻辑推导如下：

  • 全局性假说：如果这种碰撞是全基因组均匀分布的，那通常意味着文库复杂度极低（但这与 Fastp 结果矛盾）。
  • 局部性假说：如果这种碰撞集中在少数区域，那么它极有可能是某些局部因素导致的异常堆积（如 PCR 扩增偏好、组装塌缩，或是真实的拷贝数变异等）。

结论：因此，我们首先需要验证这种“碰撞”在空间上的分布情况。最直观的方法就是看这些重复 Reads 是不是在基因组的特定区域扎堆出现。

如果“坐标碰撞”是由上述局部异常引起的，那么这些区域通常会表现出两个互相关联的特征：

  1. Picard 重复率畸高（坐标重合）
  2. 覆盖深度（Depth）异常偏高（Reads 堆积）

因此，排查的第一阶段目标被简化为：寻找全基因组中是否存在对总重复率贡献巨大的“深度热点（Hotspots）”区域。

技术验证：基于空间异质性的证据链

第一阶段：宏观分布扫描（核心证据）

目的：定位对总体重复率贡献巨大的 Top 1% 深度热点窗口。

方法：1kb滑动窗口深度扫描，提取高密度区域。

代码实现（00_generate_dup_diagnosis_list.sh 核心逻辑）：

1
2
3
4
5
6
7
8
9
10
11
# 使用 mosdepth 计算全基因组和仅重复序列的深度
singularity exec ${MOSDEPTH_SIF} mosdepth -n -t ${THREADS} --fast-mode -b 1000 ${OUT_PREFIX}_all ${INPUT_BAM}
singularity exec ${MOSDEPTH_SIF} mosdepth -n -t ${THREADS} --fast-mode -b 1000 -i 1024 ${OUT_PREFIX}_dup ${INPUT_BAM}

# 合并结果，计算每个窗口的dup率
paste ${OUT_PREFIX}_all.regions.bed.gz ${OUT_PREFIX}_dup.regions.bed.gz | cut -f 1,2,3,4,8 > ${COMBINED_TSV}

# 提取Top 1%异常窗口
TOTAL_LINES=$(wc -l < ${COMBINED_TSV})
TOP_PCT=$((TOTAL_LINES / 100))
sort -k5,5nr ${COMBINED_TSV} | head -n $TOP_PCT > ${OUT_PREFIX}_hotspots.bed

关键发现（样本A10）：

  • Top 1%窗口贡献率：83.89%
  • 结论：重复率畸高为极端局部事件，非全局性建库异常。这是支持“坐标碰撞是由局部系统偏好引起”的最核心证据。

第二阶段：补充验证——TLEN 物理指纹抽样

在确定了空间分布异常后，为了完善证据链，我们引入了 TLEN（插入片段长度）作为辅助验证。需要承认的是，TLEN 分布很难独立看出来某个区域是什么原因造成的，因为往往是多个原因叠加导致了异常区域，但它可以作为抽样确认证据。

方法：提取异常区域的TLEN和MAPQ统计矩阵。

核心代码（get_region_stats.py 关键函数）：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
def calculate_stats(tlen_array, mapq_array):
    """计算TLEN和MAPQ的全面统计指标"""
    stats = {}
    
    # TLEN统计：重点关注IQR（四分位距）
    if len(tlen_array) > 0:
        q75, q25 = np.percentile(tlen_array, [75, 25])
        stats['tlen_iqr'] = q75 - q25  # 辅助诊断指标
        stats['tlen_unique_count'] = len(np.unique(tlen_array))
    
    # MAPQ统计：评估比对质量
    if len(mapq_array) > 0:
        stats['mapq_p0'] = np.sum(mapq_array == 0) / len(mapq_array)  # 无法比对比率
        stats['mapq_p30_plus'] = np.sum(mapq_array >= 30) / len(mapq_array)  # 高质比对比率
    
    return stats

物理证据矩阵（样本典型窗口抽样）：

染色体 物理区间 Dup率 文库TLEN IQR Dup Reads TLEN IQR 补充说明
chr00 33115000-33116000 78.00% 77.0 63.0 保持一定的打断离散度
chr12 4357000-4358000 83.75% 54.0 47.0 支持独立分子碰撞
chr00 15711000-15712000 81.19% 66.0 53.0 序列可能已分化但坐标重叠

辅助结论：大部分重复率畸高的区域，其 Duplicate Reads 仍保留了一定的插入片段长度离散度（$IQR > 40$），这作为辅助证据说明它们经历了相对独立的打断过程。

第三阶段：基因组版本交叉验证（最终确认）

验证设计：

  1. 测试组：使用高质量参考基因组v4.0重新比对
  2. 对照组：原始有缺陷基因组v3.6.1

结果对比：

  • v3.6.1基因组：Picard重复率 12-20%+
  • v4.0基因组：Picard重复率下降至 ~8%

正交验证：基因组组装完整度的提升直接导致重复率下降，形成完备的因果链条验证。

可视化诊断视图

基于排查经验，我开发了增强版可视化脚本（03_report_v2.py），通过生成综合的 HTML 报告，提供六个维度的诊断视图：

全样本 Fastp 与 Picard 重复率对比

目的：展示所有样本在 Fastp 和 Picard 中的全基因组 Dup 率差异对比。 诊断意义：评估基于比对物理坐标去重与基于序列去重算法的结果差异全局趋势。

单样本：各染色体平均 Duplication 率

目的：展示单样本每条主要染色体的平均 Dup 率。 诊断意义：辅助判断是否存在整条染色体级别的异常扩增或大片段组装异常。

单样本：异常热点真实 Dup 率分布堆叠图

目的：展示全基因组复制异常最严重的前 0.01%、0.1%、1% 窗口对应的真实总体重复率占比。 诊断意义：评估极端数据对全局重复率的贡献度，判断重复现象是高度集中于局部区域还是广泛存在。

单样本：Dup率沿基因组分布热图

目的：识别Dup率异常聚集的连续区域。 实现：100kb滑动窗口，颜色映射（纯白→深红）表示Dup率高低。 诊断意义：连续高 Dup 率的红色区域通常对应潜在的基因组组装缺陷区域。

单样本：深度分布热图

目的：识别绝对深度偏高的区域。 实现：100kb滑动窗口，展示绝对测序深度（最高截断为99%分位数）。 诊断意义：联合 Dup 率热图解读，若同一区域同时出现高深度与高 Dup 率，则提示可能存在真实的生物学多拷贝片段或组装塌缩。

单样本：深度-Dup率关系散点图

目的：展示全基因组各窗口深度与Dup率的关系分布情况。 诊断意义：揭示局部深度与重复率的经验分布规律，高深度往往客观伴随高碰撞概率。

使用方式：

1
2
3
4
python 03_report_v2.py \
  --sample A10 \
  --results-dir /path/to/results \
  --output-dir ./visualization_output

结论

简单来说，我们将问题锁定到了异常区域，但是目前尚不明确为什么会有这些异常区域。

后续方向

后续如果需要进一步探究原因，可以考虑以下几种可能性：

  1. PCR 偏好：某些特定序列（如高 GC 含量或特定结构区域）在文库构建的 PCR 扩增阶段更容易被扩增，导致这些区域在测序数据中出现大量重复堆积。
  2. 组装 Collapse（塌缩）：参考基因组的组装质量不佳，导致多个原本在基因组不同位置的相似或重复序列（如转座子）被错误地组装到了同一个坐标上。在比对时，来自多个位置的 Unique Reads 被迫堆积在这个单一的错误坐标点。
  3. 真实 CNV（拷贝数变异）：样本本身在某些区域发生了真实的生物学拷贝数扩增（例如某些抗性基因的天然多拷贝）。这使得该区域产生的 Reads 数量成倍增加，从而提高了比对时的坐标碰撞率。

扩展思考：重新建库改善重复率的机制分析

在参考基因组（组装缺陷）保持恒定的前提下，若通过重新建库降低了 Picard 重复率，其物理机制可能源于以下三种路径。我们采用贝叶斯分析框架，生成以下替代假说：

假说 A：文库复杂度（Library Complexity）的提升

机制：增加初始投入的 DNA 量或优化末端修复效率。 解释：Picard 判定重复的本质是坐标碰撞。在相同的测序深度下，文库中独特的物理分子越多（复杂度越高），由于随机打断恰好导致 5’ 端坐标完全一致的概率就会降低。这在统计上符合泊松分布过程，即分子多样性稀释了碰撞概率。

假说 B：PCR 扩增循环数的减少与偏好性控制

机制：减少 PCR 循环数或换用高保真、低偏好性的聚合酶。 解释：真实的 PCR 重复虽然在本次甜瓜样本中非主导因素，但在任何文库中都客观存在。优化 PCR 条件能直接降低生化层面的克隆比例。即使存在组装塌缩，总重复率也会因为 PCR 冗余的剥离而表现出数值上的下降。

假说 C：打断随机性（Fragmentation Randomness）的改善

机制：调整超声打断参数（如 Covaris 的处理时间或频率）。 解释：如果原有的打断过程存在“非随机偏好”（例如受局部 GC 含量或染色质结构影响，导致 DNA 更有可能在特定位点断裂），那么在比对时就会出现非生物意义的坐标聚集。改善打断的随机性可以使分子的 5’ 端在基因组上分布更均匀，从而避开坐标碰撞。

结论： 重新建库对 Dup 率的改善，通常是文库复杂度提升与扩增偏好降低共同作用的结果。对于存在组装缺陷的样本，提高文库复杂度是抵消“计算性碰撞假象”最有效的实验手段。这一章的加入使文档从单一的“故障排查”升级为对“实验与分析交互”的深入理解。

附录：关键代码与使用说明

诊断工具集

  • 宏观分布扫描：00_generate_dup_diagnosis_list.sh
  • 微观指纹提取：get_region_stats.py
  • 可视化报告：03_report_v2.py（增强版）

快速使用指南

1
2
3
4
5
6
7
8
# 1. 运行宏观扫描（SLURM集群）
bash 00_generate_dup_diagnosis_list.sh

# 2. 对单个样本进行详细分析
python get_region_stats.py -b sample.bam -r hotspots.bed -o sample_stats.tsv

# 3. 生成综合可视化报告
python 03_report_v2.py --sample A10 --results-dir ./results --output-dir ./report

依赖环境

  • Python 3.8+ (numpy, matplotlib, pysam)
  • Singularity (运行mosdepth)
  • SLURM集群（可选，支持分布式计算）