从达尔文到代码手把手用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_pop2, n_samples100, n_snps1000): 生成两个群体的模拟基因型数据 参数 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, sizen_snps) # 群体1可能包含选择信号 pop1 np.random.binomial(2, neutral_freq, size(n_samples, n_snps)) # 群体2在部分SNP上产生分化 selected_idx np.random.choice(n_snps, size50, replaceFalse) 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_components2) 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_size50000, step10000): 滑动窗口计算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, rho0.001, max_snps200): 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, threshold5): 绘制选择信号图谱 参数 chrom: 染色体编号 positions: SNP物理位置数组 scores: XP-CLR得分数组 threshold: 显著性阈值 plt.figure(figsize(12, 4)) plt.scatter(positions, scores, s5, csteelblue) plt.axhline(ythreshold, colorr, linestyle--) plt.title(fXP-CLR Scores on Chromosome {chrom}) plt.xlabel(Physical Position (bp)) plt.ylabel(XP-CLR Score) plt.grid(alpha0.3) # 标注候选区域 candidate_idx np.where(scores threshold)[0] for idx in candidate_idx: plt.annotate(f{positions[idx]:,}, (positions[idx], scores[idx]), textcoordsoffset points, xytext(0,5), hacenter)5. 实战案例人类群体选择信号分析5.1 乳糖耐受性的遗传印记我们模拟分析欧洲人群和东亚人群在LCT基因区域的遗传分化# 模拟LCT区域数据 np.random.seed(42) chrom 2 start, end 136000000, 136500000 # LCT基因区域 positions np.random.randint(start, end, size500) genetic_pos (positions - start) / 1e6 * 1.2 # 转换为cM # 模拟选择信号 p_europe np.random.beta(0.8, 0.8, size500) p_europe[200:250] np.random.beta(5, 1, size50) # 模拟选择区域 p_eastasia np.random.beta(0.8, 0.8, size500) # 计算XP-CLR scores xpclr_core(p_europe, p_eastasia, genetic_pos) # 可视化 plot_selection_signals(chrom, positions, scores, threshold8)5.2 高原适应的遗传基础分析藏族人群与汉族人群在EPAS1基因区域的Fst分化# 模拟EPAS1区域数据 chrom 2 start, end 46500000, 47000000 # EPAS1基因区域 positions np.random.randint(start, end, size800) 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, markersize3) plt.axhline(y0.15, colorr, linestyle--) plt.title(fFst Values on Chromosome {chrom}) plt.xlabel(Physical Position (bp)) plt.ylabel(Fst) plt.grid(alpha0.3)6. 方法比较与进阶技巧6.1 Fst与XP-CLR的对比特征FstXP-CLR检测信号类型群体分化正选择时间敏感性反映长期分化检测近期选择计算复杂度较低较高结果解释分化程度量化选择强度量化适用场景群体结构分析适应性进化研究6.2 提高分析可靠性的策略多重检验校正from statsmodels.stats.multitest import fdrcorrection is_sig, adj_p fdrcorrection(p_values, alpha0.05)群体结构控制使用PCA或混合模型校正群体分层在相对同质的亚群体内进行分析参数优化建议XP-CLR的窗口大小通常设为50-100kb重组率rho根据物种调整人类≈1e-8/bp/generationFst窗口包含至少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_sizeint(size)) results[int(size)] { mean_fst: np.mean(fst), var_fst: np.var(fst), top5: np.percentile(fst, 95) } return results在西藏牦牛群体遗传分析的实际项目中我们通过这种从公式到代码的实现方式不仅验证了已知的高原适应基因还发现了一个新的候选基因区域。这种透明化的分析方法让研究者能够深入理解每个计算步骤的生物学意义而不是仅仅得到一个显著/不显著的二元结论。
从达尔文到代码：手把手用Python复现群体遗传学经典分析（XP-CLR/Fst计算实战）
发布时间：2026/6/2 0:35:02
从达尔文到代码手把手用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_pop2, n_samples100, n_snps1000): 生成两个群体的模拟基因型数据 参数 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, sizen_snps) # 群体1可能包含选择信号 pop1 np.random.binomial(2, neutral_freq, size(n_samples, n_snps)) # 群体2在部分SNP上产生分化 selected_idx np.random.choice(n_snps, size50, replaceFalse) 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_components2) 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_size50000, step10000): 滑动窗口计算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, rho0.001, max_snps200): 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, threshold5): 绘制选择信号图谱 参数 chrom: 染色体编号 positions: SNP物理位置数组 scores: XP-CLR得分数组 threshold: 显著性阈值 plt.figure(figsize(12, 4)) plt.scatter(positions, scores, s5, csteelblue) plt.axhline(ythreshold, colorr, linestyle--) plt.title(fXP-CLR Scores on Chromosome {chrom}) plt.xlabel(Physical Position (bp)) plt.ylabel(XP-CLR Score) plt.grid(alpha0.3) # 标注候选区域 candidate_idx np.where(scores threshold)[0] for idx in candidate_idx: plt.annotate(f{positions[idx]:,}, (positions[idx], scores[idx]), textcoordsoffset points, xytext(0,5), hacenter)5. 实战案例人类群体选择信号分析5.1 乳糖耐受性的遗传印记我们模拟分析欧洲人群和东亚人群在LCT基因区域的遗传分化# 模拟LCT区域数据 np.random.seed(42) chrom 2 start, end 136000000, 136500000 # LCT基因区域 positions np.random.randint(start, end, size500) genetic_pos (positions - start) / 1e6 * 1.2 # 转换为cM # 模拟选择信号 p_europe np.random.beta(0.8, 0.8, size500) p_europe[200:250] np.random.beta(5, 1, size50) # 模拟选择区域 p_eastasia np.random.beta(0.8, 0.8, size500) # 计算XP-CLR scores xpclr_core(p_europe, p_eastasia, genetic_pos) # 可视化 plot_selection_signals(chrom, positions, scores, threshold8)5.2 高原适应的遗传基础分析藏族人群与汉族人群在EPAS1基因区域的Fst分化# 模拟EPAS1区域数据 chrom 2 start, end 46500000, 47000000 # EPAS1基因区域 positions np.random.randint(start, end, size800) 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, markersize3) plt.axhline(y0.15, colorr, linestyle--) plt.title(fFst Values on Chromosome {chrom}) plt.xlabel(Physical Position (bp)) plt.ylabel(Fst) plt.grid(alpha0.3)6. 方法比较与进阶技巧6.1 Fst与XP-CLR的对比特征FstXP-CLR检测信号类型群体分化正选择时间敏感性反映长期分化检测近期选择计算复杂度较低较高结果解释分化程度量化选择强度量化适用场景群体结构分析适应性进化研究6.2 提高分析可靠性的策略多重检验校正from statsmodels.stats.multitest import fdrcorrection is_sig, adj_p fdrcorrection(p_values, alpha0.05)群体结构控制使用PCA或混合模型校正群体分层在相对同质的亚群体内进行分析参数优化建议XP-CLR的窗口大小通常设为50-100kb重组率rho根据物种调整人类≈1e-8/bp/generationFst窗口包含至少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_sizeint(size)) results[int(size)] { mean_fst: np.mean(fst), var_fst: np.var(fst), top5: np.percentile(fst, 95) } return results在西藏牦牛群体遗传分析的实际项目中我们通过这种从公式到代码的实现方式不仅验证了已知的高原适应基因还发现了一个新的候选基因区域。这种透明化的分析方法让研究者能够深入理解每个计算步骤的生物学意义而不是仅仅得到一个显著/不显著的二元结论。
