跨物种LDSC遗传相关性计算：原理、流程与实战避坑指南

“LDSC跨物种计算”这个活儿，我从第一次听到就觉得很带感。你手里同时握着人和小鼠的两套全基因组关联分析结果，想回答“这两个物种在某个复杂性状上的遗传架构到底像不像”，这就不是简单跑个相关性就能糊弄过去的。LDSC，全称Linkage Disequilibrium Score Regression，中文一般叫连锁不平衡分数回归，它最妙的地方在于只靠GWAS汇总统计量就能估计遗传力、遗传相关性这些参数，不需要个体级基因型数据。而一旦把两个物种的GWAS结果放进同一个LDSC框架里去算，就变成了跨物种遗传相关性分析，这在进化遗传学、动物育种、疾病模型验证里都是高频需求。

这里先给刚接触的朋友交个底：LDSC跨物种计算要做的事，就是把物种A和物种B的GWAS summary statistics，经过同源位点映射、参考面板统一、质量过滤之后，喂给LDSC回归，输出一个遗传相关性（genetic correlation，rg）估计值。这个值告诉你两套遗传信号的共享程度，比表型相关性更接近“基因层面的同源性”。它适合谁用？我总结下来至少三类人：一是做模式动物研究的，想验证小鼠模型某个性状的遗传基础能不能代表人；二是做数量遗传学或育种的，想比较不同物种间同一经济性状的遗传结构是否保守；三是搞方法学或组学整合的，需要在大规模汇总数据里快速估算跨物种共享遗传架构。无论你是生信新手还是有经验的从业者，这套流程都值得认真跑一遍。

1. LDSC跨物种计算在做什么

1.1 它到底算出了一个什么数字

先把LDSC的原理用最直白的方式讲清楚。GWAS跑完之后，每个SNP会得到一个卡方统计量（chi-square statistic），它反映了这个位点与性状的关联强度。LDSC的核心观察是：一个SNP的期望卡方值，与它周围所有SNP的连锁不平衡状态有关。这个“周围连锁不平衡的量”就是LD score，它是某个SNP与其周围SNP的r²之和。回归模型可以写成这样：

E[χ²] = 1 + N·h²·l_j / M + 截距项

这里N是GWAS样本量，h²是SNP遗传力，l_j是第j个SNP的LD score，M是有效SNP数量。用LD score做回归，斜率里藏着遗传力信息，截距则反映了群体分层、样本重叠等混杂因素。

跨物种场景下，我们要把两个物种各自的GWAS信号放在一起比较。方法上是把两个GWAS在同一个SNP位点上的效应量或者z-score拿出来，计算它们的“跨物种z-score乘积”，然后用每个SNP的LD score做回归，输出一个rg值。更直观地说：如果两个物种在几乎相同的因果位点上有相似方向的效应，rg就会接近1；如果共享的因果位点很少或者效应方向不一致，rg就会显著低于1甚至接近0。所以rg估计的不是表型相关，而是“因果变异在基因组层面上的共享程度”。

1.2 为什么跨物种场景下要格外谨慎

同物种的LDSC分析，参考面板、SNP集合、等位基因频率都相对好处理，跨物种一掺和进来，麻烦就成倍增长。我总结出几个必须死盯的关键点。

第一，参考面板必须统一。LDSC的LD score文件是基于特定人群参考面板算出来的，比如欧洲人群（EUR）、东亚人群（EAS）等。跨物种计算时，无论你是拿人鼠比较，还是人猪比较，都必须给两个物种使用同一个参考面板下的LD score文件。你不可能拿人类EUR面板的LD score去匹配小鼠基因组上的SNP，因为LD结构完全不同。

第二，SNP空间要能对齐。人基因组有约千万级SNP，小鼠有大几百万，两个基因组的坐标、等位基因编码方式都不一样。跨物种计算要求我们找到一组“同源位点”——通常是把非人物种的GWAS位点映射到人类基因组坐标上，或者通过直系同源关系确定SNP对应关系。这个环节做不好，后面交集数量骤减，rg估计就会出现巨大的标准误。

