古DNA分析¶

一句话概述¶

古DNA（ancient DNA, aDNA）分析的生信流程——包括损伤模式识别、mapDamage认证、PMDtools提取真实古代reads、污染评估和群体遗传学推断，用于研究人类迁徙和进化历史。

核心知识点总览¶

知识点	关键内容	重要程度
aDNA损伤模式	脱氨基C→T/G→A末端富集	⭐⭐⭐⭐⭐
mapDamage2	损伤模式定量与认证	⭐⭐⭐⭐⭐
PMDtools	基于PMD score筛选古代reads	⭐⭐⭐⭐
污染评估	mtDNA/X染色体污染检测	⭐⭐⭐⭐
数据预处理	接头去除/合并overlap/短reads处理	⭐⭐⭐⭐
认证标准	建立aDNA数据可信度的多重证据	⭐⭐⭐⭐
群体遗传推断	PCA/ADMIXTURE/f-statistics	⭐⭐⭐
分子定年	放射性碳定年/tip-dating	⭐⭐⭐

各步骤详解¶

第一步：古DNA的特殊性质¶

白话解释： 古DNA来自数百到数万年前的生物遗骸（如骨骼、牙齿、毛发）。与现代DNA相比，aDNA有几个独特特征：(1) 极度断裂——片段通常只有30-80bp；(2) 化学损伤——最典型的是胞嘧啶脱氨基变成尿嘧啶（被测序仪读作T）；(3) 含量极低——可能只占总DNA的0.01-5%（其余是微生物DNA）。这些特性决定了aDNA分析需要特殊的处理策略。

技术细节： aDNA损伤特征（damage patterns）： - 脱氨基（deamination）：C→U(T)，5mC→T，在片段末端富集 - 碱基丢失（depurination）：purine位点断裂，导致片段短小 - 交联（cross-links）：分子间交联阻止PCR扩增 - 末端修复（end repair）：片段末端常有单链悬挂区

损伤的诊断特征： - 5'端：C→T替换频率升高（正链） - 3'端：G→A替换频率升高（互补链上的C→T） - 频率从末端向内部衰减

# aDNA损伤模式示意
# 5' ----[C→T enriched]----[normal]----[G→A enriched]---- 3'
#         ↑ position 1-3                ↑ last 3 positions
#
# 典型损伤频率：
# Position 1: C→T ~30-50% (古老样本)
# Position 2: C→T ~10-15%
# Position 3: C→T ~5-8%
# 内部: C→T ~1-2% (background)

第二步：数据预处理¶

白话解释： 古DNA测序通常使用特殊的建库方法（如单链文库、UDG处理/不处理）。由于片段很短，双端测序的reads常常有大量overlap。预处理需要去接头、合并overlap的paired-end reads，并特别注意保留短reads（不能用常规的length filter）。

技术细节：

# === aDNA数据预处理 ===

# 1. 接头去除与reads合并（AdapterRemoval2）
# aDNA reads短，PE reads大量overlap需要合并
AdapterRemoval \
    --file1 sample_R1.fq.gz \
    --file2 sample_R2.fq.gz \
    --basename sample \
    --trimns \
    --trimqualities \
    --minquality 20 \
    --minlength 30 \
    --collapse \
    --threads 8

# 输出文件：
# sample.collapsed        - 合并的reads（重要！aDNA分析主要用这个）
# sample.collapsed.truncated - 质量修剪后的合并reads
# sample.pair1.truncated  - 未合并的R1
# sample.pair2.truncated  - 未合并的R2

# 也可用fastp处理
fastp -i sample_R1.fq.gz -I sample_R2.fq.gz \
    --merge --merged_out sample.merged.fq.gz \
    --overlap_len_require 11 \
    --length_required 30 \
    -h sample_fastp.html

# 2. 比对（使用BWA aln 适合短reads）
# BWA aln比BWA-MEM更适合短aDNA片段（<70bp）
bwa aln -n 0.01 -o 2 -l 1024 -t 16 \
    reference.fasta sample.collapsed.gz > sample.sai

bwa samse reference.fasta sample.sai sample.collapsed.gz | \
    samtools sort -@ 8 -o sample.sorted.bam

# 对于较长片段（>70bp），BWA-MEM也可以
bwa mem -t 16 reference.fasta sample.collapsed.gz | \
    samtools sort -@ 8 -o sample.sorted.bam

