基因组选择性清除¶

一句话概述¶

通过计算iHS、XP-EHH、Fst和Tajima's D等统计量检测基因组中受到自然选择或人工选择作用的区域（选择性清除），是群体遗传学揭示适应性进化遗传基础的核心分析方法。

核心知识点表格¶

知识点	关键内容	常用工具
选择信号特征	单倍型频率升高、多态性降低、连锁不平衡增强	-
iHS	群体内单倍型纯合性检测正选择	selscan, rehh
XP-EHH	群体间单倍型差异检测选择	selscan, rehh
Fst	群体分化指数	VCFtools, PLINK, pixy
Tajima's D	等位基因频率谱偏移检测	VCFtools, ANGSD
CLR/SweeD	经典清除(completed sweep)	SweeD, SweepFinder2
复合指标	多指标联合提高检测力	hapbin, XP-nSL
基因功能注释	候选区域功能解读	GREAT, g:Profiler

各步骤详解¶

第一步：理解选择性清除原理¶

白话解释： 当一个有利突变出现并受到正向选择时，这个突变的频率会在群体中迅速增加。由于基因重组来不及打破连锁，这个突变附近的中性变异也会"搭便车"一起升高频率。结果就是：这个区域的遗传多样性降低、出现高频长单倍型、群体间分化增加。检测这些特征就能找到受选择的区域。

技术细节： - 经典清除(hard sweep): 新突变从头出现并被选择固定 - 软清除(soft sweep): 已有变异(standing variation)被选择 - 不完全清除(incomplete sweep): 正在进行中的选择 - 不同检测方法对应不同选择模式和时间尺度 - 群体结构(population structure)是最大混淆因素

iHS/XP-EHH的时间敏感性：

Tajima's D    → 检测较古老的选择(>10万年)
Fst           → 检测群体分化后的选择
iHS           → 检测近期不完全清除(~3万年内)
XP-EHH        → 检测近期一个群体固定的选择
nSL           → 对重组率不敏感，更稳健

第二步：数据准备与质控¶

白话解释： 选择信号检测需要高质量的基因型数据。通常使用全基因组SNP数据（来自SNP芯片或WGS），需要进行严格的质控、分群、phasing（单倍型推断）。

代码示例：

# 起始数据：VCF文件（多群体）
# 质控流程

# 1. 基本过滤
bcftools view -m2 -M2 --min-ac 1 \   # 只保留双等位基因SNP
    -i 'QUAL>30 && INFO/DP>10' \
    input.vcf.gz | \
bcftools filter -e 'F_MISSING > 0.1' \   # 缺失率<10%
    -Oz -o filtered.vcf.gz

# 2. 按群体分组
bcftools view -S population_A.samples filtered.vcf.gz -Oz -o popA.vcf.gz
bcftools view -S population_B.samples filtered.vcf.gz -Oz -o popB.vcf.gz

# 3. LD-based pruning（用于Fst等，但iHS/XP-EHH不需要）
plink --vcf popA.vcf.gz --indep-pairwise 50 5 0.5 --out ld_prune
# 注意：iHS/XP-EHH需要完整的连锁信息，不能做LD pruning!

# 4. Phasing（单倍型推断）- iHS/XP-EHH必须
# 使用SHAPEIT5（SHAPEIT4的升级替代版）或Beagle
# SHAPEIT5 phase_common（原SHAPEIT4已被替代）:
phase_common --input popA.vcf.gz \
    --map genetic_map_chr{}.txt \
    --region chr1 \
    --output popA_phased.vcf.gz \
    --thread 8

# Beagle5:
java -jar beagle.jar \
    gt=popA.vcf.gz \
    map=plink.chr1.GRCh38.map \
    out=popA_phased \
    nthreads=8

# 5. 准备遗传图谱（genetic map）
# 从SHAPEIT或HapMap下载重组率图谱
# 格式: chr position rate(cM/Mb) cM

第三步：iHS计算（群体内选择检测）¶

白话解释： iHS (integrated Haplotype Score) 比较一个SNP位点上祖先等位基因和衍生等位基因各自单倍型延伸的长度。如果衍生等位基因的单倍型异常长（说明它是最近才升高频率的，连锁还没被重组打破），iHS就会是一个极端值，提示该位点可能受到了正选择。