第三，效应方向的一致性。LDSC回归里的z-score是有方向的，如果A1/A2等位基因在两个物种的GWAS里定义相反，效应方向的判断就会直接颠倒。跨物种数据来自不同实验室、不同芯片平台，等位基因链方向错误几乎是必然会遇到的事。

第四，样本量与统计功效差异。人GWAS动不动几十万样本，小鼠GWAS常见样本量往往只有几千甚至几百。样本量小的物种，z-score噪声大，LDSC回归斜率会向0收缩，rg估计也就容易偏低。这不是生物信号弱，而是统计噪声在作怪。

把同物种分析和跨物种分析放一起对比，差异就特别明显。

所以，跨物种LDSC的结果不是一个“普适公式”输出，而是需要非常克制地解读。接下来我把跑通这个流程的完整数据链路拆开讲。

2. 完整的跨物种数据链路

2.1 两个物种的GWAS汇总数据从哪来

启动一个跨物种LDSC分析，第一步自然是备齐两个物种的GWAS summary statistics。人的数据好找，GWAS Catalog、Open Targets、各类联盟公开数据都行。非人物种就需要费点心思，小鼠可以找IMPC、Jackson Laboratory的公开数据，猪牛羊鸡可以参考Animal QTLdb和各类育种联盟。

不管数据来自哪，最终喂给LDSC的汇总统计文件必须包含几列关键信息：SNP标识（rsID或chr:pos）、效应等位基因A1、非效应等位基因A2、效应量（beta或OR）、标准误或z-score、P值，以及等位基因频率。如果有样本量N列，就直接用--N-col指定；如果没有，就得从文章里查出来写在命令行。

这里有个经常被忽略的点：文件里的染色体坐标版本。人类数据有hg19、hg38之分，小鼠有mm9、mm10等版本。如果你拿到的两个文件版本不一致，映射时就会损失大量位点。我的习惯是第一步先用awk或者R脚本检查一下染色体编号和位置分布，确认版本之后再往下走。

2.2 物种间SNP映射怎么处理

这是跨物种LDSC计算里最像“手艺活”的环节。我试过几种方案，最终最省心的是“统一到人类基因组坐标”加“直系同源映射”的组合拳。

方案A：如果非人物种的GWAS数据里有rsID，可以尝试直接和人类GWAS的rsID取交集。问题是很多非人物种芯片的rsID并不规范，甚至完全是自定义ID，直接取交集可能损失大半位点，必须谨慎。

方案B：使用Ensembl BioMart获取直系同源映射表。这个表会告诉你小鼠基因和人类基因之间的直系同源关系，但它的粒度是基因不是SNP。你需要先把非人GWAS的SNP注释到基因或转录本上，再通过同源关系带回人类坐标。这个过程比较繁琐，但对于做跨物种性状比较的研究来说是最严谨的路径。

方案C：使用UCSC的liftOver工具，把非人物种的坐标直接转成人类基因组坐标。这种方法依赖两个物种基因组之间的链比对，适合在保守基因组区域进行粗映射。实操上我一般会同时保留坐标转换结果和rsID匹配结果，合并去重后只保留那些“双策略一致”的位点，这样能显著降低假阳性对应。

做SNP映射时，我强烈建议你保留一张对照表，至少要包含：原始物种的chr、pos、rsID、A1、A2，以及映射后的人类chr、pos、rsID。后面出问题时，这张表就是排查的底图。

2.3 参考面板与LD Score文件统一

LDSC官方发布了一套基础数据文件，包括HapMap3 SNP列表（w_hm3.snplist）和分染色体的LD score文件（eur_w_ld_chr/）。这套文件里的LD score是基于千人基因组欧洲人群计算的。不同人群的LD结构有差异，但官方分析里使用EUR面板是默认做法，因为HapMap3这组SNP在多个数据库里都很好匹配。

跨物种场景下，LD score文件必须统一。人的GWAS用EUR面板，小鼠的GWAS也必须用同一个EUR面板下的LD score文件。这里不是“入乡随俗”，而是要让两个物种的回归在完全相同的SNP集和LD权重体系下比较，否则你会把LD结构的差异错误识别成遗传架构差异。有一个细节要注意：如果人的GWAS是东亚人群，你依然要选择EUR面板作为跨物种基准，而不是选EAS面板混搭，因为参考面板的不一致会直接污染rg估计。