# 3. 去重复（aDNA中PCR重复率通常很高）
# 使用picard（或sambamba/samtools markdup）
picard MarkDuplicates \
    INPUT=sample.sorted.bam \
    OUTPUT=sample.dedup.bam \
    METRICS_FILE=sample.dup_metrics.txt \
    REMOVE_DUPLICATES=true

# 4. 过滤比对质量
samtools view -b -q 30 -F 4 sample.dedup.bam > sample.filtered.bam
samtools index sample.filtered.bam

# 统计
samtools flagstat sample.filtered.bam
# 关注：总reads数、mapped%、重复率
# aDNA典型情况：内源性DNA占0.1-5%，重复率高

第三步：mapDamage2损伤认证¶

白话解释： mapDamage2是验证你的数据确实来自古DNA（而非现代污染）的关键工具。它统计reads末端的碱基替换频率，生成损伤模式图。如果看到典型的5' C→T和3' G→A末端富集模式，就能确认这确实是古代DNA。

技术细节： mapDamage2不仅生成损伤图谱，还能进行统计建模——拟合损伤参数（脱氨基概率、过悬概率），并可选择性地重新校正碱基质量（将可能的损伤位点质量下调）。

# === mapDamage2 运行 ===

# 基本模式：生成损伤图谱
mapDamage -i sample.filtered.bam \
    -r reference.fasta \
    -d mapDamage_output \
    --no-stats  # 只画图不做统计模型

# 完整模式：包含贝叶斯统计模型
mapDamage -i sample.filtered.bam \
    -r reference.fasta \
    -d mapDamage_output_full \
    --merge-reference-sequences  # 合并所有参考序列统计（布尔标记，无需数值参数）

# 带碱基质量rescaling的模式
mapDamage -i sample.filtered.bam \
    -r reference.fasta \
    -d mapDamage_rescaled \
    --rescale

# 输出文件说明：
# Fragmisincorporation_plot.pdf - 损伤模式图（最重要！）
# Length_plot.pdf               - 片段长度分布
# 5pCtoT_freq.txt              - 5'端C→T频率
# 3pGtoA_freq.txt              - 3'端G→A频率
# dnacomp.txt                  - 核苷酸组成
# Stats_out_MCMC_*             - 贝叶斯参数估计

# 解读关键图形：
# - 5'端C→T频率：真实aDNA应该position 1有明显升高（>5%）
# - 3'端G→A频率：对称的末端升高
# - 片段长度：应该呈现短片段为主的分布（peak在30-60bp）

# 编程方式解析mapDamage结果
import pandas as pd
import matplotlib.pyplot as plt

# 读取5'端C→T频率
ct_5prime = pd.read_csv("mapDamage_output/5pCtoT_freq.txt", sep="\t")
# 读取3'端G→A频率
ga_3prime = pd.read_csv("mapDamage_output/3pGtoA_freq.txt", sep="\t")

# 判断是否为真实aDNA
pos1_ct = ct_5prime.iloc[0]['C>T']  # 第1个位置的C→T频率
print(f"5' position 1 C>T frequency: {pos1_ct:.4f}")
if pos1_ct > 0.05:
    print("✓ Consistent with ancient DNA damage pattern")
elif pos1_ct > 0.01:
    print("△ Weak damage signal - may be aDNA or partially UDG-treated")
else:
    print("✗ No damage signal - likely modern DNA or fully UDG-treated")

第四步：PMDtools筛选古代reads¶

白话解释： 即使确认样本含有古DNA，其中也混有现代微生物DNA和可能的人类污染DNA。PMDtools（PostMortem Damage tools）对每条read计算一个"PMD score"——反映该read携带古DNA损伤信号的程度。通过设置阈值，只保留PMD score高的reads，可以富集真正来自古代的DNA片段。

技术细节： PMD score基于read中C→T（5'端）和G→A（3'端）替换的位置和频率计算似然比。高PMD score的reads更可能来自古代，低分的可能是现代污染。

# === PMDtools 使用 ===

# 安装
pip install pmdtools
# 或从GitHub: https://github.com/pontussk/PMDtools

# 1. 计算PMD scores
samtools view -h sample.filtered.bam | \
    pmdtools --threshold 3 --header | \
    samtools view -Sb - > sample.pmd_filtered.bam

# --threshold 3: 只保留PMD score >= 3的reads
# 阈值选择：
# 0: 保留所有reads
# 1-2: 温和过滤
# 3: 标准阈值（推荐）
# 5+: 严格过滤（可能损失过多数据）

