GWAS 与群体遗传学基础¶

原始基因型数据（VCF/BED/BGEN）
        ↓
  ① 样本级质控（call rate, 性别检查, 杂合度, 亲缘关系）
        ↓
  ② SNP 级质控（MAF, HWE, call rate, 去重复）
        ↓
  ③ 群体分层校正（PCA / GRM）
        ↓
  ④ 关联分析（GLM / 混合线性模型）
        ↓
  ⑤ 多重检验校正（Bonferroni / FDR / 基因组膨胀因子）
        ↓
  ⑥ 结果可视化与解读（Manhattan图 / QQ图 / 区域图）

# ========== 样本级质控 ==========

# 1. 计算每个样本的缺失率（--mind 0.02 表示去掉缺失率>2%的样本）
plink2 --bfile raw_data \       # 输入PLINK二进制文件（.bed/.bim/.fam）
       --mind 0.02 \            # 去掉缺失率>2%的样本
       --make-bed \             # 输出为PLINK二进制格式
       --out qc_step1           # 输出文件前缀

# 2. 性别检查（X染色体杂合度推断性别，和记录性别对比）
plink2 --bfile qc_step1 \
       --check-sex \            # 推断遗传性别
       --out sex_check          # 输出到 sex_check.sexcheck

# 3. 杂合度检查（过高=样本污染，过低=近亲交配）
plink2 --bfile qc_step1 \
       --het \                  # 计算每个样本的杂合度
       --out het_check          # 输出 het_check.het 文件
# 通常过滤掉偏离均值 ±3SD 的样本

# 4. 亲缘关系检查（去掉重复样本和一级亲属）
plink2 --bfile qc_step1 \
       --make-king-table \      # 计算 KING 亲缘系数
       --king-cutoff 0.0884 \   # 去掉亲缘系数>0.0884的样本对（二级亲属）
       --out kinship             # 输出亲缘关系表

# ========== SNP 级质控 ==========

# 5. 综合 SNP 质控
plink2 --bfile qc_step1 \
       --maf 0.01 \             # 去掉MAF<1%的SNP（太罕见，统计功效低）
       --geno 0.02 \            # 去掉缺失率>2%的SNP
       --hwe 1e-6 \             # 去掉严重偏离HWE的SNP（可能是分型错误）
       --make-bed \
       --out qc_step2

# 6. LD pruning（去掉高度连锁的冗余SNP，用于PCA等）
plink2 --bfile qc_step2 \
       --indep-pairwise 50 5 0.2 \  # 窗口50个SNP，步长5，r²阈值0.2
       --out pruned              # 输出 pruned.prune.in（保留列表）

plink2 --bfile qc_step2 \
       --extract pruned.prune.in \  # 只保留pruning后的独立SNP
       --make-bed \
       --out qc_pruned           # 用于PCA分析

# PCA 分析（用pruning后的独立SNP计算主成分）
plink2 --bfile qc_pruned \
       --pca 20 \               # 提取前20个主成分
       --out pca_result          # 输出 pca_result.eigenvec 和 .eigenval

# pca_result.eigenvec 前10列主成分通常作为协变量加入关联分析
# 用R或Python画PCA散点图，检查群体结构

# ========== PLINK2 GLM（广义线性模型） ==========

# 定量性状关联分析（如BMI、血糖值）
plink2 --bfile qc_step2 \
       --pheno pheno.txt \           # 表型文件（FID IID 表型值）
       --covar pca_result.eigenvec \ # 协变量文件（PCA前10个主成分）
       --covar-col-nums 3-12 \       # 使用第3-12列（PC1-PC10）作为协变量
       --glm hide-covar \            # 线性回归，hide-covar不输出协变量结果
       --out gwas_qt                 # 输出关联结果

# 二分类性状关联分析（如病例/对照）
plink2 --bfile qc_step2 \
       --pheno pheno_cc.txt \        # 表型文件（1=对照，2=病例）
       --covar pca_result.eigenvec \
       --covar-col-nums 3-12 \
       --1 \                         # 表型编码：0=对照，1=病例
       --glm firth-fallback hide-covar \ # 逻辑回归+Firth校正（小样本更准）
       --out gwas_cc                  # 输出关联结果

# ========== REGENIE（两步法，适合大样本如UK Biobank） ==========

# Step 1：全基因组回归拟合（构建预测模型，校正群体结构）
regenie \
  --step 1 \                    # 第一步：全基因组回归
  --bed qc_step2 \              # 输入PLINK格式
  --phenoFile pheno.txt \       # 表型文件
  --covarFile covar.txt \       # 协变量文件
  --bt \                        # 二分类性状（--qt 用于定量性状）
  --bsize 1000 \                # block大小（每次处理1000个SNP）
  --lowmem \                    # 低内存模式
  --out step1_result            # 输出预测结果

# Step 2：逐个SNP关联检验
regenie \
  --step 2 \                    # 第二步：关联检验
  --bed qc_step2 \              # 输入基因型
  --phenoFile pheno.txt \
  --covarFile covar.txt \
  --bt \                        # 二分类性状
  --firth --approx \            # Firth 近似校正（推荐用于二分类）
  --pred step1_result_pred.list \ # 第一步的预测结果
  --bsize 200 \                 # 第二步的 block 大小
  --minMAC 5 \                  # 最小等位基因计数 ≥ 5
  --out step2_result            # 输出最终关联结果

# ========== GCTA fastGWA（混合线性模型，适合大样本） ==========

# 1. 构建稀疏GRM（遗传关系矩阵）
gcta64 --bfile qc_step2 \
       --make-grm --out grm_all    # 构建完整GRM

