做物种GWAS分析的朋友大概率都遇到过同一个尴尬局面手里握着一套猪、鸡、牛或者其他农业动物的GWAS summary statistics想估算全基因组水平的SNP遗传力或者看看不同性状之间的遗传相关性结果打开标准的ldsc流程发现整个人类遗传学生态是围绕hg19/hg38、HapMap3 SNP集合、1000 Genomes参考面板搭起来的直接拿来跑非人类数据轻则结果离谱重则中途报错。今天这篇就专门拆一拆ldsc跨物种计算这件事把原理、流程、参考面板怎么选、哪些坑最容易踩一次说清楚。ldscLD Score Regression本身是一种回归方法通过对GWAS统计量在LD Score上的回归把多基因信号和混杂因素分开从而估计遗传力和遗传相关性。跨物种计算的关键变化在于三个变量参考基因组版本、LD参考面板、SNP注释信息。这三个变量一旦对不上号整个结果的生物学含义就全变了。如果你是做动物遗传育种、进化遗传学或者比较基因组学研究的这篇博文的目标非常明确让你理解为什么ldsc跨物种计算不能粗暴套用人管线和参考面板同时给你一套能够真正跑通且结果可信的完整操作框架。1. Ldsc的统计学内核为什么它能跨出人类领域在动手操作之前有必要把ldsc的底层逻辑讲透因为跨物种计算的每一处细节调整本质上都是在维护这个统计假设成立的前提。LD Score Regression的核心思路是对每一个SNP先计算它的LD Score即该SNP与周围一定范围内所有其他SNP的平方相关系数r²之和衡量这个SNP被周围位点代理的程度然后把GWAS的检验统计量χ²当作因变量LD Score当作自变量做回归。这个回归背后的生物学逻辑是一个SNP的LD Score越高说明它和周围更多位点存在连锁不平衡如果某个性状确实受大量因果位点控制那么在高LD区域的SNP就更有可能因为搭便车效应而表现出显著的关联信号。因此真实多基因遗传力会让回归斜率整体抬高而群体分层、隐性关联、样本重叠等混杂因素造成的信号则不会随LD Score成比例变化这些混杂统一体现在回归截距上。数学上最核心的期望公式可以简化为E[χ²_j] 1 N·h²_g·ℓ_j / M N·ρ其中N是GWAS样本量h²_g是SNP遗传力ℓ_j是第j个SNP的LD ScoreM是参与分析的SNP总数ρ是混杂造成的偏差项。从这个公式能直接推出几个跨物种计算的关键结论LD Score必须来自与GWAS样本相近的群体。因为不同物种、甚至同一物种不同亚群的LD模式差异巨大用人类的LD Score去回归猪的GWAS统计量相当于拿北京的房价走势去预测鹤岗的房价根本不在一个数据生成机制里。回归斜率代表真实遗传信号截距代表混杂信号。跨物种分析时如果你发现遗传力估计值出现负值第一个要怀疑的就是LD Score与真实LD模式不匹配而不是遗传力真的是负的遗传力在数学上确实可能出现负估计但通常暗示模型设定有问题。MSNP总数的选择会影响绝对数值。人类遗传学使用HapMap3 SNP集合作为M是因为这些SNP在芯片和参考面板之间覆盖较好能降低LD Score估计的噪声。跨物种计算时如果随便取全部SNPM的定义不同遗传力估计值自然不可比。理解了这些就明白跨物种ldsc不是换一个--bfile路径那么简单而是要重新搭建一整套与目标物种匹配的LD参考体系。2. 跨物种计算的第一个门槛参考基因组与SNP注释的物种一致性这是最容易被忽视、一旦出错就满盘皆输的环节。人类GWAS的summary statistics通常直接用rsID标识SNP但在猪、鸡、牛、羊、鱼等非模式物种中情况完全不同。2.1 SNP标识符rsID不是万能的人类dbSNP数据库拥有丰富的rsID注释绝大多数GWAS位点都能通过rsID直接匹配到参考面板。但农业动物的dbSNP覆盖度参差不齐许多品种的特异性变异压根没有rsID。你拿到的非人类GWAS数据SNP标识通常只有色谱位置chr:pos甚至只有内部设计的SNP编号。这里没有捷径唯一可靠的做法是自己构建比对表以目标物种参考基因组为基准将GWAS中的chr:pos与基因型参考面板的chr:pos逐一匹配。具体步骤大致如下统一参考基因组版本。GWAS数据是ARS-UCSCv2.0跑的你的基因型参考面板却是Sscrofa11.1的坐标必须先把坐标转换统一。这一步最常用liftOver或者物种专属的坐标转换工具。以chr:pos为key做inner join。因为不同版本之间可能还有染色体编号命名差异比如Ensembl的1号染色体叫1而UCSC可能叫chr1先统一命名再匹配。仔细核对等位基因。2.2 等位基因链问题的处理等位基因的匹配是所有做跨物种分析的人必经的一道坎。非人类物种的芯片设计本身可能包含不同链上的探针加上基因组版本更新时区域发生重链翻转GWAS数据里记录的A1/A2与参考面板里的REF/ALT可能完全相反。对于这种情况常用策略是先根据chr:pos比对比对上了再看等位基因是否一致一致直接通过若等位基因是互补关系比如A/C变T/G确认是否由测序链不同导致若是则整列翻转若完全对不上等位基因这通常是坐标错位或者参考基因组版本不一致建议剔除而不是强行保留。在实际操作中我见过很多人在这一步偷懒直接用只比对位置、忽略等位基因的方式跑ldsc结果遗传力估计值高得离谱或者出现正的截距异常最后又花大量时间去排查数据反而更慢。2.3 参考基因组版本和物种常用参考面板整理一张常见的物种参考基因组与可用的LD参考来源表方便不同方向的读者直接对照物种常用参考基因组可参考的LD信息来源人GRCh37/hg19, GRCh38/hg381000 Genomes Phase 3, HapMap3猪Sscrofa11.1, ARS-UCSCv2.0自建基因型面板或PigGenes数据库牛ARS-UCD1.2, Btau_5.0.11000 Bull Genomes, 自建群体面板鸡GRCg6a, Gallus_gallus-5.0自建GBS/芯片面板绵羊Oar_rambouillet_v1.0, Ovis_aries_v3.1自建芯片面板斑马鱼GRCz11, danRer11自建群体基因型面板注意对于非模式物种几乎不存在像1000 Genomes那样权威且通用的LD参考面板。你手中如果正好有一套大样本量的基因型数据芯片或者全基因组测序结果那它就是最合适的LD参考来源。3. 实操流程从基因型数据到LD Score再到遗传力估计理论铺垫完之后来看一套完整的跨物种ldsc实战流程。这里的假设是你已经拥有一份目标物种的GWAS summary statistics、一份该物种的基因型数据PLINK格式。3.1 第一步数据准备与质控先处理基因型参考面板。ldsc计算LD Score需要PLINK格式的bed/bim/fam文件。基因型数据拿到手之后建议做下面几步预处理过滤低MAF位点。具体阈值可以根据物种和芯片密度调整一般建议MAF0.01因为极低频位点的LD Score估计极不稳定。去除高缺失率个体和SNP。标准是--mind 0.1和--geno 0.1左右。检查染色体命名一致性。bim文件里的染色体编号要和GWAS数据里的编号逻辑统一。如果基因型面板包含的SNP数量特别多比如全基因组测序数据可以按--ld-wind-snps参数控制窗口内的SNP数避免LD计算时间爆炸。GWAS summary statistics需要整理成ldsc要求的格式它要求包含以下列列名含义SNPSNP标识符A1效应等位基因A2另一个等位基因ZGWAS中的z-score如果没有z-score可以用beta/se计算z beta / seNGWAS样本量若各SNP样本量不同这一列建议用每个SNP实际的N实际操作中如果原始GWAS没有z-score而只有p值和beta/se没关系z beta/se直接生成。另外如果原始GWAS是逻辑回归类型同样适用。3.2 第二步使用物种自身的基因型数据生成LD Score命令行如下核心是--l2参数python ldsc.py \ --bfile species_genotype \ --l2 \ --ld-wind-snps 200 \ --out species_LDScore这里的--ld-wind-snps 200表示以捕获得最近的200个SNP为窗口计算LD Score是ldsc默认值。人类遗传学分析中默认窗口是1 cM约1 Mb但在跨物种场景下各物种重组率差异很大如果基因型面板密度较高建议根据注释信息估算合适的物理距离窗口。计算结束后会生成后缀为.l2.ldscore.gz的文件这就是后面回归要用到的LD Score数据。3.3 第三步运行遗传力回归有了LD Score遗传力估计就一行命令python ldsc.py \ --h2 species_gwas.sumstats.gz \ --ref-ld-chr species_LDScore \ --w-ld-chr species_LDScore \ --out species_h2这里有一个容易忽略的细节--ref-ld-chr和--w-ld-chr在标准人类流程中经常分别为baseline LD和1000 Genomes权重LD是因为人类场景下回归加权的权重与LD参考可能需要不同的SNP集合。在跨物种自建参考面板时通常可以直接都用同一个LD Score文件。如果你发现回归结果极不稳定可以分别构造不同的参考与权重数据再对比。3.4 第四步检查回归诊断信息ldsc输出的log文件里会包含关键信息回归斜率、截距、遗传力估计值、观察到的平均χ²、LD Score均值等。其中要重点看的平均χ²。如果远大于1说明要么样本量非常大要么存在明显混杂需要关注截距的估计。截距。理想情况下截距接近1。如果截距明显大于1且置信区间不包含1说明存在群体分层或样本重叠等混杂遗传力估计需要谨慎解读。回归斜率。如果斜率不显著甚至为负请立即怀疑LD Score与GWAS数据的等位基因链匹配问题或者参考面板与GWAS样本群体的亲缘关系太远。4. 踩坑实录跨物种ldsc最容易翻车的五个细节这一节是全文最值钱的部分因为我个人在跨物种ldsc实战中把这些坑一个不落地全踩过一遍每条都对应一次通宵排错。4.1 窗口设置与物种重组率不匹配人类遗传学的标准ldsc通常使用1 cM窗口或等价的物理距离这个选择基于人类约1 cM/Mb的重组率。但不同物种的重组率差异可高达数倍鸡的常染色体总图距约3300 cM猪约2000 cM斑马鱼则截然不同。跨度设置不对LD Score所捕获的周围位点影响范围就错了回归斜率的基因解释也会失效。如果你的基因型面板有遗传图谱信息优先采用遗传距离定义窗口比如1 cM而不是物理距离。没有遗传图谱的情况下一个折衷方案是根据物种有效群体大小估算LD衰减距离然后再决定以多少个SNP或多大物理距离作为窗口。4.2 MAF过滤阈值不一致导致LD Score偏差LD Score的计算高度依赖群体等位基因频率。参考面板中的SNP如果MAF偏低r²值通常估计不准确抓入GWAS分析时MAF分布不匹配同样会引入系统偏差。在实际操作中我的建议是参考面板和GWAS summary statistics的MAF过滤原则要保持一致至少在数量级上不要差太远。如果GWAS数据来自低深度测序其中有大量低频变异芯片又没有覆盖那你计算出的LD Score天然就和GWAS统计量的频率范围不匹配遗传力估计会出现向下偏倚。4.3 样本量差异和截距的解读跨物种也一样适用许多跨物种GWAS实际来自多机构合并的meta分析不同SNP的N列时常不同。在某些数据中因为某些芯片位点在部分群体里缺失N缺得厉害。这时如果直接给全部SNP设一个固定总样本量回归截距的计算就会失真。ldsc的正确做法是用每个SNP实际的N在sumstats中写入样本量列。另一个情景是GWAS数据和参考面板来自同一个牛场、同一个鸡群这时样本重叠本身不是问题因为LD参考面板用于刻画LD结构混杂估计完全依赖截距判断。但如果你的meta分析中部分机构和参考面板的数据同源截距偏高不一定就是群体分层也可能是参考样本本身混进了GWAS样本需要在写论文时把这种情况说清楚。4.4 二值性状的责任尺度转换公式跨物种遗传力估计很多情况下处理的是抗病力、发病率等二值性状。ldsc对二值性状的默认遗传力是观测尺度observed scale但大家实际关心的往往是责任尺度liability scale遗传力。跨物种情况下群体发病率的数据经常缺失、不准确更容易忘记转换。责任尺度转换公式在ldsc输出结果中可以直接用--h2-ldsc参数输出的观测尺度结果乘以一个转换因子得到责任尺度。转换因子P(1-P) / [z(t)²]其中P是群体患病率t是正态分布阈值z(t)是概率密度函数在阈值处的取值。如果群体的P估计不可靠不要强行转换宁可报告观测尺度结果并在讨论中说明。4.5 注释信息与分层LD Score的扩展问题人类遗传学里的分层LD Score回归partitioned heritability把SNP分为功能类别如编码区、增强子、启动子等再估计每个类别的遗传力富集度。跨物种计算中如果你也想做这一步最大的障碍是功能注释体系的物种差异。人类的功能注释来自ENCODE、Roadmap等丰富的数据而非模式物种往往只有参考基因组上的基因注释甚至只有转录组数据。一个可行的替代方案是基于参考基因组的GENCODE或Ensembl注释自己构造基因区、基因间区、UTR区等简化的注释文件。与人类分析相比注释类别要少很多但依然能提供有价值的富集信号。python ldsc.py \ --bfile species_genotype \ --l2 \ --annot species_gene_annot \ --ld-wind-snps 200 \ --out species_LDScore_annot这里需要生成annot格式文件每一行对应一个SNP每一列对应一个注释类别值表示该SNP是否属于该类别1/0。然后运行python ldsc.py \ --h2 species_gwas.sumstats.gz \ --ref-ld-chr species_LDScore_annot \ --w-ld-chr species_LDScore_annot \ --overlap-annot \ --print-coefficients \ --out species_h2_partitioned跨物种做分层回归时由于注释类别少、相互重叠有限遗传力富集度的估计会比较粗糙但仍然可以从编码区vs非编码区谁解释了更多遗传力的角度得出有意义的生物学结论。5. 低于人类复杂度的场景遗传相关性和跨物种比较最后再聊一个常常被问到的应用场景。除了遗传力估计ldsc还能估算两个性状之间的遗传相关性。在跨物种研究中这有可能服务于不同目标同一物种不同性状的遗传相关性。这不需要跨物种计算但要注意summary statistics同样需要自建LD Score。不同物种同一性状的遗传相关性。这种情况本质上非常复杂因为不同物种的SNP集合、等位基因频率、LD结构都不一样很难找到一一对应的正交SNP集合。除非是做全基因组序列比对级别的分析否则不推荐直接跑ldsc --rg跨物种比较。在跨物种计算的世界里一个相对实际的应用是当你没有现成的人类那样的权威LD参考面板但有多个品种或品系的基因型数据时可以分别用不同品系做LD参考跑同一套GWAS数据然后对比结果。这种设计可以评估LD参考面板差异对遗传力估计的影响幅度也是审稿人常问的问题之一。下面是一张我实测中不同参考面板对同一批模拟GWAS数据的遗传力估计影响表数值为演示值参考面板来源h²估计值截距平均χ²备注目标群体自身A群体0.321.021.21基准结果同物种另一群体B群体0.281.051.22LD模式差异导致低估人类参考面板错误示范0.730.611.24结果不可信这个表格结果符合理论预期参考面板与GWAS样本的群体亲缘关系越远LD Score与真实LD模式越不匹配回归结果越不稳定甚至截距会偏离1。如果你确实需要做跨物种的遗传结构比较更稳妥的方案是先在每个物种内单独估算遗传力再比较h²估计值的相对大小而不是直接做一个跨物种的rg。毕竟遗传力本质上是一个物种内部、特定群体、特定环境下的描述性参数跨物种直接比较绝对值存在大量不可控的变数。6. 最后再分享一点个人经验做跨物种ldsc分析我一直觉得最核心的教训是不要试图用人类的工具链生态去硬套非人类数据而是理解工具背后的统计假设再根据物种特性去调整数据输入。另外在跨物种实战中如果结果异常但找不出原因我习惯先跑单条染色体做小规模测试——选一条中等大小的常染色体用该染色体上的SNP子集独立计算LD Score并回归。如果单条染色体的结果跟全基因组结果趋势一致那说明整体流程无大毛病如果结果方向相反就赶紧回头查参考面板和GWAS之间的样本对应关系。这种做法的好处是单条染色体跑得非常快迭代调试的效率远高于每次全基因组重跑。最后准备一个快速检查清单每次跑完跨物种ldsc都按这个自查一遍GWAS数据与参考面板的参考基因组版本是否一致染色体编号命名是否统一chr前缀有没有处理干净等位基因链方向是否匹配互补链是否已处理参考面板的MAF过滤阈值是否和GWAS的MAF范围在同一量级LD Score窗口选择是否与物种重组率匹配回归结果截距是否接近1斜率是否为正summary statistics中的N列是不是按SNP实际样本量填写。这些细节全部对上了跨物种ldsc的结果才值得被写进论文里。否则花了几个月整理数据最后因为一个等位基因链方向的小问题得出完全错误的遗传力结论那才叫真的欲哭无泪。