760. 线粒体DNA分析¶

核DNA vs mtDNA：
  核DNA：23对染色体，约30亿bp，双亲遗传
  mtDNA：1个环形分子，16569bp，只从母亲遗传

为什么mtDNA分析重要？
  1. 遗传病：线粒体疾病（MELAS、LHON等）
  2. 法医学：母系鉴定，降解样本（mtDNA多拷贝，不易降解）
  3. 群体遗传：追溯人类迁徙史（"线粒体夏娃"）
  4. 癌症研究：肿瘤中mtDNA突变和拷贝数变化
  5. 衰老研究：mtDNA突变累积与衰老

mtDNA的特殊之处：
  - 多拷贝：每个细胞100-10000个mtDNA分子
  - 异质性：同一细胞中可能有野生型和突变型mtDNA共存
  - 高突变率：比核DNA高10-17倍（无组蛋白保护）
  - NUMTs干扰：核DNA中有mtDNA片段的"化石"残留

# ===== 从WGS BAM中提取线粒体读段 =====

# Step 1: 提取比对到chrM的读段
samtools view -b \
  -h \                                  # 包含header
  sample.bam \                          # 输入全基因组BAM
  chrM \                                # 只提取chrM
  > sample_chrM.bam                     # 输出mtDNA BAM

samtools index sample_chrM.bam          # 建索引

# Step 2: 统计mtDNA测序深度
mosdepth \
  --by 100 \                            # 100bp窗口
  --no-per-base \                       # 不输出逐碱基深度
  sample_mito \                         # 输出前缀
  sample_chrM.bam                       # 输入BAM

# 查看平均深度
cat sample_mito.mosdepth.summary.txt
# chrM  16569  平均深度（通常是核DNA的100-1000倍）

# ===== GATK线粒体最佳实践流程 =====
# GATK推荐用Mutect2检测mtDNA变异（因为异质性类似体细胞突变）

# Step 1: 准备参考序列
# 使用移位后的参考（shifted reference），避免环形基因组的边界问题
# GATK提供了移位后的参考和区间文件

# Step 2: 提取mtDNA读段并重新比对
# 比对到标准chrM
gatk PrintReads \
  -R reference.fa \
  -I sample.bam \
  -L chrM \                             # 只取chrM区间
  -O sample_chrM_reads.bam              # 输出chrM reads

# 也比对到移位后的chrM（处理环形边界区域）
gatk PrintReads \
  -R reference_shifted.fa \             # 移位后的参考
  -I sample.bam \
  -L chrM_shifted \
  -O sample_chrM_shifted.bam

# Step 3: Mutect2检测变异（标准方向）
gatk Mutect2 \
  -R reference.fa \
  -I sample_chrM_reads.bam \            # mtDNA BAM
  -L chrM \                             # 线粒体区间
  --mitochondria-mode \                 # 关键：线粒体模式
  --max-reads-per-alignment-start 75 \  # 限制每位点最大读段数
  --max-mnp-distance 0 \               # 不合并MNP
  --annotation StrandBiasBySample \     # 链偏好注释
  -O sample_mito_raw.vcf.gz             # 输出原始变异

# Step 4: 过滤变异
gatk FilterMutectCalls \
  -R reference.fa \
  -V sample_mito_raw.vcf.gz \           # 原始变异
  --mitochondria-mode \                 # 线粒体模式
  --min-allele-fraction 0.01 \          # 最低等位基因频率1%（检测低异质性）
  -O sample_mito_filtered.vcf.gz        # 过滤后变异

# Step 5: 合并标准和移位方向的结果
# （移位方向捕获环形边界区域的变异，合并后得到完整结果）
gatk MergeVcfs \
  -I sample_mito_filtered.vcf.gz \
  -I sample_mito_shifted_filtered.vcf.gz \
  -O sample_mito_final.vcf.gz           # 最终mtDNA变异