# 也可用稀疏GRM加速
gcta64 --grm grm_all \
       --make-bK-sparse 0.05 \     # 稀疏化（阈值0.05）
       --out sp_grm                # 稀疏GRM

# 2. 运行 fastGWA 关联分析
gcta64 --bfile qc_step2 \
       --grm-sparse sp_grm \       # 使用稀疏GRM
       --fastGWA-mlm \             # 混合线性模型关联分析
       --pheno pheno.txt \
       --qcovar covar_quant.txt \  # 定量协变量（如年龄、PC1-10）
       --covar covar_cat.txt \     # 分类协变量（如性别）
       --threads 10 \              # 使用10个线程
       --out fastgwa_result        # 输出结果

# PLINK2 多重检验校正
plink2 --adjust-file gwas_cc.PHENO1.glm.logistic.hybrid \  # 输入关联结果
       --adjust-gc \         # 基因组控制校正（genomic control）
       --out adjusted        # 输出校正后的 p 值

# Manhattan图和QQ图（Python + matplotlib）
import pandas as pd                    # 读取数据
import matplotlib.pyplot as plt        # 绑图
import numpy as np                     # 数学计算

# 读取 PLINK2 关联结果
df = pd.read_csv('gwas_cc.PHENO1.glm.logistic.hybrid',
                 sep='\t',              # Tab 分隔
                 usecols=['#CHROM', 'POS', 'P'])  # 只读需要的列
df.columns = ['CHR', 'POS', 'P']       # 重命名列
df = df.dropna(subset=['P'])            # 去掉缺失 p 值
df['-log10P'] = -np.log10(df['P'])      # 计算 -log10(p)

# --- Manhattan 图 ---
fig, ax = plt.subplots(figsize=(16, 6))
colors = ['#4C72B0', '#DD8452']         # 交替颜色（奇数/偶数染色体）
x_offset = 0                            # X轴累加偏移

for chrom in range(1, 23):              # 逐染色体画图
    subset = df[df['CHR'] == chrom]
    color = colors[chrom % 2]
    ax.scatter(subset['POS'] + x_offset,     # X = 位置 + 偏移
               subset['-log10P'],            # Y = -log10(p)
               c=color, s=2, alpha=0.6)
    x_offset += subset['POS'].max() if len(subset) > 0 else 0

ax.axhline(y=-np.log10(5e-8), color='red',    # 全基因组显著线
           linestyle='--', label='p=5e-8')
ax.axhline(y=-np.log10(1e-5), color='blue',   # 提示性显著线
           linestyle=':', label='p=1e-5')
ax.set_xlabel('Chromosome')
ax.set_ylabel('-log10(P)')
ax.set_title('Manhattan Plot')
ax.legend()
plt.tight_layout()
plt.savefig('manhattan.png', dpi=300)   # 保存高分辨率图片
plt.show()

# --- QQ 图 ---
observed = np.sort(df['-log10P'])[::-1]           # 实际 p 值（降序）
expected = -np.log10(np.arange(1, len(observed)+1)
                     / (len(observed)+1))         # 期望 p 值

fig, ax = plt.subplots(figsize=(6, 6))
ax.scatter(expected, observed, s=2, alpha=0.5)     # 散点图
ax.plot([0, max(expected)], [0, max(expected)],     # 对角线
        color='red', linestyle='--')
ax.set_xlabel('Expected -log10(P)')
ax.set_ylabel('Observed -log10(P)')
ax.set_title('QQ Plot')

# 计算 genomic inflation factor λ
chisq = df['P'].apply(lambda p: (-2) * np.log(p) if p > 0 else np.nan)
# 更标准的方式：用 scipy
from scipy.stats import chi2
lambda_gc = np.median(chi2.ppf(1 - df['P'].dropna(), 1)) / chi2.ppf(0.5, 1)
ax.text(0.05, 0.95, f'λ = {lambda_gc:.3f}',
        transform=ax.transAxes, fontsize=12)
plt.tight_layout()
plt.savefig('qq_plot.png', dpi=300)
plt.show()

# PLINK2 计算 Fst
# --pheno 文件中群体标签作为表型
plink2 --bfile qc_step2 \
       --fst POP \             # 按 POP 列计算群体间 Fst
       --pheno pop_labels.txt \  # 群体标签文件
       --out fst_result         # 输出 Fst 结果

# GCTA 计算 Fst
gcta64 --bfile qc_step2 \
       --pca 20 \              # 先做PCA
       --out pca_fst            # PCA结果可辅助Fst计算

# PLINK2 PCA（推荐使用 LD pruning 后的 SNP）
plink2 --bfile qc_pruned \
       --pca 20 \              # 提取前20个主成分
       --out population_pca    # 输出 .eigenvec（主成分坐标）和 .eigenval（特征值）

# Hardy-Weinberg 平衡
p² + 2pq + q² = 1          # p, q 是两个等位基因的频率

# 基因组膨胀因子
λ = median(χ²_observed) / 0.4549  # 0.4549 = χ²(df=1) 的中位数

# Fst（Weir & Cockerham 估计）
Fst = (HT - HS) / HT       # HT=总杂合度，HS=亚群平均杂合度

# 遗传力（广义）
H² = V_G / V_P              # 遗传方差 / 表型方差

# SNP遗传力
h²_SNP = V_SNP / V_P        # SNP解释的方差 / 表型方差

1. 先理解核心概念（SNP/MAF/LD/HWE/群体分层）
2. 用 PLINK2 跑一遍完整的模拟数据质控流程
3. 学会画 Manhattan 图和 QQ 图
4. 了解混合线性模型的原理（GCTA/REGENIE）
5. 阅读 T2D GWAS 的代表性论文
6. 将 GWAS 思路和该宏基因组项目做对比