从达尔文到代码：手把手用Python复现群体遗传学经典分析（XP-CLR/Fst计算实战）
从达尔文到代码：手把手用Python复现群体遗传学经典分析（XP-CLR/Fst计算实战）
群体遗传学作为连接微观基因与宏观演化的桥梁，其分析方法正经历着从"黑箱工具"到"透明算法"的范式转变。本文将带您穿越150年的理论积淀，用Python代码重现XP-CLR和Fst这两个标志性分析方法的完整实现过程。不同于直接调用现成软件包，我们将从数学公式推导开始，逐步构建可解释、可定制的分析流程，让您真正掌握这些方法的计算本质。
1. 理论基础：从自然选择到统计量
1.1 群体遗传学的数学语言
群体遗传学的核心是将达尔文的自然选择、遗传漂变等概念量化为可计算的统计指标。以Fst为例，其本质是衡量群体间分化程度的方差分析：
Fst = (HT - HS)/HT其中HT代表总群体的基因多样性，HS代表亚群体内部的平均基因多样性。这个看似简单的公式背后蕴含着丰富的生物学意义：
- 0 ≤ Fst ≤ 1：完全无分化到完全隔离的连续谱系
- 阈值经验值：
- 0-0.05：弱分化
- 0.05-0.15：中等分化
0.15：强分化
1.2 XP-CLR的选择信号检测原理
XP-CLR方法通过比较两个群体的等位基因频率差异来检测正选择区域，其核心是构建复合似然比：
def xpclr_score(p1, p2, snp_positions): """计算XP-CLR得分的基本框架""" scores = [] for i in range(len(p1)): # 计算单个SNP的似然比 lr = likelihood_ratio(p1[i], p2[i]) scores.append(lr) return smooth_scores(scores, snp_positions)该方法特别擅长检测不完全的选择清除（incomplete sweep），即有利突变尚未在群体中达到固定的情况。
2. 数据准备与预处理
2.1 模拟数据的生成
为验证算法准确性，我们首先构建模拟群体数据：
import numpy as np from scipy.stats import beta def simulate_genotypes(n_pop=2, n_samples=100, n_snps=1000): """ 生成两个群体的模拟基因型数据 参数： n_pop: 群体数量 n_samples: 每个群体的样本数 n_snps: SNP数量 返回： genotypes: 形状为(n_pop, n_samples, n_snps)的数组 pop_labels: 群体标签 """ # 中性SNP的频率分布 neutral_freq = np.random.beta(0.5, 0.5, size=n_snps) # 群体1（可能包含选择信号） pop1 = np.random.binomial(2, neutral_freq, size=(n_samples, n_snps)) # 群体2（在部分SNP上产生分化） selected_idx = np.random.choice(n_snps, size=50, replace=False) shifted_freq = neutral_freq.copy() shifted_freq[selected_idx] = beta.ppf(0.9, 0.5, 0.5) pop2 = np.random.binomial(2, shifted_freq, size=(n_samples, n_snps)) return np.stack([pop1, pop2]), selected_idx2.2 真实数据的处理流程
对于WGS（全基因组测序）数据，标准的预处理步骤包括：
VCF文件解析：
import allel vcf = allel.read_vcf('input.vcf.gz', fields=['samples', 'variants/CHROM', 'variants/POS', 'calldata/GT']) gt = allel.GenotypeArray(vcf['calldata/GT'])质量控制：
- 缺失率过滤（<10%）
- 次要等位基因频率（MAF > 0.05）
- 哈迪-温伯格平衡检验（p > 0.001）
群体分层检测（避免假阳性）：
from sklearn.decomposition import PCA pca = PCA(n_components=2) coords = pca.fit_transform(gt.to_n_alt())
3. Fst计算的Python实现
3.1 单SNP的Fst计算
基于Weir & Cockerham的θ估计方法：
def calculate_fst(gt_pop1, gt_pop2): """ 计算两个群体间的Fst值 参数： gt_pop1: 群体1的基因型数组 (n_samples, ) gt_pop2: 群体2的基因型数组 (n_samples, ) 返回： fst: 该SNP的Fst值 """ # 计算等位基因频率 p1 = np.mean(gt_pop1) / 2 p2 = np.mean(gt_pop2) / 2 p_total = (np.sum(gt_pop1) + np.sum(gt_pop2)) / (2*(len(gt_pop1)+len(gt_pop2))) # 计算方差分量 n1, n2 = len(gt_pop1), len(gt_pop2) n_total = n1 + n2 hs = (n1*p1*(1-p1) + n2*p2*(1-p2)) / (n_total - 2) ht = 2*p_total*(1-p_total) # 避免除零错误 if ht == 0: return 0.0 return (ht - hs) / ht3.2 滑动窗口Fst分析
基因组水平的Fst通常采用滑动窗口计算：
def sliding_window_fst(pos, gt_pop1, gt_pop2, window_size=50000, step=10000): """ 滑动窗口计算Fst 参数： pos: SNP位置数组 gt_pop1: 群体1基因型矩阵 (n_samples, n_snps) gt_pop2: 群体2基因型矩阵 (n_samples, n_snps) window_size: 窗口大小（bp） step: 步长（bp） 返回： windows: 窗口位置列表 [(start, end), ...] fst_values: 各窗口Fst值 """ min_pos, max_pos = np.min(pos), np.max(pos) windows = [] fst_values = [] for start in range(int(min_pos), int(max_pos), step): end = start + window_size if end > max_pos: continue # 获取窗口内SNP索引 window_idx = np.where((pos >= start) & (pos <= end))[0] if len(window_idx) < 5: # 至少需要5个SNP continue # 计算窗口平均Fst window_fst = [] for idx in window_idx: fst = calculate_fst(gt_pop1[:, idx], gt_pop2[:, idx]) window_fst.append(fst) windows.append((start, end)) fst_values.append(np.nanmean(window_fst)) return windows, fst_values注意：实际应用中建议使用重叠窗口和加权平均，以平滑结果并减少信息丢失。
4. XP-CLR算法的完整实现
4.1 核心计算模块
def xpclr_core(p1, p2, genetic_pos, rho=0.001, max_snps=200): """ XP-CLR核心计算 参数： p1: 群体1的等位基因频率数组 p2: 群体2的等位基因频率数组 genetic_pos: 遗传距离数组（cM） rho: 重组率（每cM每代的重组事件） max_snps: 最大考虑的SNP数量 返回： scores: XP-CLR得分数组 """ scores = np.zeros(len(p1)) for i in range(len(p1)): # 选择邻近的SNP dist = np.abs(genetic_pos - genetic_pos[i]) neighbor_idx = np.argsort(dist)[:max_snps] # 计算复合似然比 lr = 0 for j in neighbor_idx: d = dist[j] * rho # 转换为遗传距离 weight = np.exp(-d) lr += weight * single_snp_lr(p1[j], p2[j]) scores[i] = lr / np.sum(np.exp(-dist[neighbor_idx]*rho)) return scores def single_snp_lr(p1, p2): """计算单个SNP的似然比""" # 中性模型下的联合概率 p_neutral = (p1 + p2) / 2 l_neutral = (p1**2 * p2**2 + (1-p1)**2 * (1-p2)**2) # 选择模型下的联合概率 l_selected = (p1 * p2 + (1-p1)*(1-p2)) return np.log(l_selected / l_neutral)4.2 结果可视化与解释
import matplotlib.pyplot as plt def plot_selection_signals(chrom, positions, scores, threshold=5): """ 绘制选择信号图谱 参数： chrom: 染色体编号 positions: SNP物理位置数组 scores: XP-CLR得分数组 threshold: 显著性阈值 """ plt.figure(figsize=(12, 4)) plt.scatter(positions, scores, s=5, c='steelblue') plt.axhline(y=threshold, color='r', linestyle='--') plt.title(f'XP-CLR Scores on Chromosome {chrom}') plt.xlabel('Physical Position (bp)') plt.ylabel('XP-CLR Score') plt.grid(alpha=0.3) # 标注候选区域 candidate_idx = np.where(scores > threshold)[0] for idx in candidate_idx: plt.annotate(f'{positions[idx]:,}', (positions[idx], scores[idx]), textcoords="offset points", xytext=(0,5), ha='center')5. 实战案例：人类群体选择信号分析
5.1 乳糖耐受性的遗传印记
我们模拟分析欧洲人群和东亚人群在LCT基因区域的遗传分化：
# 模拟LCT区域数据 np.random.seed(42) chrom = 2 start, end = 136000000, 136500000 # LCT基因区域 positions = np.random.randint(start, end, size=500) genetic_pos = (positions - start) / 1e6 * 1.2 # 转换为cM # 模拟选择信号 p_europe = np.random.beta(0.8, 0.8, size=500) p_europe[200:250] = np.random.beta(5, 1, size=50) # 模拟选择区域 p_eastasia = np.random.beta(0.8, 0.8, size=500) # 计算XP-CLR scores = xpclr_core(p_europe, p_eastasia, genetic_pos) # 可视化 plot_selection_signals(chrom, positions, scores, threshold=8)5.2 高原适应的遗传基础
分析藏族人群与汉族人群在EPAS1基因区域的Fst分化：
# 模拟EPAS1区域数据 chrom = 2 start, end = 46500000, 47000000 # EPAS1基因区域 positions = np.random.randint(start, end, size=800) genetic_pos = (positions - start) / 1e6 * 1.0 # 模拟基因型数据 gt_han = np.random.binomial(2, 0.2, size=(50, 800)) gt_tibetan = gt_han.copy() gt_tibetan[:, 350:400] = np.random.binomial(2, 0.8, size=(50, 50)) # 计算滑动窗口Fst windows, fst_values = sliding_window_fst(positions, gt_han, gt_tibetan) # 可视化 plt.figure(figsize=(12,4)) window_centers = [(w[0]+w[1])/2 for w in windows] plt.plot(window_centers, fst_values, '-o', markersize=3) plt.axhline(y=0.15, color='r', linestyle='--') plt.title(f'Fst Values on Chromosome {chrom}') plt.xlabel('Physical Position (bp)') plt.ylabel('Fst') plt.grid(alpha=0.3)6. 方法比较与进阶技巧
6.1 Fst与XP-CLR的对比
| 特征 | Fst | XP-CLR |
|---|---|---|
| 检测信号类型 | 群体分化 | 正选择 |
| 时间敏感性 | 反映长期分化 | 检测近期选择 |
| 计算复杂度 | 较低 | 较高 |
| 结果解释 | 分化程度量化 | 选择强度量化 |
| 适用场景 | 群体结构分析 | 适应性进化研究 |
6.2 提高分析可靠性的策略
多重检验校正：
from statsmodels.stats.multitest import fdrcorrection is_sig, adj_p = fdrcorrection(p_values, alpha=0.05)群体结构控制：
- 使用PCA或混合模型校正群体分层
- 在相对同质的亚群体内进行分析
参数优化建议：
- XP-CLR的窗口大小通常设为50-100kb
- 重组率（rho）根据物种调整（人类≈1e-8/bp/generation）
- Fst窗口包含至少20个SNP以保证稳定性
def optimize_window_size(gt_pop1, gt_pop2, pos, sizes=[20e3, 50e3, 100e3]): """测试不同窗口大小对结果的影响""" results = {} for size in sizes: windows, fst = sliding_window_fst(pos, gt_pop1, gt_pop2, window_size=int(size)) results[int(size)] = { 'mean_fst': np.mean(fst), 'var_fst': np.var(fst), 'top5': np.percentile(fst, 95) } return results在西藏牦牛群体遗传分析的实际项目中，我们通过这种从公式到代码的实现方式，不仅验证了已知的高原适应基因，还发现了一个新的候选基因区域。这种透明化的分析方法让研究者能够深入理解每个计算步骤的生物学意义，而不是仅仅得到一个"显著/不显著"的二元结论。