我自己跑跨物种分析时，会把w_hm3.snplist先解压，然后检查两个物种的GWAS经过映射后，有多少位点能落到HapMap3集合里。这个交集数量直接决定了结果能不能信。一般来说，交集少于5万个SNP时，rg的置信区间会宽到让你怀疑人生。

2.4 munge_sumstats.py预处理实操

预处理是整个流程里最机械但最容易埋雷的步骤。LDSC官方工具munge_sumstats.py会把任意格式的GWAS汇总统计转成标准格式，并做质量过滤。

我常用的命令长这样：

bash复制python munge_sumstats.py \
  --sumstats human_gwas.txt \
  --out munged_human \
  --merge-alleles w_hm3.snplist \
  --a1 A1 --a2 A2 --snp SNP \
  --info INFO --info-min 0.9 \
  --maf MAF --maf-min 0.01 \
  --N-col N

解释一下几个关键参数：--merge-alleles会把输入的SNP列表和HapMap3列表合并，只保留能对上的位点，同时检查等位基因是否与参考面板一致；--info-min和--maf-min是常规质量过滤；--N-col直接指定样本量列。

小鼠等非人物种的汇总文件做munge时也一样，只是SNP列可能不是rsID而是chr:pos。那就需要先在映射环节已经把转换成人类rsID的列准备好，或者直接按chr:pos与参考面板匹配。实际跑的时候，munge日志会输出类似“SNPs with chi-square > 99999 were removed”“SNPs in sumstats but not in reference panel”这类提示，每一行都要看，别直接关掉终端。

3. 跑通两物种rg计算

3.1 命令行实操

等到两个物种各自生成.sumstats.gz文件之后，就可以进入核心计算了。LDSC的--rg参数接收至少两个gzip文件，格式是逗号分隔，第一个文件是“base”，第二个是“target”。命令长这样：

bash复制python ldsc.py \
  --rg munged_human.sumstats.gz,munged_mouse.sumstats.gz \
  --ref-ld-chr eur_w_ld_chr/ \
  --w-ld-chr eur_w_ld_chr/ \
  --out human_mouse_rg \
  --n-blocks 200

--ref-ld-chr指定的是用来算LD score的参考面板，--w-ld-chr指定的是回归时加权用的LD score。官方建议这两个都用同一套，跨物种分析更是绝对不能混用。

跑完之后，你会在human_mouse_rg.log里看到完整回归过程，在human_mouse_rg.rg里看到核心结果。.rg文件是下面这种格式：

code复制p1 p2 rg se z p h2_obs h2_obs_se
munged_human.sumstats.gz munged_mouse.sumstats.gz 0.6345 0.0812 7.82 4.5e-15 0.2512 0.0133

这种结果是我比较喜欢看到的样子。p1是第一个物种，p2是第二个物种，rg就是我们要的遗传相关性，se是标准误，z是z分数，p是显著性，后面两列是第一个物种的观测遗传力估计值及其标准误。

3.2 关键参数背后的逻辑

关于--n-blocks，很多教程没重点讲。LDSC在回归时为了获得稳健的标准误，会把基因组划分成一些区块，通过删除每一块重新估计来校正标准误。默认是200块，但如果你跨物种分析里有效SNP数比较少，我建议适当减少到100或者50，不然有些区块里的SNP太少，jackknife结果会变得非常不稳定。

另一个值得留意的参数是--chisq-max。某些GWAS文件里会有个别SNP的卡方值异常高，这些位点通常是拷贝数变异区域或者比对错误的伪信号，会把回归斜率拽偏。实操中我会按照数据分布设一个上限，比如卡方最大值控制在80左右。跨物种场景下我更建议设置这个阈值，因为映射带来的错位可能制造假的极端值。

还有一个容易忽略的点是--intercept-h2。如果两个物种之间可能有样本重叠（比如同一批个体测了两个表型），或者存在群体分层，回归截距就会偏高。默认情况下LDSC会同时估计截距并输出，你不需要手动指定。只有当截距远低于1时，比如0.9以下，才需要考虑是不是过滤太狠或者映射错误杀掉了真实信号。