# 查看结果
bcftools stats sample_mito_final.vcf.gz | head -30  # 基本统计

# ===== Python分析mtDNA异质性 =====
import pysam  # 导入pysam读取VCF
import pandas as pd  # 导入pandas
import matplotlib.pyplot as plt  # 导入绑定图

# 读取mtDNA VCF
vcf = pysam.VariantFile("sample_mito_final.vcf.gz")  # 打开VCF

variants = []  # 存储变异信息
for record in vcf:
    if record.filter.keys() == ["PASS"]:  # 只要PASS的变异
        for sample in record.samples:
            ad = record.samples[sample]["AD"]  # 等位基因深度
            if ad and sum(ad) > 0:
                ref_depth = ad[0]       # 参考等位基因深度
                alt_depth = sum(ad[1:]) # 替代等位基因深度
                total = ref_depth + alt_depth
                vaf = alt_depth / total if total > 0 else 0  # 变异等位基因频率

                variants.append({
                    "pos": record.pos,           # 位置
                    "ref": record.ref,            # 参考碱基
                    "alt": str(record.alts[0]),   # 变异碱基
                    "vaf": vaf,                   # VAF(=异质性水平)
                    "depth": total,               # 总深度
                    "gene": "unknown"             # 基因（后面注释）
                })

df = pd.DataFrame(variants)  # 转为DataFrame
print(f"检测到 {len(df)} 个mtDNA变异")

# 分类异质性水平
df["hetero_level"] = pd.cut(
    df["vaf"],
    bins=[0, 0.01, 0.1, 0.5, 0.99, 1.0],
    labels=["极低(<1%)", "低(1-10%)", "中(10-50%)", "高(50-99%)", "同质(>99%)"]
)
print(df["hetero_level"].value_counts())

# 可视化：异质性分布
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 左图：沿mtDNA基因组的变异分布
axes[0].scatter(df["pos"], df["vaf"], alpha=0.6, s=20)  # 散点图
axes[0].set_xlabel("mtDNA Position (bp)")  # x轴
axes[0].set_ylabel("VAF (Heteroplasmy Level)")  # y轴
axes[0].set_title("Variant Distribution along mtDNA")
axes[0].set_xlim(0, 16569)  # mtDNA全长

# 右图：VAF分布直方图
axes[1].hist(df["vaf"], bins=50, edgecolor="black", alpha=0.7)
axes[1].set_xlabel("VAF")
axes[1].set_ylabel("Count")
axes[1].set_title("Heteroplasmy Level Distribution")
axes[1].axvline(x=0.01, color="red", linestyle="--", label="1% threshold")
axes[1].legend()

plt.tight_layout()
plt.savefig("mtdna_heteroplasmy.png", dpi=300)
plt.show()

# ===== HaploGrep 3：单倍群分型 =====
# 安装：https://haplogrep.i-med.ac.at/ 或命令行版本
# pip install haplogrep  # 或下载jar文件

# 方法一：在线版（最简单）
# 上传VCF到 https://haplogrep.i-med.ac.at/
# 自动返回单倍群分类

# 方法二：命令行
java -jar haplogrep3.jar classify \
  --input sample_mito_final.vcf \       # 输入mtDNA VCF
  --output sample_haplogroup.txt \      # 输出单倍群结果
  --format vcf \                        # 输入格式
  --phylotree phylotree17 \             # 系统发育树版本
  --metric kulczynski                   # 距离度量

# 输出结果解释：
# SampleID  Haplogroup  Quality  Range
# sample1   H2a2a1      0.97     1-16569
# → 单倍群H2a2a1，质量分数0.97，覆盖全长

# 常见东亚单倍群：
# A, B, C, D, F, G, M, N, R
# 常见欧洲单倍群：
# H, U, K, J, T, V, W, X

# ===== 从WGS估算mtDNA拷贝数 =====

