LDSC遗传相关性分析教程

🧬 新手导读：LDSC和MR是什么关系？

LDSC不是MR替代品，而是前置评估工具：先看两性状是否存在共享遗传结构，再决定MR优先级。

背景介绍

LDSC可估计遗传相关性（rg）和截距，帮助判断性状间是否可能存在共同遗传基础，以及结果是否受群体结构偏差影响。

备注（统一三段）

适用人群：想在MR前做共享遗传结构评估、提升研究设计严谨性的同学。

常见错误：把rg显著当作因果证据是常见误区，LDSC仅支持相关性层面判断。

论文写法：论文可在“辅助分析”写LDSC目的、截距结果及其对MR策略的影响。

什么时候用LDSC？

🧬 在做 MR 前，先看暴露与结局是否存在全基因组层面的共享遗传基础。
📉 当多个性状有相关表型时，用 rg 帮助筛选优先进入 MR 的候选。
🧪 用截距（intercept）和回归比值判断群体分层或模型偏差是否明显。

ldsc.py 核心用法（单对性状）

先 munge_sumstats.py 统一格式，再用 ldsc.py --rg 计算遗传相关性。

# 0) 建议环境
conda create -n ldsc python=3.10 -y
conda activate ldsc
git clone https://github.com/bulik/ldsc.git
cd ldsc
pip install -r requirements.txt

# 1) 暴露GWAS标准化（munge）
python ./ldsc/munge_sumstats.py \
 --sumstats ./input/exposure_gwas.txt.gz \
 --N 350000 \
 --snp rsid \
 --a1 effect_allele \
 --a2 other_allele \
 --p p_value \
 --signed-sumstats beta,0 \
 --merge-alleles ./ref/w_hm3.snplist \
 --out ./ldsc_tmp/exposure

# 2) 结局GWAS标准化（munge）
python ./ldsc/munge_sumstats.py \
 --sumstats ./input/outcome_gwas.txt.gz \
 --N 420000 \
 --snp rsid \
 --a1 effect_allele \
 --a2 other_allele \
 --p p_value \
 --signed-sumstats beta,0 \
 --merge-alleles ./ref/w_hm3.snplist \
 --out ./ldsc_tmp/outcome

# 3) 计算遗传相关性 rg
python ./ldsc/ldsc.py \
 --rg ./ldsc_tmp/exposure.sumstats.gz,./ldsc_tmp/outcome.sumstats.gz \
 --ref-ld-chr ./ref/eur_w_ld_chr/ \
 --w-ld-chr ./ref/eur_w_ld_chr/ \
 --out ./ldsc_out/exposure_outcome_rg

# 4) 可选：分别估计h2，检查性状本身可解释遗传度
python ./ldsc/ldsc.py \
 --h2 ./ldsc_tmp/exposure.sumstats.gz \
 --ref-ld-chr ./ref/eur_w_ld_chr/ \
 --w-ld-chr ./ref/eur_w_ld_chr/ \
 --out ./ldsc_out/exposure_h2

批量rg模板（多个代谢物 vs 同一结局）

#!/usr/bin/env bash
set -euo pipefail

mkdir -p ./ldsc_out/batch_rg

OUTCOME="./ldsc_tmp/outcome.sumstats.gz"
REF="./ref/eur_w_ld_chr/"
MET_DIR="./ldsc_tmp/metabolites"  # 已munge完成的代谢物sumstats目录

for FILE in "${MET_DIR}"/*.sumstats.gz; do
 TRAIT=$(basename "${FILE}" .sumstats.gz)
 echo "[RUN] ${TRAIT} vs outcome"

 python ./ldsc/ldsc.py \
  --rg "${FILE},${OUTCOME}" \
  --ref-ld-chr "${REF}" \
  --w-ld-chr "${REF}" \
  --out "./ldsc_out/batch_rg/${TRAIT}_vs_outcome"
done

echo "Batch LDSC rg finished."

结果解读（重点看这几列）

rg：遗传相关性系数，范围约在 -1 到 1；正值表示同向遗传效应。
se：rg 的标准误；越小表示估计越稳定。
z / p：显著性统计量；批量分析建议做 FDR 校正。
h2：若某性状 h2 很低，rg 结果通常不稳定，需谨慎解释。
intercept：明显偏离 1 时，提示可能存在群体分层或其他混杂。
✅ 实战建议：优先关注 |rg|较大 + FDR显著 + intercept合理的性状对。

R脚本：批量读取日志并自动注释

library(data.table)
library(stringr)
library(dplyr)

log_files <- list.files("./ldsc_out/batch_rg", pattern = "\\.log$", full.names = TRUE)

parse_one_log <- function(log_file) {
 lines <- readLines(log_file, warn = FALSE)
 rg_line_id <- grep("^\\./ldsc_tmp/.+\\.sumstats\\.gz\\s+\\./ldsc_tmp/.+\\.sumstats\\.gz", lines)[1]
 if (is.na(rg_line_id)) return(NULL)

 vals <- str_split(str_squish(lines[rg_line_id]), "\\s+")[[1]]
 if (length(vals) < 6) return(NULL)

 data.frame(
  p1 = vals[1],
  p2 = vals[2],
  rg = as.numeric(vals[3]),
  se = as.numeric(vals[4]),
  z = as.numeric(vals[5]),
  p = as.numeric(vals[6]),
  logfile = basename(log_file),
  stringsAsFactors = FALSE
 )
}

rg_table <- bind_rows(lapply(log_files, parse_one_log))
rg_table <- rg_table %>%
 mutate(
  p_fdr = p.adjust(p, method = "fdr"),
  direction = case_when(
   rg > 0 ~ "正相关",
   rg < 0 ~ "负相关",
   TRUE ~ "接近0"
  ),
  strength = case_when(
   abs(rg) >= 0.30 ~ "较强",
   abs(rg) >= 0.10 ~ "中等",
   TRUE ~ "较弱"
  ),
  interpretation = paste0(direction, " / ", strength)
 ) %>%
 arrange(p_fdr)

fwrite(rg_table, "./ldsc_out/batch_rg_summary.tsv", sep = "\t")