3.3 判读结果与特殊情况

拿到rg值之后，别急着写结论。跨物种rg的解读逻辑与同物种分析有些不同。

rg接近1，且置信区间不跨0，说明两个物种在共享的因果位点上表现出高度一致的效应方向与幅度，这是“遗传架构保守”的强证据。rg在0.4到0.7之间，说明存在一定程度的共享成分，但两个物种也各自演化出了独特的遗传结构，这种情形在人类与小鼠的某些复杂性状上很常见。rg接近0，意味着共享因果变异很少或者效应方向高度不一致。这里要特别提醒一点：rg不显著不等于“没有共享遗传架构”，也可能是某个物种的GWAS功效不足。

我在实际项目中见过一个典型场景：人骨密度GWAS和小鼠骨密度GWAS的rg跑出来只有0.2左右，但置信区间跨0。后来排查发现小鼠那套GWAS样本量只有2000多，有效SNP数也少得可怜，继续把鼠的汇总数据做个更严格的MAF过滤后，rg升到了0.58。这个例子说明，跨物种rg对比的是两个样本各自的估计精度，样本功效差异掩盖真实共享信号的情况非常常见。

4. 常见问题与排查方案实录

4.1 输出NaN或交集SNP太少

我第一次跑跨物种LDSC时，出来的.rg文件里rg直接是NaN，log里一大片“SNPs with chi-square > 99999 were removed”以及“Warning: fewer than 50 SNPs remain”。多数情况下是交集SNP数量太少导致的。排查时我一般会拆成三步走。

第一步，检查munge之后的.sumstats.gz有多少行。如果只有几千行，说明映射或MAF过滤太严，需要回退。
第二步，检查参考面板是否匹配。--merge-alleles会把输入数据卡到HapMap3列表，如果你用的GWAS是低密度芯片数据，很多位点本来就不在HapMap3集合里，交集会惨不忍睹。这时候需要考虑换用更宽泛的SNP集合，或者调整映射策略。
第三步，检查两个物种的坐标版本。我曾经因为人类数据用的hg19、小鼠数据用的是GRCm38，而构建映射表时没注意版本差异，导致几千个位点全部偏移，后来统一版本后交集数量直接翻倍。

4.2 strand flip与等位基因错位

这是跨物种分析里最隐蔽的坑。人类和小鼠的芯片不同，等位基因编码方式也不同，有些位点在两个文件里A1和A2的链方向是反的。munge阶段虽然会在同一碱基上检查A1/A2是否匹配，但遇到互补链时它可能无法自动识别。

出现这个问题的典型表现是：munge日志里大量出现“Attempting to fix allele mismatches”但fix率很低，最终rg标准误巨大，或者结果明显不合理。

我的解决办法是，在映射阶段就对每个位点做等位基因标准化。利用dbSNP的链信息，把两个物种的等位基因都转换到正链上，然后再做匹配。如果条件不允许，至少也要在munge命令里使用--merge-alleles并检查输出文件里的A1/A2是否与参考面板一致。实际操作中，我会用PLINK的--flip思路写一个小脚本，先统计A1/A2与参考面板的匹配率，匹配率低于95%的位点自动剔除或翻转，这比让munge自己判断省心得多。

4.3 参考面板与群体背景不匹配

跨物种计算中，把所有物种的GWAS都映射到人类参考面板上是正确的做法，但这里有一个微妙问题：如果人类那侧的GWAS是东亚人群，而你统一用了EUR面板，人的LD结构匹配就会变差。LD score的偏差会在回归里引入噪声。

实际操作中我会做两个敏感性检验来评估这个影响。第一，分别在EUR和EAS面板下跑同物种LDSC，看看h²估计值差多大。第二，在跨物种分析里分别用两个面板跑一次，如果rg结果差异较大，说明LD结构不匹配带来的影响不可忽视。跨物种流程本质上是一个标准化过程，统一基准比追求单侧最优更重要，所以EUR面板往往成为默认基准。但也别把“默认”当真理，每次都要检查log里的回归斜率稳定性。