# 方法一：mosdepth计算深度比
# mtDNA-CN ≈ 2 × (mtDNA平均深度 / 核DNA平均深度)
# 乘以2是因为核DNA是二倍体

# 计算核DNA平均深度
mosdepth --no-per-base nuclear sample.bam
NUCLEAR_DEPTH=$(grep "total" sample.mosdepth.summary.txt | awk '{print $4}')

# 计算mtDNA平均深度
mosdepth --no-per-base mito sample_chrM.bam
MITO_DEPTH=$(grep "chrM" mito.mosdepth.summary.txt | awk '{print $4}')

# 估算mtDNA拷贝数
echo "scale=1; 2 * ${MITO_DEPTH} / ${NUCLEAR_DEPTH}" | bc
# 正常范围：100-1000（因组织类型不同）
# 心肌/肌肉: ~1000+
# 血液: ~100-500
# 肝脏: ~500-2000

# ===== Python批量计算mtDNA-CN =====
import subprocess  # 导入subprocess运行命令
import pandas as pd  # 导入pandas

samples = ["sample1", "sample2", "sample3"]  # 样本列表
results = []  # 存储结果

for sample in samples:
    # 用samtools计算深度
    # 核DNA深度（取chr1作为代表）
    cmd_nuclear = f"samtools depth -r chr1 {sample}.bam | awk '{{sum+=$3; n++}} END {{print sum/n}}'"
    nuclear_depth = float(subprocess.getoutput(cmd_nuclear))  # 核DNA深度

    # mtDNA深度
    cmd_mito = f"samtools depth -r chrM {sample}.bam | awk '{{sum+=$3; n++}} END {{print sum/n}}'"
    mito_depth = float(subprocess.getoutput(cmd_mito))  # mtDNA深度

    # 计算拷贝数
    mtdna_cn = 2 * mito_depth / nuclear_depth  # mtDNA-CN公式

    results.append({
        "sample": sample,
        "nuclear_depth": nuclear_depth,
        "mito_depth": mito_depth,
        "mtDNA_CN": mtdna_cn
    })

df = pd.DataFrame(results)
print(df.to_string(index=False))

# ===== NUMTs：核基因组中的mtDNA假象 =====
# NUMTs(Nuclear Mitochondrial DNA segments)是进化过程中
# 从线粒体转移到核基因组的mtDNA片段

# 问题：WGS中来自NUMTs的读段可能被错误比对到chrM
# 导致假的"异质性"或"变异"

# 解决方案一：使用GATK --mitochondria-mode
# 自动处理NUMTs干扰

# 解决方案二：过滤低VAF变异
# 真正的低异质性变异 vs NUMTs假象很难区分
# 一般建议：VAF < 1%的变异需谨慎解释

# 解决方案三：用专门的工具（如MitoHPC）
# MitoHPC先重新比对到环形化chrM，再回收NUMTs区域的误比对读段

# 解决方案四：检查变异是否在已知NUMTs区域
# 已知NUMTs位置可从UCSC Genome Browser下载

# ===== 线粒体DNA分析速查 =====

# 提取mtDNA读段
samtools view -b sample.bam chrM > mito.bam

# GATK Mutect2检测变异
gatk Mutect2 -R ref.fa -I mito.bam \
    -L chrM --mitochondria-mode \
    -O mito_raw.vcf.gz
gatk FilterMutectCalls -R ref.fa \
    -V mito_raw.vcf.gz --mitochondria-mode \
    --min-allele-fraction 0.01 \
    -O mito_filtered.vcf.gz

# 单倍群分型
java -jar haplogrep3.jar classify \
    --input mito.vcf --output haplo.txt

# mtDNA拷贝数
mtDNA-CN = 2 × (mtDNA深度 / 核DNA深度)

# 关键注意：
# chrM vs MT → 检查BAM的命名
# NUMTs → 用--mitochondria-mode
# 环形边界 → 用shifted reference
# 异质性 → VAF < 1%需谨慎