全基因组关联分析（GWAS）

**全基因组关联分析（Genome-Wide Association Study, GWAS）**是一种利用统计学方法研究基因变异（通常是单核苷酸多态性，SNPs）与特定性状或疾病之间关联的分析方法。以下是对GWAS的介绍，以及分析流程和示例代码。

---

GWAS的基本原理

GWAS是一种无假设驱动的研究方法，主要目标是通过分析整个基因组范围内的遗传变异与表型之间的关系，识别与疾病或性状相关的基因位点。研究通常涉及以下步骤：

1. 基因型数据和表型数据收集： 获得个体的基因型数据（例如从全基因组测序或基因芯片）和表型数据。

2. 数据质量控制： 对基因型数据和样本进行质控（如排除低质量SNP和样本）。

3. 单位点关联分析： 对每个SNP进行统计测试，通常使用线性或逻辑回归分析。

4. 多重检验校正： 控制假阳性率（如Bonferroni校正或FDR）。

5. 结果可视化： 生成曼哈顿图和QQ图。

6. 功能注释和生物学解释： 分析显著关联的SNP所在的基因或区域，并探索潜在机制。

---

GWAS分析步骤

1. 数据准备

• 基因型数据：PLINK格式（bed/bim/fam）。

• 表型数据：通常是一个包含样本ID和性状值的表格文件。

2. 数据质控

• 去除低质量的SNP和样本。

• 确保基因型呼叫率 > 95%。

• 剔除等位基因频率（MAF）过低的SNP。

• 检查哈迪-温伯格平衡（HWE）。

3. 单位点关联分析

• 对连续性状：线性回归。

• 对二分类性状：逻辑回归。

4. 多重检验校正

• Bonferroni校正：显著性水平 =α/总SNP数\alpha / \text{总SNP数}α/总SNP数。

• FDR：调整p值以控制假发现率。

5. 结果可视化

• 曼哈顿图：展示每个染色体的SNP显著性。

• QQ图：检查p值的分布，评估系统性偏倚。

GWAS分析通常涉及从数据处理到可视化的一整套工具链。以下是主流的GWAS工具分类及其功能，每种工具附有简单的分析代码。

---

1. 数据质控与基因型操作工具

1.1 PLINK

PLINK 是经典的工具，广泛用于基因型数据的质控、格式转换和初步分析。

• 功能： 数据质控（过滤MAF、HWE、呼叫率）、关联分析、数据转换等。

安装：

# 下载和安装PLINKwget http://zzz.bwh.harvard.edu/plink/plink_linux_x86_64.zipunzip plink_linux_x86_64.zipchmod +x plinkmv plink /usr/local/bin/plink

质控代码：​​​​​​​

# Step 1: 过滤呼叫率低的SNP和样本plink --bfile input_data --geno 0.05 --mind 0.05 --make-bed --out qc_step1
# Step 2: 过滤MAF低于0.01的SNPplink --bfile qc_step1 --maf 0.01 --make-bed --out qc_step2
# Step 3: 过滤HWE显著偏离的SNP（P < 1e-6）plink --bfile qc_step2 --hwe 1e-6 --make-bed --out qc_final

关联分析：​​​​​​​

# 连续性状分析plink --bfile qc_final --pheno phenotype.txt --linear --out gwas_linear
# 二分类性状分析plink --bfile qc_final --pheno phenotype.txt --logistic --out gwas_logistic

1.2 BCFtools

BCFtools 专注于处理和操作VCF/BCF格式的基因型数据。

• 功能： VCF格式转换、样本筛选、变异过滤等。

安装：​​​​​​​

# 安装sudo apt-get install bcftools

代码示例：​​​​​​​

# 筛选MAF > 0.01的变异bcftools view -i 'MAF>0.01' input.vcf > filtered.vcf
# 转换VCF为PLINK格式bcftools convert --plink input.vcf --out converted

2. 高效关联分析工具

2.1 GEMMA

GEMMA 是针对混合线性模型（MLM）优化的工具，适合分析复杂性状和控制遗传结构。

功能：

• 混合线性模型（MLM）。

• 支持遗传相关矩阵（GRM）和环境因子。

安装：​​​​​​​

# 下载和安装GEMMAwget https://github.com/genetics-statistics/GEMMA/releases/download/v0.98.1/gemma-0.98.1-linux-static.zipunzip gemma-0.98.1-linux-static.zipchmod +x gemmamv gemma /usr/local/bin/

分析代码：​​​​​​​

# Step 1: 生成遗传相关矩阵（GRM）gemma -bfile input_data -gk 1 -o kinship
# Step 2: 运行GWAS关联分析（线性混合模型）gemma -bfile input_data -k output/kinship.sXX.txt -lm 4 -outdir output

