第一次跑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分析,基本上可以归成三大类,每类的目的和数据要求差别很大。
第一类叫保守性富集分析。思路是:先定义一组跨物种保守区域(比如人类和小鼠的同源编码区、保守调控元件、受纯化选择的区间),然后用分层LDSC(stratified 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环境跑LDSC |
4.5 两个容易忽略的实操细节
除了上面的问题,还有两个细节我在实操里反复踩,值得单独拿出来说。
第一个是--w-ld-chr参数千万不能漏。LDSC的权重矩阵依赖SNP的LD Score和等位基因频率,如果不提供权重文件,程序会用默认值,结果在遗传力估计上会出现系统性偏差。很多人以为跨物种分析需要自己生成权重文件,其实不需要:直接用官方提供的那套权重即可,因为权重只和人类参考面板有关,和GWAS数据来源无关。
第二个是--out文件名不要用中文和空格。LDSC内部解析路径时对这类字符处理不太友好,报错信息又不容易看懂,我有一次排查半天,就是文件名里多了一个空格。
5. 后续还可以怎么扩展
跨物种LDSC分析做完基础版之后,有几个方向非常值得延伸。
一个是结合S-LDSC做功能注释富集。你把人和小鼠GWAS分别跑一遍分层LDSC,然后把各个注释类别的富集倍数放在一起比较,比如增强子、启动子、保守非编码区、3'UTR等等。只要两组数据都映射到同一套人类基因组注释上,这个比较在计算上是完全可行的。实际做出来以后你会发现,有些注释类别在人里富集、在小鼠里不富集,这正是进化分化的信号,论文里很好讲故事。
另一个方向是利用跨物种r_g做因果推断辅助。比如你同时有抑郁相关的人类GWAS和小鼠应激行为GWAS,算出r_g显著为正之后,可以继续用跨物种共定位分析找共享因果位点。LDSC解决的是“有没有共享遗传结构”的问题,共定位解决的是“具体共享哪个位点”的问题,两者互补。
最后还有一个偏方法学的点:如果你有多个物种的数据,比如人、小鼠、猪、狗都能拿到GWAS,可以考虑把LDSC的交叉分析推广成一个两两矩阵。矩阵里每个元素都是r_g,看起来就像一棵热图树,能直观展示性状遗传基础的种间亲疏关系。这个做法计算量不大,但信息量很高。
我在实际跑跨物种LDSC的时候,最大的体会是:这个分析真正难的不是命令怎么写,而是你时刻要清楚每一步在算什么。坐标映射改变了SNP的位置标签,等位基因对齐改变了效应量的符号,LD Score参考决定了回归的自变量,任何一个环节的理解偏差,都会让结果在数值上“没问题”但生物学上完全误导人。
最后再分享一个小技巧:跑完整套流程后,记得把两组的summary statistics、映射中间文件、LD Score子集文件、输出日志都归档到一个目录里,最好连脚本参数也存成README。跨物种分析的链条很长,三个月后你回头再看结果,如果只有最终r_g那一个数字,你会完全想不起来它怎么来的。留下来一套完整记录,不管是写论文的方法部分还是应对审稿人的复核意见,都会轻松得多。