第一次跑LDSC跨物种计算的时候我以为就是把人的GWAS数据换成一个物种的summary statistics然后按老流程走一遍而已。真正上手才发现这个分析的本质根本不是“换个输入文件”而是要把两套完全不同坐标系里的遗传信号投影到同一个框架里比较。做人和小鼠的跨物种LDSC最核心的一件事是你手里的小鼠GWAS结果落在小鼠基因组坐标上而LDSC的参考模型、LD分数、基因注释全部是围绕人类基因组构建的这两套数据不能直接合并。你需要通过同源信息做坐标映射再针对映射后的SNP集合重新计算LD Score之后才能做遗传力估计和跨物种遗传相关分析。这篇文章会把这套流程的整体设计、每一步的细节和坑都拆开讲。主要面向手里已经有人或动物GWAS数据、想用LDSC做跨物种遗传结构比较的研究者也适合刚接触LDSC但一直没搞明白“为什么跨物种不是简单换数据”的人。1. LDSC跨物种计算到底在算什么1.1 不只换物种是换坐标和模型先简单回顾一下LDSC的核心逻辑因为跨物种场景里的很多坑都源于对模型假设的误解。LDSC全称是LD Score Regression它的出发点是一个非常简单的关系对某一个SNP j它的GWAS检验统计量 χ²_j 的期望值是E[χ²_j] ≈ 1 N·h²·ℓ_j/M N·a其中N是GWAS样本量h²是所有SNP贡献的总遗传力M是参与分析的SNP总数ℓ_j是这个SNP的LD Score定义为它和周围连锁不平衡区域内所有SNP的 r² 之和a则代表群体分层或隐匿相关性带来的混淆因子。这个模型能够成立依赖一个关键前提我们用来计算LD Score的参考面板必须和GWAS效应量的坐标体系一致。这里的“坐标一致”不仅仅是基因组位置一致还包括等位基因方向一致、SNP ID一致、LD结构一致。放到跨物种场景下问题立刻来了。假设人GWAS里有SNP rs123456位于人类chr1:1000000小鼠对应的同源区域在那段基因的3端但小鼠芯片和测序分型的位点可能是另一个rsID、另一个等位基因组合甚至位置坐标在小鼠基因组的chr4上。你如果不做任何处理直接把人和小鼠的SNP合并跑LDSC模型里的ℓ_j是按人类LD结构算出来的而小鼠GWAS的效应量却来自小鼠群体的LD结构两者根本不匹配回归自然失去意义。所以“跨物种计算”真正的含义是利用人类参考模型作为桥梁通过同源映射把另一种物种的GWAS信号“翻译”到人类坐标体系然后在统一的注释和LD框架下分析。1.2 跨物种计算的三条常见路线我实际接触到的跨物种LDSC分析基本上可以归成三大类每类的目的和数据要求差别很大。第一类叫保守性富集分析。思路是先定义一组跨物种保守区域比如人类和小鼠的同源编码区、保守调控元件、受纯化选择的区间然后用分层LDSCstratified LDSC简称S-LDSC看这些区间里的遗传力富集程度。如果人GWAS的遗传力在保守区域显著富集说明这些区域的变异对性状有不成比例的影响再在小鼠GWAS上做同样的分析如果富集模式一致就可以说“该性状的遗传结构在物种间是保守的”。这种分析不一定算跨物种遗传相关但对理解保守性非常重要。第二类是单物种遗传力估计的对照分析。比如你在小鼠身上做了某项行为学实验的GWAS样本量不大但你想知道这些位点的贡献是否集中在某些功能类别里。你可以把小鼠GWAS映射到人类坐标后跑S-LDSC比较不同注释类别上的富集差异。这类分析个体户也能跑通只要映射这一步做扎实。第三类是跨物种遗传相关分析也就是利用cross-trait LDSC的模式计算人类性状和小鼠性状之间的遗传相关r_g。这个分析在思路上和普通cross-trait LDSC没有本质区别仍然是两组GWAS summary statistics加上一个reference panel通过协方差和方差的比值得到r_g。但特殊之处在于两组GWAS的SNP集合必须映射到同一个参考坐标下否则跨物种r_g就是无源之水。三条路线可以单独做也可以组合先用路线一判断保守性再用路线三量化相关强度。我这次主要写的是路线三的完整流程因为它的步骤最全、坑也最多。2. 数据准备映射这一步决定成败2.1 同源基因映射别直接按rsID合并我第一次踩坑就是在这一步想省事直接用小鼠GWAS的rsID去合并人类LD Score文件结果发现合并率不到20%。原因其实很简单——小鼠的rsID是人类dbSNP里的ID同一段序列在小鼠基因组上多数没有对应的rsID或者因为芯片设计不同分型的是另一个位点。正确做法是基于同源基因做坐标映射。核心思路是把小鼠GWAS的每个SNP注释到它所在的小鼠基因上用同源映射表找到这个基因对应的人类基因把SNP位置“平移”到该人类基因的转录本坐标上得到它在人类基因组中的位置。具体实现时我个人推荐直接用UCSC的liftOver工具配合链文件chain file把小鼠基因组坐标如mm10直接转换到人类基因组坐标如hg19转换后再根据位置匹配到人类LD Score文件的SNP。这样可以避免手动做基因映射时的外显子/内含子边界偏差。不过liftOver并不是万能的。它本质上做的是基于全基因组比对坐标的转换对于同源基因、保守非编码区转换准确率很高但对于物种特异的插入/缺失区域、高度重复区域liftOver可能返回空值或错误位置。实操中我一般会设一个最低保留率如果映射后的SNP数量不到原始SNP数的60%我会考虑换一个链文件版本或改用基于同源基因表的映射方式。2.2 等位基因方向和strand问题跨物种映射里最隐蔽的坑是等位基因方向不一致。人在参考基因组上记录的是正链等位基因而小鼠芯片的探针设计可能是在负链上检测的两者的A/C等位基因在互补链上会颠倒。如果这一步不对齐后面LDSC回归算出来的遗传力会严重偏低遗传相关甚至可能变成负值。判断方向的方法其实不复杂。LDSC官方建议用等位基因频率做参考找一组映射后的SNP比较它们在小鼠GWAS里的effect allele频率和人类参考面板里的ref allele频率如果两者的等位基因定义一致那么频率应该是正相关如果大多是负相关说明整体strand搞反了需要把效应量和等位基因同时取互补。这一步我最常用的工具是LDSC官方仓库里的munge_sumstats.py它在清理summary statistics时会自动做allele alignment但前提是输入文件里必须包含A1、A2两列并且A1是效应等位基因A2是另一个等位基因。如果原始GWAS只给了effect和other没有明确哪个是A1、哪个是A2一定要先处理好再进munge。2.3 清洗summary statistics的实操规范和普通LDSC一样跨物种分析的第一步是整理summary statistics。这里特别强调几个针对跨物种场景的额外要求SNP ID列在munge之前先做好“物种坐标→人类坐标”的映射并为每个SNP指定人类rsID。不要试图在munge之后再做映射因为munge会过滤掉一部分SNP你事后合并不回原始坐标。样本量列不同物种GWAS的样本量差异可能很大小鼠GWAS可能只有几千个体人的却有几十万。LDSC回归里N是已知量直接写在summary statistics里没问题但要清楚N差异大时权重计算会受影响后面结果解释要谨慎。去除重复SNP映射到人类坐标后可能出现多个小鼠SNP对应同一个人类SNP位点的情况此时需要保留P值最小或样本量最大的那个否则LDSC会报错。清理完成后建议单独做一个QC报告映射率、保留率、等位基因方向一致性、有效SNP数量这些数字先记录下来后面分析出问题好排查。3. 完整实操人-小鼠跨物种LDSC分析流程3.1 环境与工具准备LDSC的运行环境不需要我多说主要是用Python2/3兼容的版本跑官方脚本配套的还有ldsc.py、munge_sumstats.py两个核心入口以及baseline-LD模型文件。做跨物种分析时我额外建议准备以下资源人类参考LD Score文件推荐使用官方提供的eur_w_ld_chr/欧洲人群参考或baselineLD_v2.2目录按你分析的性状人群选择人类-SNP位置注释文件用于把映射后的位置转换为rsID人-鼠同源映射表推荐使用MGI的HOM表HomoloGene也可以但MGI的注释信息更全链文件用于liftOver坐标转换hg19ToMm10或mm10ToHg19根据你的参考版本来选。工具版本上我要提醒一句LDSC官方脚本对Python版本比较敏感很多人在Python3.8环境里跑会报module相关的错。我自己一般用conda单独建一个Python3.6环境来跑LDSC省得和日常环境互相干扰。3.2 从映射坐标到LD Score计算流程上第一个正式步骤是坐标映射。假设你对小鼠GWAS已经完成了基本的格式整理比如每行是一个SNP包含chr、pos_mm10、A1、A2、Z或P。现在要做的是# 使用liftOver进行坐标转换 liftOver mouse_snps.bed mm10ToHg19.over.chain.gz mouse_snps_hg19.bed unmapped.txt这一步会得到一个小鼠SNP在人类基因组上的坐标区间。接下来我一般会用bedtools intersect或awk脚本把这些坐标和人类dbSNP的坐标文件做匹配拿到对应的rsID。这里有一个非常容易被忽略的点liftOver输出的坐标区间可能不止一个因为基因组比对在旁侧同源序列上可能发生连锁映射。我建议只保留第一个匹配且长度最短的区间并过滤掉映射后的坐标落在chrX/chrY等特殊染色体上的SNP。拿到人类rsID之后下一步是把这些SNP集合对应的LD Score单独提出来。这个操作看似简单但有个细节要注意官方LD Score文件里SNP数量远大于你映射得到的SNP数量如果你直接把ldsc.py里的--out设为全部人类SNP的LD Score文件回归时会因为多数SNP没有表型信息而被自动排除。正确做法是准备一个--merge-alleles文件或者直接用--ref-ld-chr指定一个只包含映射SNP的LD Score子集。为什么非要重新算LD Score原因我之前提过LDSC的模型是以LD Score作为自变量而LD Score本身是从参考面板里估计出来的。如果你分析的目标SNP集合只是全基因组的一小部分直接用全LD Score会引入一个偏差目标SNP集合内的经验LD Score分布和全基因组分布不同回归截距会漂移。最稳妥的方式是先用PLINK在参考面板上针对你映射后的SNP集合计算LD矩阵再按LDSC格式生成这个子集的LD Score文件。如果不想自己算也可以退而求其次直接用官方LD Score文件但通过--keep-snps参数指定SNP集合。这个参数在ldsc.py里是支持的它会保留你指定的SNP。实测下来只要保留的SNP数量在8万以上结果和重新算LD Score差别不大低于这个数量我强烈建议别省这一步。真实项目里小鼠GWAS往往覆盖5万到20万SNP映射到人类后能匹配到的有效SNP可能只有3万到8万这就比较危险了。3.3 单性状遗传力估计先验证再交叉在做跨物种遗传相关之前我强烈建议先把两个物种各自的遗传力都单独跑一遍。这一步有三个作用验证数据清洗是否正确、确认映射后的SNP集合是否足够支撑回归、为后续解释r_g提供参照。对每一组GWAS命令大致是python ldsc.py \ --h2 mousedata.sumstats.gz \ --ref-ld-chr eur_w_ld_chr/ \ --w-ld-chr eur_w_ld_chr/ \ --out mouse_h2_firstrun运行结束后观察输出文件里的几个关键指标截距Intercept应该接近1如果明显大于1说明GWAS存在人群分层或样本重叠遗传力估计值Total Observed h2应该是一个正数且标准误不要太大。如果h2是负数或者标准误比估计值还大一般说明有效SNP太少、样本量太小、或者等位基因方向没有对齐。跨物种场景里这一步经常出现“h2为0”的情况。我遇到十次有八次是因为munge的时候等位基因方向没对齐导致Z分数的符号和参考面板不匹配。还有一种情况是映射后的SNP太少模型拟合不出来。这时候不要急着跑交叉流程先回到上面的QC报告去查。只有两个物种的单性状估计都合理再继续下一步。3.4 cross-trait LDSC计算跨物种遗传相关交叉分析的本质是同时读入两组GWAS的summary statistics利用同一个reference面板估计两组效应量之间的协方差和各自的方差最后得到遗传相关r_g。命令结构如下python ldsc.py \ --rg human.sumstats.gz,mouse_hg19.sumstats.gz \ --ref-ld-chr eur_w_ld_chr/ \ --w-ld-chr eur_w_ld_chr/ \ --out cross_species_rg你需要确认两组文件都已经munge过并且都在同一套人类坐标体系下。很多人在这一步犯的错误是人的GWAS没有做任何跨物种映射用的是全基因组所有SNP而小鼠GWAS只映射了一部分SNP两组SNP集合差异巨大。LDSC处理这种不平衡集合时交叉项方差会被大量“缺失”的SNP稀释导致r_g被严重低估。解决办法是准备一个共同的SNP集合叫人鼠同源SNP集合就是映射后能同时出现在两个GWAS里的那些位点。在跑--rg前用--keep-snps把这个交集文件传进去强制LDSC只在共同位点上分析。这个方法实操效果很好但要注意交集SNP明显变少的时候标准误会变大结果解释要谨慎。输出文件里除了r_g还要注意看Genetic Covariance和p值。如果r_g为正但p值不显著可能是样本量不足不要强行解释成“有相关趋势”如果r_g为负且置信区间很宽先别急着下结论优先排查等位基因方向。4. 常见问题与排查技巧实录4.1 映射率太低该怎么办你最可能遇到的第一个问题就是映射率低。小鼠GWAS原始SNP数5万liftOver后只剩2万这种情况我见过不止一次。排查思路分成三路先看链文件版本是否匹配。小鼠参考基因组有mm9、mm10、mm39人类有hg19、hg38不同版本之间的链文件互不兼容使用错误的链文件会直接把大量坐标drop掉。再看原始数据的chr命名。有些GWAS输出是chr1有些是1在liftOver前统一成标准格式可以提升命中率。最后看重复区域。如果丢失的SNP主要在着丝粒、端粒等重复区说明坐标本身没有错只是这些区域没有可靠的跨物种比对这时可以接受。如果是基因层面的映射还有一种补救方式不直接做坐标转换而是用同源基因表把小鼠SNP归到基因再取对应人类基因上下游一定范围比如±50kb内的SNP作为匹配集。这个方法会引入更多噪音但在映射率极低时可以作为备选分析时配合--keep-snps限制位点集合就好。4.2 截距异常、h2估计为0或负值这是LDSC跨物种分析里最让人头疼的一类问题而且往往从头到尾都找不到明显报错。我总结出三个高频原因。第一个是等位基因方向没对齐。现象是h2趋近0且出现大量负的遗传力组件。解法是回到munge输出文件手动检查都50个映射SNP的A1和A2频率是否与参考面板一致如果一半以上反了整体取补再重新跑。第二个是样本量策略不对。小鼠GWAS总样本量N可能只有2000映射后有效SNP也少回归中N会作为权重的一部分过小的N会让所有SNP的统计量都偏低截距和斜率都无法估计。此时可以考虑在munge时采用--N-cas和--N-con如果是个案-对照设计把样本量结构明确写清楚或者直接换用更大的小鼠GWAS数据集。第三个是参考面板不匹配。我做过一个跨物种分析发现人GWAS是非洲人群但用的LD Score参考是欧洲人群结果截距和h2都偏高。跨物种场景里你能选的参考面板本来就有限基本原则是尽可能选择与小鼠近交系遗传背景相似的参考设计。当然LDSC的reference本身来自人类参考基因组这个环节没法替换只能在解释时说明局限性。4.3 跨物种r_g与预期不符的排查顺序假设你预期人与小鼠在某个性状上有比较高的正遗传相关结果r_g接近0甚至为负。我先说结论这个结果不一定是错的跨物种遗传结构本身就可能存在差异比如同样是一个行为性状人类和小鼠的调控机制可能已经分化到相关度很低。但如果要排查技术问题我建议按这个顺序来查共同SNP数量。计算--keep-snps后的SNP交集大小如果少于3万结果基本不可靠查基因型方向。对比两组GWAS的效应量符号是否与LD Score参考一致查GWAS质量的指标。比如两组GWAS的λGC是否异常高、截距是否远离1查遗传力估计是否正常如果单性状h2本身就没跑出来r_g自然失真查表型定义。人和小鼠的表型在测量尺度上往往不对齐比如“焦虑样行为”在不同物种里用不同测试范式衡量即使遗传相关高也未必反映为高r_g。这些排查都做完之后如果r_g还是接近0那更可能是生物学事实而非技术错误。这个结论同样有价值它说明该性状的遗传基础在进化中发生了改变保守的信号可能只集中在少数区域而不是全基因组尺度。4.4 问题速查表现象最可能原因快速解法映射率低于40%链文件版本不匹配核对参考基因组版本使用正确链文件munge后SNP过少原始数据chr命名不规范统一为chr前缀格式后再映射h2估计为负等位基因方向未对齐检查A1/A2频率方向整体取补截距远大于1样本重叠或人群分层用--intercept-h2固定为1重跑测试r_g结果与预期不符共同SNP集合过小用--keep-snps限制分析集合运行时module错误Python版本不兼容用Python3.6环境跑LDSC4.5 两个容易忽略的实操细节除了上面的问题还有两个细节我在实操里反复踩值得单独拿出来说。第一个是--w-ld-chr参数千万不能漏。LDSC的权重矩阵依赖SNP的LD Score和等位基因频率如果不提供权重文件程序会用默认值结果在遗传力估计上会出现系统性偏差。很多人以为跨物种分析需要自己生成权重文件其实不需要直接用官方提供的那套权重即可因为权重只和人类参考面板有关和GWAS数据来源无关。第二个是--out文件名不要用中文和空格。LDSC内部解析路径时对这类字符处理不太友好报错信息又不容易看懂我有一次排查半天就是文件名里多了一个空格。5. 后续还可以怎么扩展跨物种LDSC分析做完基础版之后有几个方向非常值得延伸。一个是结合S-LDSC做功能注释富集。你把人和小鼠GWAS分别跑一遍分层LDSC然后把各个注释类别的富集倍数放在一起比较比如增强子、启动子、保守非编码区、3UTR等等。只要两组数据都映射到同一套人类基因组注释上这个比较在计算上是完全可行的。实际做出来以后你会发现有些注释类别在人里富集、在小鼠里不富集这正是进化分化的信号论文里很好讲故事。另一个方向是利用跨物种r_g做因果推断辅助。比如你同时有抑郁相关的人类GWAS和小鼠应激行为GWAS算出r_g显著为正之后可以继续用跨物种共定位分析找共享因果位点。LDSC解决的是“有没有共享遗传结构”的问题共定位解决的是“具体共享哪个位点”的问题两者互补。最后还有一个偏方法学的点如果你有多个物种的数据比如人、小鼠、猪、狗都能拿到GWAS可以考虑把LDSC的交叉分析推广成一个两两矩阵。矩阵里每个元素都是r_g看起来就像一棵热图树能直观展示性状遗传基础的种间亲疏关系。这个做法计算量不大但信息量很高。我在实际跑跨物种LDSC的时候最大的体会是这个分析真正难的不是命令怎么写而是你时刻要清楚每一步在算什么。坐标映射改变了SNP的位置标签等位基因对齐改变了效应量的符号LD Score参考决定了回归的自变量任何一个环节的理解偏差都会让结果在数值上“没问题”但生物学上完全误导人。最后再分享一个小技巧跑完整套流程后记得把两组的summary statistics、映射中间文件、LD Score子集文件、输出日志都归档到一个目录里最好连脚本参数也存成README。跨物种分析的链条很长三个月后你回头再看结果如果只有最终r_g那一个数字你会完全想不起来它怎么来的。留下来一套完整记录不管是写论文的方法部分还是应对审稿人的复核意见都会轻松得多。