还有一种不匹配来自等位基因频率差异。人类和小鼠在某个同源位点上的MAF可能一个是0.35另一个是0.03。MAF极低的位点，LD score估计误差大，z-score噪声高，留它在分析里只会制造干扰。我会在实际流程中增加一个“交叉频率过滤”的步骤：只保留两个物种MAF都大于0.05的位点。这个阈值可以根据有效SNP数动态调整，如果交集太少就放宽到0.03。

4.4 多组合并行计算的小窍门

如果你不只是比较人和小鼠，而是要同时比较多组物种，比如人-鼠、人-猪、人-犬，那就不适合一组一组手动敲命令。我在项目里会把munge和LDSC封装成两个bash脚本，用物种ID做循环。

bash复制for sp in mouse pig dog;
do
  python munge_sumstats.py \
    --sumstats ${sp}_gwas.txt \
    --out munged_${sp} \
    --merge-alleles w_hm3.snplist \
    --a1 A1 --a2 A2 --snp SNP \
    --info INFO --info-min 0.9 \
    --N-col N
done

python ldsc.py \
  --rg munged_human.sumstats.gz,munged_mouse.sumstats.gz,munged_pig.sumstats.gz,munged_dog.sumstats.gz \
  --ref-ld-chr eur_w_ld_chr/ \
  --w-ld-chr eur_w_ld_chr/ \
  --out multi_species_rg \
  --n-blocks 200

LDSC的--rg支持多于两个文件，它会给出所有两两组合的相关性矩阵。这个功能很实用，但输出文件会变得比较大。建议每次跑完都记下总SNP数和交集情况，不然同一套数据跑三个物种，结果文件混在一起时，分析哪个组合有效哪个组合是噪声就成了新的麻烦。

4.5 我踩过的几个具体的坑

最后分享几个让我印象深刻的现场事故，希望能帮大家少走弯路。

第一个坑：小鼠GWAS用了错误的染色体命名。小鼠有19条常染色体和性染色体，官方文件里有时是“1”，有时是“chr1”，如果你的脚本里没有统一前缀，映射表会丢一半位点。这个问题的排查其实很简单，只要在munge日志里看一眼染色体分布就露馅了。

第二个坑：混用了hg19和hg38。我有一次跑人犬比较时收到一份人的GWAS坐标是hg19，另一份参考表是hg38，直接按位置取交集后大量位点错位。后来我写了统一坐标版本的脚本，先把所有数据liftOver到hg38，再重新映射，才把问题解决。跨物种分析里坐标版本是事故高发区，拿到数据第一件事必须是确认版本。

第三个坑：munge阶段对小鼠数据用了过高的--info-min。小鼠很多公开数据集的imputation质量评分体系与人类不同，直接按人类0.9的阈值过滤后，小鼠GWAS的有效SNP只剩不到四分之一，rg估计直接崩了。后来我针对小鼠数据把--info-min下调到0.6，交集数恢复到了可接受范围。不同物种的数据质量指标不能一概而论，阈值设定要单独看分布。

第四个坑：把.sumstats.gz文件交给LDSC之后再重复映射。LDSC只认位置和等位基因格式，它不会帮你重新做物种映射。如果前一步映射有误，后一步再努力也只是把错误保持一致。所以我会在映射阶段反复校验，先把物种内GWAS和参考面板的等位基因匹配率拉高到95%以上，再进LDSC。

写在实际操作之后

跨物种LDSC跑得多了，我最大的体会是：这个工具的输出并不复杂，复杂的是让两个物种的数据在同一个标准化流程下真正“可比”。参考面板、坐标版本、等位基因方向、质量过滤阈值，任何一个环节偷懒，都会让rg变成一堆不知所云的数字。我个人习惯在做跨物种分析之前，先对每个物种单独跑一遍同物种LDSC的h2估计，确认单侧数据质量没有硬伤，再进入跨物种比较。这样就算rg结果不理想，你也能清楚问题出在数据源而不是方法。

另外建议所有做跨物种分析的团队，把munge后的.sumstats.gz文件和映射表作为标准产出物保存下来。这些中间文件不仅方便复现，还能在之后加入新物种时直接复用。我自己就是靠这套保存下来的中间文件，把原本要跑两天的数据链路缩短到半天以内。希望这篇文章能帮你把LDSC跨物种计算的每一步都踩稳。