生信碱移
PCMR分析
PCMR 是一种基于遗传变异的聚类方法，用于检测相关水平多效性并增强因果推断。其核心思想是将遗传工具变量根据多效性模式分类，通过混合高斯模型区分因果效应与多效性效应，避免假阳性结果。
探究复杂性状之间的因果关系并确定疾病的因果风险因素，对于揭示各种疾病的病因学至关重要。孟德尔随机化 (Mendelian Randomization, MR) 是一种利用遗传变异作为工具变量（IVs）推断暴露因素（如血脂）与结局（如心血管疾病）间因果关系的统计方法。MR 本质的目标是去判断暴露是否会影响结局，这与一些队列观察性研究的目标是一致的（探索血脂是否会影响心血管疾病）。其核心假设是工具变量需满足以下条件：
关联性：工具变量与暴露因素强相关；
独立性：工具变量与混杂因素无关（如环境因素）；
排他性：工具变量仅通过暴露因素影响结局（无水平多效性）。
如果看不懂上方文章，强烈推荐小编之前的推文，一文读懂孟德尔随机化，点击此处跳转查看。
比如，广泛使用的逆方差加权（IVW）方法假设所有显著的遗传变异都遵循这些核心原则，从而成为有效的工具变量。工具变量的引入以及限制是为了消除混淆因素介导的因果影响，常规使用的工具变量来自 GWAS 研究中被鉴定到的遗传位点（SNP）。值得一提的点在于，当遗传变异通过暴露以外的路径直接影响结局时，这种现象被称为水平多效性。具体分为两类：
相关水平多效性（Correlated Horizontal Pleiotropy）：IVs 通过共享的混杂因素（如代谢通路）同时影响暴露（血脂水平）和结局（心血管疾病），导致因果效应估计偏倚。
无关水平多效性（Uncorrelated Horizontal Pleiotropy）：IVs 直接作用于结局，与暴露无关。SNP直接与心血管疾病关联，但是我们的研究目的是血脂与结局，SNP 只是我们需要的工具。
▲ 水平多效性示意图：图中标出了两种水平多效性类型。
传统MR方法（如IVW、Egger回归）假设水平多效性效应可忽略或服从特定分布，但在相关水平多效性变体占比较高时（如30-40%），这些方法易产生假阳性结果。还是前面的例子，我们通过 MR 分析推断高密度脂蛋白与冠心病存在因果关系，但是实际是两者共享类似的代谢通路紊乱导致的多效性干扰，并不存在因果关系。来自中山大学遗传与生物信息学系的李淼新课题组开发了一种升级的孟德尔随机化方法 PCMR，于2025年3月21日发表于 Nature Communications [IF: 14.7]，基于遗传变异的聚类来检测相关水平多效性并增强因果推断。其核心思想是将工具变量根据多效性模式分类，通过混合高斯模型区分因果效应与多效性效应，避免假阳性结果。作者对 48 种暴露-常见疾病组成的数据集进行了分析，它们在 7 对暴露-疾病中检测到了水平相关多效性，并发现使用 PCMR 避免了检测到假阳性因果联系。
▲ DOI: 10.1038/s41467-025-57912-5。
PCMR 通过三个流程进行MR分析：①首先，将遗传工具变量按多效性模式进行聚类；②随后，异质性检验判断是否存在共享混淆因子导致的多效性；③最后，在存在多效性时，PCMR增强因果效应的稳健性。本文简要介绍 PCMR 的使用，感兴趣的同学可以通过下方链接进一步了解：
https://github.com/856tangbin/PCMR
0. 安装R包
通过以下代码安装PCMR包：
devtools::install_github("856tangbin/PCMR")
library(PCMR)
1. 执行PCMR因果推断
作者使用精神分裂症（SCZ，暴露）到重度抑郁症（MDD，结局）的进行了PCMR的应用演示。
① PCMR 需要两个输入文件：
IVs文件：包含与暴露或结局显著关联的IVs信息
随机样本文件：初始化参数估计的随机SNP集合
上述随机样本变体是通过 R 包 cause 中的函数获得的，参见github仓的./demo/IVs_filter.R代码。这里，我们可以看看示例输入数据的样子：
# 加载示例数据（内置数据集）
data(IVs_and_InitEst, package = "PCMR")
# 查看数据结构
head(IVs_and_InitEst$X_clump) # IVs数据
head(IVs_and_InitEst$X_clump1) # 随机样本数据
② 参数初始化。使用随机样本估计初始值（sigma2为不相关多效性方差，rho为样本重叠相关性）：
set.seed(0)
init <- PCMR_initEst(
beta_hat_1 = X_clump1$beta_hat_1, # 暴露效应量（随机样本）
seb1 = X_clump1$seb1, # 暴露标准误
beta_hat_2 = X_clump1$beta_hat_2, # 结局效应量
seb2 = X_clump1$seb2 # 结局标准误
)
# 查看初始值
print(init)
# $sigma2 = [0.0, 1.99e-05]（不相关多效性方差）
# $rho = 0.0197（样本重叠相关性）
③ 执行多效性聚类分析，将IVs分为2类（默认模型为随机效应模型）
result_random <- PCMR(
beta_hat_1 = X_clump$beta_hat_1, # IVs的暴露效应量
seb1 = X_clump$seb1, # IVs的暴露标准误
beta_hat_2 = X_clump$beta_hat_2, # IVs的结局效应量
seb2 = X_clump$seb2, # IVs的结局标准误
num_gamma = 2, # 聚类类别数（默认就是2类）
model = "1", # 模型选择（1=随机效应，2=固定效应）
isIntact = TRUE,
rho = init$rho,
sigma2 = init$sigma2
)
# 查看聚类结果
print(result_random$gamma) # 各类效应值（θ = β + α）
#33.33333% 66.66667%
#0.01703264 0.18722993
print(result_random$pi_gamma) # 各类别占比
#0.498 0.502
PCMR_plot(result_random)
如上所示，聚类出的两类工具变量效应值分别为0.017和0.187，占比接近1:1。gamma 的分类估计值是相关水平多效性效应和因果效应的总和，称为相关水平多效性效应； gSigma2 和 pi_gamma 分别测量随机相关水平多效性效应和不同类别的比例。默认模型是 PCMR 的随机效应模型（model=1），而 PCMR 的固定效应模型（model=2）将 gSigma2 设为零。
④ 使用异质性检验，判断是否存在显著相关水平多效性：
# 使用估计因子c衡量效应异质性
result_random <- PCMR_cEst(
result = result_random,
ref_beta_outcome = X_clump1$beta_hat_2, # 随机样本的结局效应量
ref_se_outcome = X_clump1$seb2, # 随机样本的结局标准误
cores = 10
)
# Bootstrap估计差异显著性（耗时约30分钟）
result_random <- PCMR_testCausal_bootstrap(result_random, cores = 10)
result_random <- PCMR_testCorPlei(result_random)
print(result_random$CHVP_test) # 输出P值
#[1] 1.13e-07
可以看到现在是存在水平多效性的。
⑤ 存在水平多效性时（作者推荐是P < 0.20），可以通过以下函数进行因果效应评估：
result_random <- PCMR_testCausal(result_random)
print(result_random$Pvalue)
#[1] 0.329 # 因果效应不显著（P > 0.05）
print(result_random$effect)
#[1] 0.1872299 # 最大类别的效应估计
print(result_random$discernable_prob)
#[1] 0.791 # 最大类别为真实因果群组的概率为79.1%
综合上述示例分析的三个流程，可以获得以下结果/结论：
① 多效性聚类：PCMR 将工具变量聚类为两个不同的类别，其相关的水平效应分别为 0.017 和 0.187；
② 异质性检验分析：相关的水平效应
result_random$CHVP_test差异为1.13e-07非常显著，意味着存在水平多效性；③ 因果分析：在考虑水平效应的情况下，PCMR 的因果性评估结果表明，SCZ 对 MDD 的关系不显著（P=0.329）。由最大 IV 类支持的估计因果效应为 0.187(即上方第三步的效应，包含因果效应β和多效性α值)，但可辨别的概率仅为 0.791。
2. 整合生物信息增强因果推断
PCMR 鉴定的不同类工具变量可以帮助整合通路注释信息用来进行机制分析。
① 首先，可以使用下面的代码提取不同类别的工具变量：
prb_thrd <- 0.5 # 概率阈值
probability_1 <- rowSums(result_random$W[, 1, ]) # 属于第1类的概率
# 提取两类IVs的rsID
gray_rsid <- X_clump$rsid[probability_1 > prb_thrd] # 第1类（HVP效应较小）
blue_rsid <- X_clump$rsid[probability_1 < 1 - prb_thrd] # 第2类（HVP效应较大）
# 输出结果（以逗号分隔）
cat("Gray类IVs（HVP效应=0.017）:\n", paste(gray_rsid, collapse = ", "))
##[1] "rs6715326 rs12422904 rs12891289 rs113213 rs422118 rs1123263861 rs1231032367
cat("\nBlue类IVs（HVP效应=0.187）:\n", paste(blue_rsid, collapse = ", "))
#[1] "rs6715366 rs12432904 rs12892189 rs11687313 rs3739118 rs11263861 rs12310367 rs10777956 rs4702 rs67439964 rs11136325 rs4129585 rs3824451 rs7927176 rs12293670 rs12652777 rs149165 rs7830315 rs11941714 rs13145415 rs117145318 rs2333321 rs61405217 rs12498839 rs6943762 rs13233308 rs16851048 rs4636654 rs11027839 rs2456020 rs198806 rs12199613 rs1796518 rs6909479 rs4712938 rs2531806 rs13195636 rs35531336 rs1265099 rs9272446 rs9461856 rs116416291 rs9461916 rs570263 rs9274623 rs209474 rs1144708 rs2967 rs11693094 rs12129573 rs6762456 rs12489270 rs5995756 rs9623320 rs699318 rs4812325 rs60135207 rs13016542 rs7900775 rs11191580 rs79780963 rs2815731 rs9910403 rs8055219 rs187557 rs10860960 rs1604060 rs2241033 rs56393513 rs4936277 rs4987094 rs2514218 rs1881046 rs79210963 rs634940 rs9304548 rs2710323 rs13080668 rs7634476 rs11638554 rs708228 rs35351411 rs3814883 rs2190864 rs7191183 rs72692857 rs13107325 rs7403630 rs11807834 rs72769124 rs145071536 rs9428966 rs4658559 rs2489213 rs7099380 rs6974218 rs74480281 rs11682175 rs12969453 rs74914300 rs72934602 rs4632195 rs9636107 rs505061 rs9318627 rs13407231 rs6701322 rs56205728 rs2255663 rs117472063 rs11664298 rs10894308 rs4936215 rs7001340 rs61786047 rs6520064 rs2238057 rs740417 rs12712510 rs6125656 rs926288"
② 随后，将上面输出的gray_rsid和blue_rsid上传至SNPense，获取对应的基因名称。然后使用g:Profiler在线工具对这些基因做富集分析，结果示例如下：
如上所示，Blue类基因富集于“化学突触传递”，提示可能与神经信号通路相关（潜在多效性的来源）。
这个工具可以试试
欢迎各位佬哥佬姐关注