古DNA（aDNA）数据分析¶

一句话概述：古DNA分析从古代标本中提取高度降解的 DNA 片段，通过专门的生信流程（损伤校正、低深度处理、污染评估）重建古代群体的遗传历史。

核心知识点速查表¶

概念	白话解释
古DNA（aDNA）	从考古遗骸（骨骼、牙齿）中提取的DNA
DNA损伤	古DNA特有的化学修饰（C→T/G→A脱氨基）
mapDamage	评估和可视化 DNA 损伤模式的工具
PMD	死后损伤（Post-Mortem Damage），区分内源/污染
捕获测序	用探针富集目标区域（节省测序成本）
EAGER	古DNA分析的一站式流程
contamination	现代DNA混入古DNA样本的污染
基因型似然值	低深度数据不"确定"基因型，用概率表示

一、古DNA的特殊性（白话版）¶

比喻：古DNA就像一本被水泡过、被虫蛀过、碎成几千片的古书。你需要： 1. 把碎片拼回去（比对到参考基因组） 2. 判断哪些碎片是原书的（内源DNA），哪些是后来混进去的（污染） 3. 修复水渍造成的字迹变形（DNA损伤校正）

古DNA三大挑战： - 高度降解：片段通常只有 30~80bp（现代DNA是几百kb） - DNA损伤：C→U（测序读为T），在片段末端最严重 - 高污染：样本中 >90% 可能是细菌或现代人的DNA

二、分析流程¶

整体流程¶

原始测序 → 接头去除 → 比对(短片段模式) → 去重 → 损伤评估 → 污染估计 → 基因型推断 → 下游分析

第 1 步：安装核心工具¶

# 安装古DNA分析核心工具
conda create -n ancient_dna python=3.10                   # 创建环境
conda activate ancient_dna                                 # 激活环境

conda install -c bioconda \
    fastp \                   # 质控和接头去除
    bwa \                     # 比对（BWA-aln模式，适合短片段）
    samtools \                # BAM文件操作
    picard \                  # 去除PCR重复
    mapdamage2 \              # DNA损伤评估
    angsd \                   # 低深度基因型分析
    schmutzi \                # 污染估计
    eager \                   # 一站式流程（可选）
    htslib                    # VCF/BCF操作

第 2 步：质控与接头去除¶

# 古DNA文库通常是双端测序，但片段很短可能有 read-through
# 需要做接头修剪和 overlap 合并

# 用 fastp 做质控和接头去除
fastp -i raw_R1.fastq.gz \          # 正向读
      -I raw_R2.fastq.gz \          # 反向读
      -o clean_R1.fastq.gz \        # 清理后的正向读
      -O clean_R2.fastq.gz \        # 清理后的反向读
      --merge \                      # 合并重叠的双端读
      --merged_out merged.fastq.gz \ # 合并后的读
      --overlap_len_require 11 \     # 最小重叠长度11bp
      --adapter_sequence AGATCGGAAGAGCACACGTCTGAACTCCAGTCA \   # Illumina接头
      --adapter_sequence_r2 AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT \ # 反向接头
      -h qc_report.html \           # HTML质控报告
      -j qc_report.json             # JSON质控报告

# 或用 AdapterRemoval（古DNA领域更常用）
AdapterRemoval --file1 raw_R1.fastq.gz \
               --file2 raw_R2.fastq.gz \
               --collapse \                  # 合并重叠读
               --trimns \                    # 修剪N碱基
               --trimqualities \             # 修剪低质量碱基
               --minquality 20 \             # 最低质量20
               --minlength 30 \              # 最短长度30bp
               --basename sample_clean       # 输出前缀

第 3 步：比对到参考基因组¶

# 古DNA片段很短（30-80bp），需要用 BWA-aln 而不是 BWA-MEM
# BWA-aln 对短序列更准确

# 建索引
bwa index reference.fa              # 建立参考基因组索引

