两个性状在同一区域都出现关联信号,不代表它们由同一个 causal variant 驱动。它们可能共享变异,也可能只是各自的 causal variants 处于 linkage disequilibrium(LD)中。
coloc 使用区域内的 summary statistics 比较这些解释。常见场景包括 GWAS 与 eQTL、pQTL 或 mQTL 的比较,也可以比较两个疾病或连续性状。
library(coloc)
五个假设
经典 coloc.abf() 为一个区域计算五种假设的 posterior probability:
| 假设 | 含义 |
|---|---|
| H0 | 两个性状在区域内都没有关联 |
| H1 | 只有性状 1 有关联 |
| H2 | 只有性状 2 有关联 |
| H3 | 两个性状都有信号,但由不同 causal variants 驱动 |
| H4 | 两个性状的信号由同一个 causal variant 驱动 |
分析重点通常是比较 PP.H3 与 PP.H4,而不是只问 PP.H4 是否超过某个固定阈值。若区域对其中一个性状缺乏足够关联证据,H0、H1 或 H2 可能占主导,此时不能把较低的 H4 简单解释为“不共定位”。
输入数据结构
每个性状分别整理为一个 named list。优先使用效应值与方差:
| 字段 | 含义 | 说明 |
|---|---|---|
beta |
SNP effect estimate | 与 varbeta 配套 |
varbeta |
effect estimate variance | 通常为 se^2 |
snp |
variant identifier | 两个数据集应匹配同一批变异 |
position |
genomic position | 检查区域和绘图时使用 |
type |
"quant" 或 "cc" |
定量或病例对照性状 |
sdY |
outcome standard deviation | 定量性状需要,部分情况下可估计 |
s |
case fraction | 病例对照性状需要 |
N |
sample size | 缺少部分参数或使用 p 值输入时需要 |
MAF |
minor allele frequency | 参数估计及 p 值输入时需要 |
有 beta 和标准误时:
dataset1 <- list(
beta = gwas$beta,
varbeta = gwas$se^2,
snp = gwas$snp,
position = gwas$position,
type = "cc",
s = gwas$n_case / gwas$n_total,
N = gwas$n_total,
MAF = gwas$maf
)
dataset2 <- list(
beta = eqtl$beta,
varbeta = eqtl$se^2,
snp = eqtl$snp,
position = eqtl$position,
type = "quant",
N = eqtl$n,
MAF = eqtl$maf,
sdY = eqtl$sd_y
)
缺少 beta 与 variance 时,coloc 也可以使用 pvalues、N、MAF 等信息进行近似,但完整效应估计通常更合适。不要把标准误直接传给 varbeta:
varbeta <- se^2
先对齐区域和变异
两个数据集必须使用相同 genome build,并在同一分析区域中取交集:
common_snps <- intersect(gwas$snp, eqtl$snp)
gwas_region <- gwas[gwas$snp %in% common_snps, ]
eqtl_region <- eqtl[eqtl$snp %in% common_snps, ]
随后按同一 SNP 顺序排列,并检查:
- rsID 或
chr:position:alleles是否指向同一变异; - effect allele 与 other allele 是否已经 harmonize;
- 是否混用了 genome builds;
- palindromic variants 的方向能否可靠判断;
- MAF、样本量和 case fraction 是否来自相应数据集;
- 区域是否完整覆盖两边的关联信号,而不是只保留 lead SNP。
coloc 使用区域内完整的关联模式。只拿各自显著 SNP 的交集,会破坏这个模式并产生误导。
检查数据集
内置测试数据可以展示最小结构:
data(coloc_test_data)
D1 <- coloc_test_data$D1
D2 <- coloc_test_data$D2
minimum_data <- D1[c(
"beta", "varbeta", "snp", "position", "type", "sdY"
)]
check_dataset(minimum_data)
plot_dataset(minimum_data)
对真实数据分别运行检查:
check_dataset(dataset1)
check_dataset(dataset2)
检查通过只说明结构满足函数要求,不证明 allele alignment、研究人群或生物学区域选择正确。
单信号分析:coloc.abf()
经典分析假设每个性状在该区域最多有一个 causal signal:
result <- coloc.abf(
dataset1 = dataset1,
dataset2 = dataset2
)
result$summary
查看对 H4 贡献较高的 variants:
head(
result$results[order(result$results$SNP.PP.H4, decreasing = TRUE), ],
10
)
SNP.PP.H4 是在 H4 成立条件下各 SNP 成为共享 causal variant 的相对支持,不应脱离区域级 PP.H4 单独报告。
Prior sensitivity
H4 的 posterior probability 会受 priors 影响,尤其是 p12,即一个 variant 同时关联两个性状的先验概率。
result <- coloc.abf(
dataset1 = dataset1,
dataset2 = dataset2,
p12 = 1e-6
)
sensitivity(result, rule = "H4 > 0.5")
规则应在分析计划中说明,并结合 H3、区域信号强度和研究目的解释。若结论只在很窄的 prior 范围内成立,应将这种不稳定性报告出来。
多信号区域:coloc.susie()
同一区域可能存在多个独立信号。此时单 causal variant 假设不合适,可以提供 LD matrix,先用 SuSiE 分解信号,再进行共定位。
check_dataset(dataset1, req = "LD")
check_dataset(dataset2, req = "LD")
susie1 <- runsusie(dataset1)
susie2 <- runsusie(dataset2)
result_susie <- coloc.susie(susie1, susie2)
result_susie$summary
LD matrix 必须:
- 为 SNP × SNP 方阵;
- 行名、列名和 dataset 中的
snp完全对应; - SNP 顺序一致;
- 来自与 summary statistics 祖源尽量匹配的参考人群。
identical(rownames(dataset1$LD), dataset1$snp)
identical(colnames(dataset1$LD), dataset1$snp)
参考 LD 不匹配会改变信号分解。coloc.susie() 能处理多个信号,但不能修复错误的 LD 或 allele alignment。
结果边界
共定位支持两个关联信号与共享 causal variant 相容,但它本身不证明:
- 一个性状导致另一个性状;
- 对应基因就是疾病机制中的 effector gene;
- 共享变异只通过所测分子性状发挥作用;
- 不同组织、细胞类型或祖源中的信号相同。
报告时应写清楚区域、build、样本量、性状类型、variant 数量、priors、LD 来源、H0–H4 posterior probabilities,以及是否采用多信号模型。