液体活检 ctDNA 分析¶

血液采样 → cfDNA提取 → 文库构建(含UMI) → 超深度测序 → UMI去重 → 低频突变检测 → 临床解读

# UMI（Unique Molecular Identifiers）是液体活检的核心技术
# 每个原始DNA分子被贴上唯一标签，去除PCR重复后保留真实突变

# 安装 fgbio（UMI处理的标准工具）
conda install -c bioconda fgbio  # 安装fgbio

# 第 1 步：提取 UMI
fgbio ExtractUmisFromBam \
    --input raw.bam \               # 原始BAM（含UMI序列）
    --read-structure 3M2S+T \       # 读结构：3bp UMI + 2bp skip + 模板
    --output extracted.bam           # 输出含UMI标记的BAM

# 第 2 步：比对
# 将提取UMI后的reads比对到参考基因组
samtools fastq extracted.bam | \     # BAM转FASTQ
    bwa mem -t 8 reference.fa - | \  # 比对
    samtools sort -o aligned.bam     # 排序

# 第 3 步：分组（Group Reads by UMI）
fgbio GroupReadsByUmi \
    --input aligned.bam \            # 比对后的BAM
    --strategy adjacency \           # UMI分组策略（允许1bp错误）
    --edits 1 \                      # 允许的UMI编辑距离
    --output grouped.bam             # 分组后的BAM

# 第 4 步：构建共识序列（Consensus）
fgbio CallMolecularConsensusReads \
    --input grouped.bam \            # 分组后的BAM
    --min-reads 3 \                  # 至少3条reads形成共识（重要！）
    --output consensus.bam           # 共识序列BAM

# 第 5 步：重新比对共识序列
samtools fastq consensus.bam | \
    bwa mem -t 8 reference.fa - | \
    samtools sort -o consensus_aligned.bam

# 第 6 步：过滤
fgbio FilterConsensusReads \
    --input consensus_aligned.bam \
    --ref reference.fa \
    --min-reads 3 \                  # 最少支持reads
    --min-base-quality 30 \          # 最低碱基质量
    --max-base-error-rate 0.1 \      # 最大碱基错误率
    --output filtered_consensus.bam  # 过滤后的共识BAM

# 标准变异检测工具在低VAF时假阳性太高
# 需要专门的低频突变检测工具

# 方法1：VarDict（支持低至0.1% VAF）
conda install -c bioconda vardict-java  # 安装VarDict

vardict-java \
    -G reference.fa \                # 参考基因组
    -f 0.001 \                       # 最低VAF阈值0.1%
    -b filtered_consensus.bam \      # 输入BAM
    -R chr7:55019017-55211628 \      # 目标区域（如EGFR）
    -c 1 -S 2 -E 3 -g 4 \           # BED文件列索引
    | teststrandbias.R | \           # 链偏好检验
    var2vcf_valid.pl \
    -N sample_name \
    -f 0.001 \                       # VAF阈值
    > variants.vcf                   # 输出VCF

# 方法2：Mutect2 mitochondrial mode（适合极低VAF）
gatk Mutect2 \
    -R reference.fa \
    -I filtered_consensus.bam \
    -L targets.bed \
    --min-base-quality-score 30 \
    --f1r2-median-mq 30 \
    --min-allele-fraction 0.001 \    # 最低VAF 0.1%
    -O ctdna_variants.vcf.gz

# 方法3：专用ctDNA工具 — CAPP-Seq分析用 iDES
# 集成数字错误抑制（integrated Digital Error Suppression）

# ========== Python脚本：MRD追踪分析 ==========
import pandas as pd
import matplotlib.pyplot as plt

# 读取多个时间点的突变VAF
timepoints = {
    "基线": {"TP53_R175H": 5.2, "KRAS_G12D": 3.8, "EGFR_T790M": 0.0},
    "化疗2周期": {"TP53_R175H": 1.1, "KRAS_G12D": 0.8, "EGFR_T790M": 0.0},
    "化疗4周期": {"TP53_R175H": 0.05, "KRAS_G12D": 0.02, "EGFR_T790M": 0.0},
    "随访3月": {"TP53_R175H": 0.01, "KRAS_G12D": 0.0, "EGFR_T790M": 0.3},
    "随访6月": {"TP53_R175H": 0.8, "KRAS_G12D": 0.1, "EGFR_T790M": 2.1},
}

# 转换为DataFrame
df = pd.DataFrame(timepoints).T                     # 转置
df.index.name = "Timepoint"                           # 设置索引名

# 画VAF变化曲线
plt.figure(figsize=(12, 6))
for mutation in df.columns:
    plt.plot(range(len(df)), df[mutation],             # 画每个突变的VAF变化
             marker='o', label=mutation, linewidth=2)

plt.xticks(range(len(df)), df.index, rotation=45)     # X轴标签
plt.ylabel("VAF (%)")                                  # Y轴标签
plt.xlabel("Time Point")                                # X轴标签
plt.title("ctDNA Monitoring - MRD Tracking")           # 标题
plt.legend()
plt.axhline(y=0.1, color='red', linestyle='--',        # 检测限
            label='Detection limit')
plt.yscale('log')                                       # Y轴对数刻度
plt.tight_layout()
plt.savefig("mrd_tracking.pdf")                         # 保存图片

# 利用cfDNA片段大小模式检测癌症
# 肿瘤来源的ctDNA片段通常比正常cfDNA更短

import pysam
import numpy as np

def fragment_size_distribution(bam_file):
    """计算cfDNA片段大小分布"""
    bam = pysam.AlignmentFile(bam_file, "rb")
    sizes = []

    for read in bam.fetch():
        if read.is_proper_pair and not read.is_secondary:
            frag_size = abs(read.template_length)      # 片段大小
            if 50 < frag_size < 500:                   # 合理范围
                sizes.append(frag_size)

    sizes = np.array(sizes)

    # 计算关键指标
    short_frag = np.sum((sizes >= 100) & (sizes <= 150))  # 短片段(100-150bp)
    long_frag = np.sum((sizes >= 151) & (sizes <= 250))   # 长片段(151-250bp)
    ratio = short_frag / long_frag if long_frag > 0 else 0

    print(f"中位片段大小: {np.median(sizes):.0f} bp")
    print(f"短/长片段比: {ratio:.3f}")
    print(f"短片段比例高 → 可能有肿瘤信号")

    return sizes

sizes = fragment_size_distribution("cfDNA_sample.bam")

# === ctDNA 分析速查 ===

# UMI处理（fgbio）
fgbio ExtractUmisFromBam --input raw.bam --read-structure 3M2S+T --output umi.bam
fgbio GroupReadsByUmi --input aligned.bam --strategy adjacency --output grouped.bam
fgbio CallMolecularConsensusReads --input grouped.bam --min-reads 3 --output consensus.bam

# 低频突变检测
vardict-java -G ref.fa -f 0.001 -b consensus.bam -R region > variants.vcf

# 关键参数：
# UMI最小支持reads: ≥3
# VAF检测下限: ~0.1% (UMI) 或 ~1% (无UMI)
# 推荐测序深度: >10,000X (原始) → ~2,000X (去重后)
# cfDNA正常片段大小峰: ~167bp (1个核小体)