2.2 SAIGE

SAIGE 是针对稀有变异优化的大规模GWAS工具，适合处理大规模样本。

功能：

• 稀有变异关联分析。

• 支持二分类性状和连续性状。

安装与代码：SAIGE需要配置R环境和依赖项，请参考官方文档：SAIGE GitHub。

---

2.3 BOLT-LMM

BOLT-LMM 是大规模数据集下的线性混合模型工具，适合处理人类表型与基因型数据。

功能：

• 快速线性混合模型分析。

• 支持多线程和大规模计算。

安装与代码：请参考：BOLT-LMM官方页面。分析代码类似GEMMA。

---

3. 可视化与结果解释工具

3.1 LocusZoom

LocusZoom 是专门用于绘制区域关联图（regional association plot）的工具。

安装：​​​​​​​

# 安装Python版LocusZoompip install locuszoom

代码示例：​​​​​​​

from locuszoom import plot
# 创建一个局部关联图plot("gwas_results.assoc.linear", region="chr1:1000000-2000000", pval_col="P")

3.2 R语言+ggplot2

R语言是灵活的可视化工具，可以绘制曼哈顿图和QQ图。

安装：​​​​​​​

# 安装ggplot2install.packages("ggplot2")

代码示例：​​​​​​​

# 加载数据gwas_data <- read.table("gwas_results.assoc.linear", header=TRUE)
# 曼哈顿图library(ggplot2)ggplot(gwas_data, aes(x=BP, y=-log10(P), color=as.factor(CHR))) +    geom_point() +    theme_minimal() +    geom_hline(yintercept=-log10(5e-8), color="red", linetype="dashed") +    labs(x="Genomic Position", y="-log10(p-value)", title="Manhattan Plot")

3.3 FUMA

FUMA 是一个在线工具，用于解释GWAS结果，包括基因注释、功能富集分析和路径分析。

• 官网：FUMA

---

4. 功能注释与后续分析工具

4.1 ANNOVAR

ANNOVAR 是基因功能注释工具，用于关联SNP与基因功能。

安装与代码：​​​​​​​

# 注释VCF文件annotate_variation.pl -out output -build hg19 input.vcf humandb/

4.2 Ensembl VEP

Ensembl VEP 是另一个注释工具，可预测SNP的功能影响。

安装与代码：参考：VEP官网

---

完整分析流程工具链

1. 数据质控： PLINK 或 BCFtools。

2. 关联分析： GEMMA, SAIGE, 或 BOLT-LMM。

3. 可视化： LocusZoom, ggplot2。

4. 功能注释： ANNOVAR, VEP。

5. 结果解释： FUMA, DAVID（功能富集分析）。

使用Python生成曼哈顿图和QQ图。​​​​​​​

import pandas as pdimport numpy as npimport matplotlib.pyplot as pltimport seaborn as sns
# 读取GWAS结果gwas_data = pd.read_csv("gwas_results.assoc.linear", delim_whitespace=True)
# 曼哈顿图def plot_manhattan(data, threshold=5e-8):    data['-log10(p)'] = -np.log10(data['P'])    data['Chromosome'] = data['CHR'].astype('category')    data['Chromosome'] = data['Chromosome'].cat.set_categories(        sorted(data['Chromosome'].unique(), key=lambda x: int(x)))
    plt.figure(figsize=(12, 6))    sns.scatterplot(data=data, x='BP', y='-log10(p)', hue='Chromosome', palette='tab20', legend=False)    plt.axhline(-np.log10(threshold), color='red', linestyle='--', label='Genome-wide significance')    plt.xlabel('Genomic Position')    plt.ylabel('-log10(p-value)')    plt.title('Manhattan Plot')    plt.legend()    plt.show()
plot_manhattan(gwas_data)
# QQ图def plot_qq(data):    observed = -np.log10(data['P'].sort_values())    expected = -np.log10(np.linspace(1/len(data), 1, len(data)))    plt.figure(figsize=(6, 6))    plt.scatter(expected, observed, edgecolor="k", alpha=0.7)    plt.plot([0, max(expected)], [0, max(expected)], linestyle='--', color='red')    plt.xlabel('Expected -log10(p)')    plt.ylabel('Observed -log10(p)')    plt.title('QQ Plot')    plt.show()
plot_qq(gwas_data)

生信大白记第36记，就到这里，关注我！

下一记，持续更新学习生物信息学的内容！

生信大白记邮箱账号：shengxindabaiji@163.com

生信大白记简书账号：生信大白记

生信大白记CSDN账号：生信大白记

生信大白记微信公众号：生信大白记

加入生信大白记交流群938339543