# 2. 只计算PMD scores（不过滤）
samtools view sample.filtered.bam | \
    pmdtools --platypus > pmd_scores.txt

# 3. 查看PMD score分布
# 真实aDNA样本应有明显的高分尾部
awk '{print $NF}' pmd_scores.txt | sort -n | uniq -c | \
    awk '{print $2"\t"$1}' > pmd_score_distribution.txt

# 4. 比较过滤前后的数据量
echo "Before PMD filter:"
samtools view -c sample.filtered.bam
echo "After PMD filter (threshold=3):"
samtools view -c sample.pmd_filtered.bam

第五步：污染评估¶

白话解释： 古DNA实验室最头疼的问题就是污染——来自实验人员、其他样本或环境的现代DNA。污染评估是aDNA分析的必经步骤。主要通过线粒体DNA和X染色体（男性样本）的异质性来估计污染水平。

技术细节： 污染评估方法： 1. mtDNA contamination：线粒体是单倍体，真实样本应只有一个线粒体谱系。如果检测到多个线粒体变体，说明有污染。 2. X chromosome contamination（男性）：男性只有1条X，应为纯合。X上的杂合位点提示污染。 3. ANGSD contamMix：基于等位基因频率估计污染比例。

# === 污染评估 ===

# 方法1：schmutzi - 线粒体污染估计
# 首先提取线粒体reads
samtools view -b sample.filtered.bam chrM > sample.mt.bam
samtools index sample.mt.bam

# 运行schmutzi（注意：BAM需要先加MD tags）
samtools calmd -b sample.mt.bam reference_mt.fa > sample.mt.MD.bam
samtools index sample.mt.MD.bam

# contDeam估计污染（--library double用于双链文库，single用于单链文库）
contDeam --library double --out schmutzi_out/sample \
    --uselength --ref reference_mt.fa sample.mt.MD.bam

# schmutzi联合估计（需要contaminant数据库）
schmutzi --ref reference_mt.fa \
    --out schmutzi_result/sample \
    schmutzi_out/sample sample.mt.MD.bam

# 输出: schmutzi_result_final.cont.est 包含污染比例估计

# 方法2：ANGSD contamination estimate (X chromosome, males)
angsd -i sample.filtered.bam \
    -r chrX:5000000-154900000 \
    -doCounts 1 -iCounts 1 \
    -minMapQ 30 -minQ 20 \
    -out angsd_xcont

# 使用hapmap数据估计污染
Rscript angsd_contMix.R angsd_xcont

# 方法3：简单的mtDNA异质性检查
# 在高覆盖位点统计minor allele频率
samtools mpileup -f reference_mt.fa sample.mt.bam | \
    awk '$4>=50' | \  # 只看覆盖>=50的位点
    python calc_heteroplasmy.py > mt_heteroplasmy.txt

# 真实aDNA：大多数位点minor allele freq < 1-2%
# 高污染：多个位点minor allele freq > 5%

第六步：性别判定与单倍群分型¶

白话解释： 从aDNA数据可以判定个体性别（通过X/Y染色体覆盖度比值）和线粒体/Y染色体单倍群（揭示母系/父系血统）。这些是古基因组学研究的基本信息。

技术细节：

# === 性别判定 ===

# 方法：比较X染色体和常染色体的覆盖度比值
# 女性(XX): ratio ≈ 1.0
# 男性(XY): ratio ≈ 0.5

# 使用Ry方法（Y染色体reads比例）
samtools idxstats sample.filtered.bam | \
    awk 'BEGIN{x=0;y=0;a=0} 
         $1~/^chr[0-9]/{a+=$3}
         $1=="chrX"{x=$3}
         $1=="chrY"{y=$3}
         END{
           total=a+x+y;
           ry=y/total;
           rx=x/total;
           if(ry>0.05) print "Male (Ry="ry")";
           else if(ry<0.01) print "Female (Ry="ry")";
           else print "Undetermined (Ry="ry")"
         }'

# === 线粒体单倍群分型 ===
# 1. 生成线粒体共识序列
samtools mpileup -uf reference_mt.fa sample.mt.bam | \
    bcftools call -c --ploidy 1 -O v -o sample.mt.vcf

# 2. 使用HaploGrep2分型
# 在线提交VCF: https://haplogrep.i-med.ac.at/
# 或命令行：
java -jar haplogrep.jar classify \
    --input sample.mt.vcf \
    --format vcf \
    --output haplogroup_result.txt