# BWA-aln 比对（适合 <70bp 的片段）
bwa aln -l 1024 \                   # 禁用种子（短片段不适合种子策略）
        -n 0.01 \                    # 允许较高的错误率（古DNA质量差）
        -o 2 \                       # 允许2个gap
        -t 8 \                       # 8个线程
        reference.fa \
        merged.fastq.gz \
        > sample.sai                 # 输出中间文件

bwa samse reference.fa \             # 单端模式（合并后的读）
          sample.sai \
          merged.fastq.gz \
          | samtools sort -@ 8 \     # 排序
          -o sample_sorted.bam       # 输出BAM

samtools index sample_sorted.bam     # 建BAM索引

# 如果用BWA-MEM2（>70bp的片段）
bwa-mem2 mem -t 8 reference.fa merged.fastq.gz | samtools sort -o sample.bam

第 4 步：去除 PCR 重复¶

# 古DNA文库中PCR重复率通常很高（>50%）
# 必须去重

# 用 Picard MarkDuplicates
picard MarkDuplicates \
    INPUT=sample_sorted.bam \        # 输入BAM
    OUTPUT=sample_dedup.bam \        # 输出去重BAM
    METRICS_FILE=dedup_metrics.txt \ # 重复率统计
    REMOVE_DUPLICATES=true \         # 删除重复读
    VALIDATION_STRINGENCY=LENIENT    # 宽松验证

samtools index sample_dedup.bam      # 重新建索引

# 查看重复率
grep "PERCENT_DUPLICATION" dedup_metrics.txt
# 古DNA重复率 30-80% 很常见

第 5 步：DNA 损伤评估（mapDamage2）¶

# mapDamage2 评估古DNA特有的损伤模式
mapDamage -i sample_dedup.bam \      # 输入BAM
          -r reference.fa \           # 参考基因组
          -d mapDamage_results \      # 输出目录
          --merge-reference-sequences # 合并染色体

# 输出文件：
# mapDamage_results/Fragmisincorporation_plot.pdf — 错配率图
# 典型古DNA模式：5'端 C→T 升高，3'端 G→A 升高
# 5'端第1个碱基 C→T 率 > 10% → 说明确实是古DNA

# mapDamage 还可以做损伤校正（rescaling）
mapDamage -i sample_dedup.bam \
          -r reference.fa \
          -d mapDamage_results \
          --rescale                    # 校正损伤碱基的质量值
# 输出校正后的BAM文件

第 6 步：污染估计¶

# 方法1：用线粒体DNA估计污染率
# schmutzi 工具
schmutzi --library single \          # 单链文库
         -t 4 \                       # 4个线程
         output_prefix \              # 输出前缀
         reference_mt.fa \            # 线粒体参考序列
         sample_dedup.bam             # 输入BAM

# 输出：污染率估计值
# 污染率 < 5% 通常可以接受

# 方法2：用X染色体（男性样本）
# 男性只有1条X染色体，如果X上出现杂合位点 → 污染
# ANGSD 可以做这个分析
angsd -i sample_dedup.bam \
      -r X: \                         # 只分析X染色体
      -doCounts 1 \
      -iCounts 1 \
      -minMapQ 30 \
      -minQ 20

第 7 步：低深度基因型分析（ANGSD）¶

# 古DNA深度通常很低（0.1~5X），不能可靠地"确定"基因型
# 用 ANGSD 基于基因型似然值分析

# 计算基因型似然值（genotype likelihoods）
angsd -bam bam_list.txt \           # BAM文件列表（每行一个样本）
      -GL 1 \                        # GATK模型计算GL
      -doGlf 2 \                    # 输出beagle格式GL
      -doMajorMinor 1 \             # 从数据推断major/minor等位基因
      -doMaf 1 \                     # 计算等位基因频率
      -SNP_pval 1e-6 \              # SNP显著性阈值
      -minMapQ 30 \                  # 最低比对质量
      -minQ 20 \                     # 最低碱基质量
      -nThreads 8 \                  # 线程数
      -out angsd_output              # 输出前缀