技术细节： - 基于EHH (Extended Haplotype Homozygosity) 概念 - iHH = EHH曲线下面积（积分） - iHS = ln(iHH_ancestral / iHH_derived)，标准化后使用 - 适合检测中等频率(MAF 20-80%)的不完全清除 - 需要祖先等位基因信息（ancestral allele）

代码示例：

# 使用selscan计算iHS
conda install -c bioconda selscan

# 准备输入文件
# 需要：phased VCF + genetic map + ancestral allele annotation

# 标注祖先等位基因（使用Ensembl ancestral allele）
# 在VCF中添加AA (Ancestral Allele) INFO字段

# 运行selscan计算iHS
selscan --ihs \
    --vcf popA_phased.vcf.gz \
    --map genetic_map_chr1.txt \
    --out popA_chr1 \
    --threads 8

# 输出: popA_chr1.ihs.out
# 格式(6列): locus_id, physical_pos, freq_1, ihh_1, ihh_0, unstd_ihs

# 标准化（按频率分bin标准化）
norm --ihs --files popA_chr1.ihs.out \
    --bins 100 \
    --qbins 10

# 输出: popA_chr1.ihs.out.100bins.norm
# 最后一列是标准化的iHS值

# 使用R的rehh包
Rscript << 'EOF'
library(rehh)

# 读取phased数据
hap <- data2haplohh(hap_file = "popA_phased.vcf.gz",
    map_file = "snp_map.txt",
    polarize_vcf = TRUE,  # 使用AA字段极化
    chr.name = "1")

# 计算EHH和iHS
scan <- scan_hh(hap, threads = 8)
ihs <- ihh2ihs(scan, freqbin = 0.05)

# 可视化Manhattan plot
manhattanplot(ihs, pval = TRUE, threshold = 4,
    main = "iHS genome scan")

# 筛选候选区域
candidates <- ihs$ihs[abs(ihs$ihs$IHS) > 4, ]
EOF

第四步：XP-EHH计算（群体间选择检测）¶

白话解释： XP-EHH比较同一位点在两个群体中单倍型的延伸长度。如果某个等位基因在群体A中已经被选择到高频或固定，但在群体B中没有，那么该位点在群体A中的单倍型会更长。XP-EHH适合检测已经（接近）完成的选择清除。

代码示例：

# 使用selscan计算XP-EHH
selscan --xpehh \
    --vcf popA_phased.vcf.gz \
    --vcf-ref popB_phased.vcf.gz \
    --map genetic_map_chr1.txt \
    --out popA_vs_popB_chr1 \
    --threads 8

# 标准化
norm --xpehh --files popA_vs_popB_chr1.xpehh.out \
    --bins 100 --qbins 10

# 正值: popA中受选择（单倍型更长）
# 负值: popB中受选择

# 使用rehh包
Rscript << 'EOF'
library(rehh)

# 读取两个群体数据
hapA <- data2haplohh(hap_file="popA_phased.vcf.gz", chr.name="1")
hapB <- data2haplohh(hap_file="popB_phased.vcf.gz", chr.name="1")

# 计算EHH
scanA <- scan_hh(hapA, threads=8)
scanB <- scan_hh(hapB, threads=8)

# 计算XP-EHH
xpehh <- ies2xpehh(scanA, scanB, popname1="PopA", popname2="PopB")

# 可视化
manhattanplot(xpehh, pval=TRUE, threshold=4)

# 候选区域
sig_xpehh <- xpehh$xpehh[abs(xpehh$xpehh$XPEHH_PopA_PopB) > 4, ]
EOF

第五步：Fst和Tajima's D计算¶

白话解释： Fst衡量两个群体间的遗传分化程度。高Fst区域意味着这些位点在两个群体中的等位基因频率差异大，可能是方向性选择的结果。Tajima's D检测等位基因频率谱的偏移——负值表示低频变异过多（提示选择性清除后恢复），正值表示中频变异过多（提示平衡选择）。

代码示例：

# === Fst计算 ===

# 方法1: VCFtools (Weir & Cockerham Fst)
vcftools --gzvcf filtered.vcf.gz \
    --weir-fst-pop population_A.txt \
    --weir-fst-pop population_B.txt \
    --fst-window-size 50000 \
    --fst-window-step 25000 \
    --out fst_results

