615_全基因组关联分析GWAS¶

# PLINK 1.9（经典GWAS工具）
conda install -c bioconda plink

# PLINK 2.0（更快，推荐新项目用）
conda install -c bioconda plink2

# GCTA（遗传力估计+GWAS）
conda install -c bioconda gcta

# 验证安装
plink --version
plink2 --version
gcta --version

#!/bin/bash
# GWAS质控流程——每步都不能跳！

INPUT="raw_gwas"   # 原始PLINK文件前缀
QC="gwas_qc"       # QC输出目录
mkdir -p ${QC}

# === 第1步：个体级QC ===
# 去除缺失率高的个体
plink --bfile ${INPUT} \
    --mind 0.02 \              # 个体缺失率>2%的移除
    --make-bed \
    --out ${QC}/step1

# === 第2步：SNP级QC ===
plink --bfile ${QC}/step1 \
    --geno 0.02 \              # SNP缺失率>2%的移除
    --maf 0.01 \               # MAF<1%的移除
    --hwe 1e-6 \               # 偏离HWE的移除（用对照组检验）
    --make-bed \
    --out ${QC}/step2

# === 第3步：检查性别 ===
plink --bfile ${QC}/step2 \
    --check-sex \              # 检查遗传性别与记录是否一致
    --out ${QC}/sex_check
# 查看不一致的个体
awk '$3 != $4 && $4 != 0' ${QC}/sex_check.sexcheck

# === 第4步：去除亲缘关系个体 ===
plink --bfile ${QC}/step2 \
    --genome \                 # 计算IBD
    --min 0.2 \                # 只输出PI_HAT>0.2的对
    --out ${QC}/ibd
# 去除一对中缺失率高的那个
awk '$10 > 0.2 {print $3, $4}' ${QC}/ibd.genome > ${QC}/remove_related.txt

plink --bfile ${QC}/step2 \
    --remove ${QC}/remove_related.txt \
    --make-bed \
    --out ${QC}/step3

# === 第5步：群体分层校正（PCA）===
# 先LD pruning
plink --bfile ${QC}/step3 \
    --indep-pairwise 50 5 0.2 \
    --out ${QC}/ld
plink --bfile ${QC}/step3 \
    --extract ${QC}/ld.prune.in \
    --pca 10 \                 # 计算前10个PC
    --out ${QC}/pca

# 最终QC完成的数据
cp ${QC}/step3.* ${QC}/gwas_ready.*
echo "QC完成！最终样本数: $(wc -l < ${QC}/gwas_ready.fam)"
echo "最终SNP数: $(wc -l < ${QC}/gwas_ready.bim)"

# === 病例-对照研究（二元表型）===
plink --bfile ${QC}/gwas_ready \
    --logistic \               # logistic回归（二元表型用）
    --covar ${QC}/pca.eigenvec \  # 前10个PC作为协变量
    --covar-number 3-12 \      # 使用第3-12列（PC1-PC10）
    --adjust \                 # 多重检验校正
    --ci 0.95 \                # 95%置信区间
    --out gwas_result          # 输出前缀

# === 数量性状（连续表型）===
plink --bfile ${QC}/gwas_ready \
    --linear \                 # 线性回归（连续表型用）
    --pheno phenotype.txt \    # 外部表型文件
    --covar ${QC}/pca.eigenvec \
    --covar-number 3-12 \
    --out gwas_quantitative

# === PLINK2更快的方式 ===
plink2 --bfile ${QC}/gwas_ready \
    --glm \                    # 广义线性模型
    --covar ${QC}/pca.eigenvec \
    --covar-col-nums 3-12 \
    --out gwas_plink2

# === 第1步：构建GRM（遗传关系矩阵）===
gcta --bfile ${QC}/gwas_ready \
    --make-grm \               # 构建GRM
    --thread-num 16 \          # 线程数
    --out grm                  # 输出前缀

# === 第2步：稀疏GRM（大样本更快）===
gcta --grm grm \
    --make-bK-sparse 0.05 \    # 稀疏化
    --out sparse_grm

# === 第3步：运行fastGWA ===
gcta --bfile ${QC}/gwas_ready \
    --grm-sparse sparse_grm \  # 使用稀疏GRM
    --fastGWA-mlm \            # 混合线性模型
    --pheno phenotype.txt \    # 表型文件
    --qcovar ${QC}/pca.eigenvec \  # 数量协变量
    --thread-num 16 \
    --out gwas_gcta

# === 安装和加载包 ===
# install.packages("qqman")
library(qqman)  # GWAS可视化专用包

# === 读取GWAS结果 ===
gwas <- read.table("gwas_result.assoc.logistic", header = TRUE)
# 只保留ADD（加性模型）的结果
gwas <- gwas[gwas$TEST == "ADD", ]

# === Manhattan Plot ===
pdf("manhattan_plot.pdf", width = 14, height = 6)
manhattan(
    gwas,
    chr = "CHR",       # 染色体列
    bp = "BP",         # 位置列
    p = "P",           # p值列
    snp = "SNP",       # SNP名列
    col = c("steelblue", "orange"),  # 交替颜色
    suggestiveline = -log10(1e-5),   # 暗示性阈值
    genomewideline = -log10(5e-8),   # 全基因组显著阈值
    main = "GWAS Manhattan Plot"     # 标题
)
dev.off()

# === QQ Plot ===
pdf("qq_plot.pdf", width = 7, height = 7)
qq(gwas$P, main = "GWAS QQ Plot")  # QQ图
dev.off()

# === 计算lambda值（基因组膨胀因子）===
chisq <- qchisq(1 - gwas$P, 1)  # 转为卡方统计量
lambda <- median(chisq) / qchisq(0.5, 1)  # lambda值
cat("Genomic inflation factor (lambda):", lambda, "\n")
# lambda应接近1.0，>1.1说明有群体分层问题

import gwaslab as gl  # 导入gwaslab
import pandas as pd   # 导入pandas

# 读取GWAS结果
df = pd.read_csv("gwas_result.assoc.logistic", sep="\s+")  # 读取
df = df[df["TEST"] == "ADD"]  # 只保留加性模型

# 创建gwaslab对象
mysumstats = gl.Sumstats(
    df,
    snpid="SNP",   # SNP ID列
    chrom="CHR",   # 染色体列
    pos="BP",      # 位置列
    p="P",         # p值列
    beta="BETA",   # 效应值列
    verbose=False
)

# Manhattan plot + QQ plot
mysumstats.plot_mqq()  # 一键画曼哈顿图+QQ图

# 用GCTA估计SNP遗传力
gcta --grm grm \
    --reml \                   # REML方法
    --pheno phenotype.txt \    # 表型文件
    --qcovar ${QC}/pca.eigenvec \
    --thread-num 16 \
    --out heritability

# 查看结果
cat heritability.hsq
# V(G)/Vp 就是SNP遗传力（h²_SNP）

# 安装ldsc
pip install ldsc

# 准备summary statistics
# 用ldsc格式化GWAS结果
python munge_sumstats.py \
    --sumstats gwas_result.assoc.logistic \
    --out gwas_formatted \
    --merge-alleles w_hm3.snplist  # HapMap3 SNP列表

# 运行LD score regression
python ldsc.py \
    --h2 gwas_formatted.sumstats.gz \  # GWAS结果
    --ref-ld-chr eur_w_ld_chr/ \       # LD分数参考
    --w-ld-chr eur_w_ld_chr/ \
    --out ldsc_result

# intercept接近1.0说明没有膨胀
# h2就是遗传力估计