# === Y染色体单倍群（男性）===
# 提取Y染色体变异
samtools mpileup -uf reference.fasta -r chrY sample.filtered.bam | \
    bcftools call -c --ploidy 1 -O v -o sample.chrY.vcf

# 使用yhaplo/Y-LineageTracker分型

第七步：群体遗传学分析¶

白话解释： 古基因组学的终极目标是将古代个体放入人类群体遗传学的大图景中——他们与哪些现代人群最相近？经历了怎样的迁徙和混合？这需要将aDNA数据与现代人群数据整合进行PCA、ADMIXTURE和f-statistics分析。

技术细节：

# === 群体遗传分析 ===

# 1. 随机抽样基因型（pseudo-haploid calling）
# aDNA覆盖度低，通常随机抽取一条reads作为该位点基因型
# 使用ANGSD
angsd -i sample.filtered.bam \
    -sites 1240k_sites.txt \
    -doHaploCall 1 \
    -doCounts 1 \
    -minMapQ 30 -minQ 20 \
    -out sample_haplocall

# 2. 或使用pileupCaller
samtools mpileup -R -B -q 30 -Q 20 \
    -l 1240k_snps.bed \
    sample.filtered.bam | \
    pileupCaller --randomHaploid --sampleNames sample1 \
    --samplePopName Ancient \
    --eigenstratOut sample_eigen

# 3. 与参考数据集合并
mergeit -p merge_par.txt
# 合并ancient + modern reference (如Human Origins / 1240K)

# 4. PCA投影
# 用smartpca，将古代样本投影到现代人群PCA空间
cat > pca_par.txt << 'EOF'
genotypename: merged.geno
snpname: merged.snp
indivname: merged.ind
evecoutname: pca_results.evec
evaloutname: pca_results.eval
numoutevec: 10
lsqproject: YES
numthreads: 8
EOF
smartpca -p pca_par.txt

# 5. ADMIXTURE分析
# 先转PLINK格式
convertf -p convert_par.txt
# 运行ADMIXTURE
for K in 2 3 4 5 6 7 8; do
    admixture --cv merged.bed $K -j16 | tee admixture_K${K}.log
done

# 6. f-statistics (AdmixTools2)
# f3(Ancient; Pop1, Pop2) - 混合检验
# f4(Ancient, Ref; Pop1, Pop2) - 亲缘检验
# D-statistics
qp3Pop -p f3_par.txt   # f3-statistics
qpDstat -p d_par.txt   # D-statistics

实战命令速查¶

# aDNA完整流程
# 1.预处理
AdapterRemoval --file1 R1.fq.gz --file2 R2.fq.gz --collapse --basename sample
# 2.比对
bwa aln -n 0.01 -o 2 -l 1024 ref.fa sample.collapsed.gz > sample.sai
bwa samse ref.fa sample.sai sample.collapsed.gz | samtools sort -o sample.bam
# 3.去重
picard MarkDuplicates INPUT=sample.bam OUTPUT=sample.dedup.bam REMOVE_DUPLICATES=true
# 4.损伤认证
mapDamage -i sample.dedup.bam -r ref.fa -d mapdamage_out
# 5.PMD过滤
samtools view -h sample.dedup.bam | pmdtools --threshold 3 --header | samtools view -Sb - > sample.pmd.bam
# 6.污染评估
schmutzi ... / angsd contamination

面试常问点¶

Q1: 如何认证一个样本确实含有古DNA而非现代污染？¶

A: 多重证据标准：(1) mapDamage显示典型的5' C→T末端损伤模式；(2) 片段长度分布符合aDNA特征（短片段为主）；(3) 污染估计低（<5%）；(4) 线粒体为单一单倍群（无异质性）；(5) 如有多个提取，不同提取/建库批次结果一致；(6) 空白对照无扩增产物。

Q2: UDG处理的文库为什么看不到损伤模式？这样的数据还能用吗？¶

A: UDG（Uracil-DNA-glycosylase）酶会去除脱氨基产生的尿嘧啶，消除C→T损伤信号。UDG处理的文库无法通过损伤模式认证（mapDamage图形平坦），但好处是消除了损伤导致的假SNP，提高变异检测准确性。现代做法是对同一样本做两种文库：非UDG库用于认证，UDG库用于基因分型。

Q3: 为什么aDNA比对推荐BWA aln而非BWA-MEM？¶