# 用 PCAngsd 做PCA（适合低深度数据）
pip install pcangsd                   # 安装PCAngsd
pcangsd -b angsd_output.beagle.gz \  # 输入beagle格式GL
        -o pca_output \               # 输出前缀
        -t 8                           # 线程数

# 用 ngsAdmix 做群体结构（类似ADMIXTURE但用GL）
NGSadmix -likes angsd_output.beagle.gz \  # 输入GL
         -K 3 \                            # K值
         -outfiles admix_K3 \              # 输出前缀
         -P 8                               # 线程数

三、一站式流程 EAGER¶

# EAGER（Efficient Ancient GEnome Reconstruction）
# 是古DNA分析的一站式Nextflow流程

# 安装
pip install eager                     # 安装EAGER

# 或用 nf-core/eager（推荐）
nextflow run nf-core/eager \
    -profile conda \
    --input samplesheet.tsv \          # 样本信息表
    --fasta reference.fa \             # 参考基因组
    --run_bam_filtering \              # BAM过滤
    --run_mapdamage_rescaling \        # 损伤校正
    --run_genotyping                   # 基因型推断

四、常见报错与解决¶

报错信息	原因	解决方法
比对率 < 1%	大部分reads是细菌DNA	正常（古DNA内源率通常很低）
mapDamage 无损伤模式	可能不是古DNA或文库做了UDG处理	检查文库制备方法
`too few reads`	深度太低	做捕获测序富集目标区域
污染率 > 10%	现代DNA污染严重	排除样本或用PMDtools过滤
ANGSD 报内存错误	样本或区域太大	分区域处理或减少样本

五、面试高频问题¶

Q1：古DNA分析和现代DNA分析的主要区别？

三个主要区别：(1) 需要评估和校正DNA损伤（C→T）；(2) 深度极低，必须用基因型似然值而不是确定基因型；(3) 需要评估和过滤污染（现代DNA混入）。

Q2：如何判断一个样本确实是古DNA？

看mapDamage的损伤模式：真正的古DNA在5'端有明显的C→T替换升高（通常>10%），3'端有G→A升高。这是古DNA的"指纹"，现代DNA污染不会有这种模式。

Q3：为什么古DNA分析不能用GATK等标准变异检测流程？

因为：(1) GATK等工具假设较高的覆盖度（>20X），古DNA通常只有0.1~5X；(2) DNA损伤会导致大量假阳性SNP（C→T被误判为突变）；(3) 需要用基因型似然值的方法（如ANGSD）替代硬基因型调用。

六、速查表¶

# === 古DNA分析速查 ===

# 1. 质控和接头去除
AdapterRemoval --file1 R1.fq.gz --file2 R2.fq.gz --collapse --basename clean

# 2. 比对（短片段用BWA-aln）
bwa aln -l 1024 -n 0.01 ref.fa reads.fq.gz > aln.sai
bwa samse ref.fa aln.sai reads.fq.gz | samtools sort -o sorted.bam

# 3. 去重
picard MarkDuplicates INPUT=sorted.bam OUTPUT=dedup.bam REMOVE_DUPLICATES=true

# 4. 损伤评估
mapDamage -i dedup.bam -r ref.fa -d results/ --rescale

# 5. 污染估计
schmutzi --library single output ref_mt.fa dedup.bam

# 6. 低深度分析（ANGSD）
angsd -bam list.txt -GL 1 -doGlf 2 -doMajorMinor 1 -doMaf 1 -out angsd
pcangsd -b angsd.beagle.gz -o pca
NGSadmix -likes angsd.beagle.gz -K 3 -outfiles K3

# 古DNA质量指标参考：
# 内源率：>1% 可用，>10% 较好
# 损伤率：5'C→T >10% 确认古DNA
# 污染率：<5% 可接受
# 平均深度：>0.5X 可做PCA/ADMIXTURE
# 平均深度：>5X 可做基因型确定

参考资料：EAGER pipeline (Yates et al. 2021)、mapDamage2 (Jonsson et al. 2013)、ANGSD (Korneliussen et al. 2014)、nf-core/eager