# 输出: fst_results.windowed.weir.fst
# 格式: CHROM BIN_START BIN_END N_VARIANTS WEIGHTED_FST MEAN_FST

# 方法2: pixy (处理缺失数据更好)
conda install -c bioconda pixy

# pixy需要全位点VCF（含invariant sites）
pixy --stats fst \
    --vcf all_sites.vcf.gz \
    --populations populations.txt \
    --window_size 50000 \
    --n_cores 8 \
    --output_prefix pixy_fst

# === Tajima's D计算 ===

# VCFtools
vcftools --gzvcf popA.vcf.gz \
    --TajimaD 50000 \
    --out tajima_popA

# ANGSD (适合低覆盖度数据)
angsd -bam bamlist_popA.txt \
    -doThetas 1 -doSaf 1 \
    -pest popA.sfs \
    -anc ancestral.fa \
    -GL 1 -nThreads 16 \
    -out angsd_popA

# 计算Tajima's D
thetaStat do_stat angsd_popA.thetas.idx \
    -win 50000 -step 25000 \
    -outnames tajima_angsd

# === 合并多指标筛选 ===
Rscript << 'EOF'
library(tidyverse)

# 读取各指标
fst <- read.table("fst_results.windowed.weir.fst", header=TRUE)
tajd <- read.table("tajima_popA.Tajima.D", header=TRUE)
ihs_data <- read.table("ihs_normalized.txt", header=TRUE)

# 标准化各指标并取交集
# 候选区域标准: top 1% Fst AND (iHS > 4 OR Tajima's D < -2)
fst_threshold <- quantile(fst$WEIGHTED_FST, 0.99, na.rm=TRUE)
candidates <- fst %>%
    filter(WEIGHTED_FST > fst_threshold) %>%
    inner_join(tajd, by=c("CHROM"="CHROM")) %>%
    filter(TajimaD < -2)

write.csv(candidates, "selection_candidates.csv", row.names=FALSE)
EOF

第六步：候选区域功能注释¶

白话解释： 找到受选择的基因组区域后，需要知道这些区域包含什么基因，这些基因的功能是什么，是否与已知的适应性表型相关。

代码示例：

# 将候选区域转为BED格式
awk 'NR>1 {print $1"\t"$2"\t"$3}' selection_candidates.csv > candidates.bed

# 查找区域内的基因
bedtools intersect -a candidates.bed -b genes.bed -wa -wb > candidate_genes.txt

# GO/KEGG富集分析
Rscript << 'EOF'
library(clusterProfiler)
library(org.Hs.eg.db)

genes <- read.table("candidate_genes.txt")$V7  # 基因名列
gene_ids <- bitr(genes, fromType="SYMBOL", toType="ENTREZID", OrgDb=org.Hs.eg.db)

# GO富集
ego <- enrichGO(gene = gene_ids$ENTREZID,
    OrgDb = org.Hs.eg.db, ont = "BP",
    pAdjustMethod = "BH", pvalueCutoff = 0.05)

# KEGG富集
ekegg <- enrichKEGG(gene = gene_ids$ENTREZID,
    organism = "hsa", pvalueCutoff = 0.05)

# 可视化
dotplot(ego, showCategory=20)
EOF

实战命令¶

#!/bin/bash
# === 全基因组选择信号扫描流程 ===
set -e
THREADS=16

# 对每条染色体并行处理
for chr in $(seq 1 22); do
    # Phasing（SHAPEIT5 phase_common，原SHAPEIT4的替代）
    phase_common --input popA_chr${chr}.vcf.gz \
        --map genetic_map_chr${chr}.txt \
        --output popA_chr${chr}_phased.vcf.gz --thread $THREADS

    # iHS
    selscan --ihs --vcf popA_chr${chr}_phased.vcf.gz \
        --map genetic_map_chr${chr}.txt --out ihs_chr${chr} --threads $THREADS

    # XP-EHH
    selscan --xpehh --vcf popA_chr${chr}_phased.vcf.gz \
        --vcf-ref popB_chr${chr}_phased.vcf.gz \
        --map genetic_map_chr${chr}.txt --out xpehh_chr${chr} --threads $THREADS