A: BWA-MEM对短reads（<70bp）效果不佳——它的种子长度默认19bp，对30-50bp的aDNA片段可能只能设一个种子，降低灵敏度。BWA aln的backtracking策略对短reads更有效。但对>70bp的reads（如UDG处理后保留的较长片段），BWA-MEM也是可接受的。

Q4: 什么是pseudo-haploid calling？为什么aDNA分析常用？¶

A: 由于aDNA覆盖度通常很低（0.5-5x），大多数位点只有1-2条reads覆盖。在这种情况下无法可靠地进行二倍体基因型调用（无法区分杂合与纯合）。pseudo-haploid calling策略是：每个位点随机抽取一条read的碱基作为该位点的"等位基因"。这样引入了随机噪声但不引入系统偏差，适合下游的群体遗传学分析。

Q5: 古DNA研究中如何处理群体遗传分析的低覆盖问题？¶

A: (1) 使用pseudo-haploid calling避免基因型调用偏差；(2) 将古代样本投影（projection）到由高质量现代样本定义的PCA空间，而非联合计算；(3) 使用f-statistics等对缺失数据鲁棒的统计量；(4) ADMIXTURE分析时用投影模式（--supervised）。

易错点¶

1. 对UDG处理的样本强行做损伤认证¶

UDG处理的文库看不到C→T信号，不代表不是古DNA。需要配合非UDG文库做认证，或使用其他认证标准。

2. 使用现代DNA比对参数¶

aDNA比对不能用默认参数。需要：关闭seeding（-l 1024）、允许更多mismatch（-n 0.01）、保留短reads（不设length filter或设极低如30bp）。

3. 忽略PCR重复的严重程度¶

aDNA文库PCR重复率常达50-90%（因为起始DNA极少需大量PCR扩增）。不去重会严重高估覆盖度。

4. 损伤位点被当作真实SNP¶

未经UDG处理的数据中，损伤产生的C→T会被误认为真实变异。解决方案：(1) 使用mapDamage rescale降低损伤位点质量；(2) 过滤transition变异只用transversion；(3) 使用UDG处理的文库做基因分型。

5. 污染评估不充分就进行群体遗传分析¶

污染>5%的样本如果不处理直接做PCA/ADMIXTURE，会导致古代样本向现代人群偏移，产生虚假的"混合"信号。

补充知识¶

aDNA建库策略¶

单链文库（ssDNA library）：效率最高，适合降解严重的样本
双链文库（dsDNA library）：传统方法，较成熟
UDG处理文库：消除损伤，适合基因分型
Half-UDG（部分UDG）：保留末端几个碱基的损伤信号用于认证

常用aDNA资源¶

1240K capture panel：120万个人类群体遗传学信息SNP的靶向捕获
Allen Ancient DNA Resource (AADR)：已发表古基因组数据集
Human Origins dataset：参考现代人群基因型

引用推荐¶

mapDamage2: Jónsson et al., Bioinformatics, 2013
PMDtools: Skoglund et al., PNAS, 2014
EAGER pipeline: Peltzer et al., Genome Biology, 2016
综述: Orlando et al., Nature Reviews Genetics, 2021

古DNA分析¶

一句话概述¶

核心知识点总览¶

各步骤详解¶

第一步：古DNA的特殊性质¶

第二步：数据预处理¶

第三步：mapDamage2损伤认证¶

第四步：PMDtools筛选古代reads¶

第五步：污染评估¶

第六步：性别判定与单倍群分型¶

第七步：群体遗传学分析¶

实战命令速查¶

面试常问点¶

Q1: 如何认证一个样本确实含有古DNA而非现代污染？¶

Q2: UDG处理的文库为什么看不到损伤模式？这样的数据还能用吗？¶

Q3: 为什么aDNA比对推荐BWA aln而非BWA-MEM？¶

Q4: 什么是pseudo-haploid calling？为什么aDNA分析常用？¶

Q5: 古DNA研究中如何处理群体遗传分析的低覆盖问题？¶

易错点¶

1. 对UDG处理的样本强行做损伤认证¶

2. 使用现代DNA比对参数¶

3. 忽略PCR重复的严重程度¶

4. 损伤位点被当作真实SNP¶

5. 污染评估不充分就进行群体遗传分析¶

补充知识¶

aDNA建库策略¶

常用aDNA资源¶

引用推荐¶

📚 相关文章推荐