连锁不平衡（LD）分析与可视化¶
一句话概述：连锁不平衡（LD）衡量两个遗传位点等位基因间的非随机关联程度，通过 PLINK 和 PopLDdecay 计算 r² 和 D'，绘制 LD 衰减曲线揭示群体历史。
核心知识点速查表¶
| 概念 | 白话解释 |
|---|---|
| LD（连锁不平衡） | 两个位点的等位基因"抱团"出现，不是随机组合 |
| r² | LD 的常用度量（0~1），越大表示两个位点越"绑定" |
| D' | 另一种 LD 度量（0~1），更敏感但受样本量影响大 |
| LD 衰减 | 随距离增大 r² 逐渐降到背景水平的过程 |
| LD pruning | 去掉高 LD 的 SNP（用于 PCA、ADMIXTURE 等） |
| LD block | 一段 LD 很高的基因组区域（单倍型块） |
| PopLDdecay | 专门计算 LD 衰减的工具，直接读 VCF |
| 半衰距离 | r² 衰减到一半时的物理距离，反映群体大小 |
一、LD 的原理（白话版）¶
比喻：想象两副扑克牌（两个基因位点），如果你总是同时抽到红心A和黑桃K（特定等位基因组合出现得比随机多），这就是"连锁不平衡"。
LD 产生的原因： 1. 物理连锁：两个位点在染色体上离得近，很少被重组打断 2. 群体事件：瓶颈效应、建立者效应会在无关位点间产生 LD 3. 自然选择：正向选择会在选择位点周围产生高 LD 4. 群体结构：不同群体混合会产生 LD
LD 衰减的含义：自交物种（如水稻）LD 衰减慢（LD block 大），异交物种（如玉米）LD 衰减快（LD block 小）。
二、方法 1：PopLDdecay（推荐）¶
# PopLDdecay：直接从VCF计算LD衰减，速度快、存储省
# 安装
git clone https://github.com/BGI-shenzhen/PopLDdecay.git # 克隆仓库
cd PopLDdecay
make # 编译
# 或 conda 安装
conda install -c bioconda popldecay # conda安装
# 基本用法：计算全基因组 LD 衰减
PopLDdecay -InVCF input.vcf.gz \ # 输入VCF文件（支持gz压缩）
-OutStat LDdecay_all \ # 输出统计结果
-MaxDist 500 \ # 最大距离500kb（默认300kb）
-MAF 0.05 \ # 最低等位基因频率过滤
-Miss 0.1 # 最大缺失率
# 输出文件：LDdecay_all.stat.gz — 压缩的LD统计结果
# 按亚群分别计算（需要样本列表）
# 准备样本列表文件（每行一个样本名）
echo -e "sample1\nsample2\nsample3" > group_A.list
PopLDdecay -InVCF input.vcf.gz \
-OutStat LDdecay_groupA \
-SubPop group_A.list \ # 指定亚群样本列表
-MaxDist 500
# 可视化：单群体 LD 衰减曲线
perl PopLDdecay/bin/Plot_OnePop.pl \
-inFile LDdecay_all.stat.gz \ # 输入统计文件
-output LD_figure # 输出图片前缀
# 输出 LD_figure.png 和 LD_figure.pdf
# 可视化：多群体比较
# 准备群体列表文件（每行：群体名 统计文件路径）
cat << 'EOF' > pop_list.txt
GroupA LDdecay_groupA.stat.gz
GroupB LDdecay_groupB.stat.gz
GroupC LDdecay_groupC.stat.gz
EOF
perl PopLDdecay/bin/Plot_MutiPop.pl \
-inList pop_list.txt \ # 输入群体列表
-output LD_comparison # 输出图片
三、方法 2：PLINK 计算 LD¶
# PLINK 适合计算特定区域的 LD 或 LD block
# 3.1 计算成对 SNP 的 r²
plink --vcf input.vcf.gz \ # 输入VCF
--r2 \ # 计算r²
--ld-window 100 \ # 最多比较100个SNP的窗口
--ld-window-kb 1000 \ # 最大物理距离1Mb
--ld-window-r2 0 \ # 输出所有r²（包括0）
--out ld_pairwise # 输出前缀
# 输出文件：ld_pairwise.ld
# 格式：CHR_A BP_A SNP_A CHR_B BP_B SNP_B R2
# 3.2 LD pruning（为其他分析准备独立SNP集）
# 标记高LD的SNP对
plink --vcf input.vcf.gz \
--indep-pairwise 50 10 0.2 \ # 窗口50 SNP，步长10，r²阈值0.2
--out pruned # 输出前缀
# 输出：pruned.prune.in（保留的SNP）、pruned.prune.out（移除的SNP）
# 提取独立SNP
plink --vcf input.vcf.gz \
--extract pruned.prune.in \ # 只保留通过筛选的SNP
--make-bed \
--out data_ldpruned
# 3.3 计算 LD block（Haploview 算法）
plink --bfile mydata \
--blocks no-pheno-req \ # 不要求表型数据
--blocks-max-kb 200 \ # 最大block大小200kb
--out ld_blocks # 输出前缀
# 输出：ld_blocks.blocks.det — LD block 详细信息
四、用 R 自定义 LD 衰减图¶
# ========== R脚本：LD衰减曲线 ==========
# 读取PLINK的LD输出
ld <- read.table("ld_pairwise.ld", header=TRUE) # 读取LD结果
# 计算SNP对之间的距离
ld$dist_kb <- abs(ld$BP_B - ld$BP_A) / 1000 # 距离转为kb
# 按距离分组计算平均r²
breaks <- seq(0, 500, by=10) # 0~500kb，每10kb一组
ld$bin <- cut(ld$dist_kb, breaks=breaks) # 分组
avg_r2 <- tapply(ld$R2, ld$bin, mean, na.rm=TRUE) # 每组平均r²
mid_dist <- (breaks[-length(breaks)] + breaks[-1]) / 2 # 组中值
# 画 LD 衰减曲线
pdf("LD_decay_custom.pdf", width=10, height=6) # 创建PDF
plot(mid_dist, avg_r2, # 画散点
type="l", lwd=2, col="steelblue", # 蓝色线
xlab="Distance (kb)", # X轴
ylab=expression(r^2), # Y轴（r²上标）
main="LD Decay Curve", # 标题
ylim=c(0, max(avg_r2, na.rm=TRUE))) # Y轴范围
# 标注半衰距离
half_r2 <- max(avg_r2, na.rm=TRUE) / 2 # r²最大值的一半
abline(h=half_r2, col="red", lty=2) # 水平虚线
# 找到半衰距离
half_dist <- mid_dist[which.min(abs(avg_r2 - half_r2))]
abline(v=half_dist, col="red", lty=2) # 垂直虚线
text(half_dist + 20, half_r2 + 0.02, # 标注文字
paste0("Half-decay: ", half_dist, " kb"),
col="red")
dev.off()
五、LD 热图（Heatmap）¶
# 用 LDheatmap R包画特定区域的LD热图
# install.packages("LDheatmap")
library(LDheatmap) # 加载包
# 从PLINK输出构建LD矩阵
# 先提取特定区域的LD
# plink --bfile mydata --chr 1 --from-bp 1000000 --to-bp 2000000 --r2 square --out region_ld
# 读取方形LD矩阵
ld_matrix <- as.matrix(read.table("region_ld.ld")) # 读取r²矩阵
snp_pos <- read.table("region_ld.bim")$V4 # SNP位置
# 画LD热图
pdf("LD_heatmap.pdf", width=10, height=10)
LDheatmap(ld_matrix, # LD矩阵
genetic.distances=snp_pos/1e6, # 位置(Mb)
color=heat.colors(20), # 颜色方案
title="LD Heatmap (r²)", # 标题
add.map=TRUE) # 添加基因组位置
dev.off()
六、常见报错与解决¶
| 报错信息 | 原因 | 解决方法 |
|---|---|---|
PopLDdecay Segmentation fault | 内存不足 | 分染色体计算 |
| PLINK LD 输出文件太大 | 成对比较数量爆炸 | 限制 --ld-window-kb 和 --ld-window |
| r² 在短距离就接近 0 | MAF 过滤太严或样本太少 | 放宽 MAF 阈值，增加样本 |
| LD 衰减曲线不平滑 | SNP 密度太低或分组太粗 | 增加 SNP 密度或调整 bin 大小 |
No variants remaining | 过滤条件太严 | 检查 MAF、缺失率等过滤参数 |
七、面试高频问题¶
Q1：r² 和 D' 的区别？
r² 范围 0~1，受等位基因频率影响，适合 GWAS 设计和 LD 衰减分析。D' 范围 0~1，在样本量小时容易被高估（向 1 偏倚），适合 LD block 定义。一般 LD 衰减用 r²，LD block 用 D'。
Q2：LD 衰减快慢说明什么？
LD 衰减快 → 有效群体大、重组率高（如玉米，衰减距离 ~1kb）。LD 衰减慢 → 有效群体小或自交（如水稻，衰减距离 ~100-200kb）。GWAS 中 LD 衰减快需要更密的标记。
Q3：为什么做 PCA/ADMIXTURE 前要 LD pruning？
因为这些方法假设 SNP 独立。如果不做 LD pruning，基因组中某个高 LD 区域（如着丝粒附近）的 SNP 权重会被放大，导致群体结构被某些区域的信号主导而非全基因组信号。
八、速查表¶
# === LD 分析速查 ===
# PopLDdecay：全基因组LD衰减（推荐）
PopLDdecay -InVCF input.vcf.gz -OutStat LDdecay -MaxDist 500 -MAF 0.05
perl Plot_OnePop.pl -inFile LDdecay.stat.gz -output fig
# PLINK：LD计算与pruning
plink --vcf in.vcf.gz --r2 --ld-window-kb 1000 --out ld_result # 计算r²
plink --vcf in.vcf.gz --indep-pairwise 50 10 0.2 --out pruned # LD pruning
plink --vcf in.vcf.gz --extract pruned.prune.in --make-bed --out data_pruned
# PLINK：LD block
plink --bfile mydata --blocks no-pheno-req --out blocks
# r² 经验阈值：
# r² > 0.8 → 强LD（两个SNP几乎等价）
# r² > 0.2 → 有意义的LD
# r² < 0.1 → 基本独立
参考资料：PopLDdecay (Zhang et al. 2019, Bioinformatics)、PLINK 1.9 LD文档、GWLD R包 (2023)、Speciation Genomics Tutorial