done

# 合并和标准化
norm --ihs --files ihs_chr*.ihs.out --bins 100
norm --xpehh --files xpehh_chr*.xpehh.out --bins 100

# Fst (全基因组滑动窗口)
vcftools --gzvcf merged.vcf.gz --weir-fst-pop popA.txt --weir-fst-pop popB.txt \
    --fst-window-size 50000 --fst-window-step 25000 --out genome_fst

面试常问点¶

Q1: iHS和XP-EHH分别适合检测什么类型的选择？ A: iHS检测群体内正在进行的不完全清除（衍生等位基因频率20-80%），敏感于最近~3万年内的选择。XP-EHH检测一个群体中已完成或接近完成的清除（频率>80%），适合检测群体分化后的差异性选择。

Q2: 为什么Fst不能单独用来检测选择？ A: Fst高值可能由随机遗传漂变（尤其是小群体）、群体瓶颈、迁移模式等非选择因素导致。需要结合其他指标（如单倍型延伸长度）来区分选择和漂变。

Q3: soft sweep和hard sweep的检测差异？ A: Hard sweep产生经典的选择信号（多态性剧烈下降、长单倍型），iHS/XP-EHH能有效检测。Soft sweep信号较弱（多个单倍型同时升频），传统方法敏感性低，需要用nSL、H12/H2H1等新方法。

Q4: 群体结构如何影响选择信号检测？ A: 未校正的群体结构会产生大量假阳性Fst异常值。解决：1）PCA确认无亚结构或校正；2）使用基于似然比的方法（如BayeScan）内置结构校正；3）排除近亲个体。

Q5: 如何区分正选择和背景选择(background selection)？ A: 背景选择（清除有害突变连带降低附近中性多态性）也降低diversity，但不产生高频长单倍型。iHS/XP-EHH对背景选择不敏感。结合多指标：真正的正选择应同时显示高Fst + 极端iHS + 低diversity。

易错点¶

未做Phasing就算iHS/XP-EHH：这些方法依赖单倍型信息，必须先进行statistical phasing。
遗传图谱不匹配：iHS需要遗传距离（cM），用物理距离会产生错误结果。确保genetic map与参考基因组版本一致。
祖先等位基因标注错误：iHS需要正确极化（区分祖先/衍生等位基因），错误极化会反转信号方向。
窗口大小选择不当：Fst/Tajima's D的窗口太小则噪声大，太大则信号被稀释。建议50kb窗口、25kb步进。
忽略局部重组率变异：低重组率区域自然产生长单倍型，可能产生假阳性。需要用遗传距离而非物理距离。
多重检验未校正：全基因组数百万SNP测试，需要严格的多重检验校正或使用经验分位数阈值。

补充知识¶

选择信号检测方法时间尺度¶

方法	检测时间尺度	选择类型	数据需求
Tajima's D / π	古老选择(>10^5年)	完成的清除	单群体
CLR/SweeD	完成的经典清除	Hard sweep	单群体SFS
Fst	群体分化后	方向性选择	两群体+
iHS/nSL	近期不完全(~3万年)	进行中选择	单群体+phased
XP-EHH	近期一方固定	完成的选择	两群体+phased
H12/H2H1	近期	Hard+soft sweep	单群体+phased

常用选择信号经验阈值¶

|iHS| > 2.0    → 候选(top ~5%)
|iHS| > 3.0    → 强信号(top ~0.5%)
|iHS| > 4.0    → 极显著
|XP-EHH| > 2.0 → 候选
Fst > top 1%    → 候选
Tajima's D < -2 → 显著偏离中性
PBS > top 0.1%  → 强候选(三群体设计)

基因组选择性清除¶

一句话概述¶

核心知识点表格¶

各步骤详解¶

第一步：理解选择性清除原理¶

第二步：数据准备与质控¶

第三步：iHS计算（群体内选择检测）¶

第四步：XP-EHH计算（群体间选择检测）¶

第五步：Fst和Tajima's D计算¶

第六步：候选区域功能注释¶

实战命令¶

面试常问点¶

易错点¶

补充知识¶

选择信号检测方法时间尺度¶

常用选择信号经验阈值¶

📚 相关文章推荐