Gemma自动化脚本

GEMMA 自动化 GWAS 分析流程

目录

简介

本项目是一个自动化的GWAS（全基因组关联分析）流程脚本，它将从原始的VCF文件和表型数据开始，一直到生成最终的显著性位点列表，极大地简化了使用PLINK和GEMMA进行分析的复杂步骤。

该脚本特别设计用于处理多个性状，并支持二分类性状（如抗病/感病）和数量性状（如株高、产量）。它内置了数据质控、群体结构校正（PCA）、多种GWAS模型选择和P值校正等关键功能。

核心功能

  • 端到端自动化: 从VCF文件一键运行到GWAS结果，无需手动分步操作。
  • 多表型支持: 自动为表型文件中的每一个性状（每一列）独立运行完整的GWAS分析。
  • 智能质控 (QC): 使用PLINK对基因型数据进行标准化质控，包括MAF、HWE和缺失率过滤。
  • 灵活的模型选择:
    • LMM (线性混合模型): 适用于校正群体结构和个体间亲缘关系，是GWAS的标准模型。
    • BSLMM (贝叶斯稀疏线性混合模型): 一种更高级的模型，可以同时估计SNP的固定效应和随机效应，适用于复杂性状。
  • 群体结构校正: 可选地使用PCA（主成分分析）的结果作为协变量，以校正在LMM中由群体分层引起的假阳性。
  • 多种P值校正: 支持严格的Bonferroni校正、控制假阳性的FDR (Benjamini-Hochberg)校正，或不进行校正（使用名义P值）。
  • 智能表型处理: 自动将二分类表型（0/1编码）转换为PLINK和GEMMA所需的1/2编码（对照/病例），避免常见错误。
  • 详细日志与报告: 生成详细的运行日志 (gwas_analysis.log) 和一份最终的分析总结报告 (analysis_summary.txt)，便于溯源和检查。

环境要求

在运行此脚本前，请确保您的系统环境中已安装并配置好以下软件和库：

  1. Python 3: 脚本运行环境。

  2. Python库:

    • pandas: 用于数据处理。
    • statsmodels: 用于FDR校正。如果未安装，脚本会提示并退出。
    1
    pip install pandas statsmodels

  3. PLINK (v1.9或更高版本): 用于VCF格式转换和数据质控。必须将其可执行文件路径添加到系统的PATH环境变量中。

  4. GEMMA: 核心GWAS分析软件。必须将其可执行文件路径添加到系统的PATH环境变量中。

安装与准备

软件安装

请根据上述要求，确保plink和gemma命令可以在您的终端中直接运行。

输入文件准备

您需要准备两个核心输入文件：

  1. VCF文件 (--vcf):

    • 标准的VCF格式文件，包含了您群体的基因型数据。

  2. 表型文件 (--pheno):

    • 必须是制表符分隔（tab-separated）的文本文件。
    • 必须包含表头 (header)。
    • 第一列必须是样本ID，且该列的表头名称不限（脚本会自动识别）。样本ID必须与VCF文件中的样本ID一致。
    • 从第二列开始，每一列代表一个性状。
    • 对于二分类性状: 请使用 0代表对照组(control/healthy)，1代表病例组(case/affected)。脚本会自动将其转换为PLINK所需的 1 (对照) 和 2 (病例)。
    • 缺失值可以留空或使用NA。

表型文件示例 (phenotypes.txt):

1
2
3
4
5
sample	resistance_A	height_cm	yield_kg
sample01	1	120.5	2.5
sample02	0	98.2	1.8
sample03	1	115.0	NA
...

使用方法

快速开始

以下是一个典型的运行示例，分析一个二分类性状，并使用PCA校正群体结构和FDR进行P值校正。

1
2
3
4
5
6
7
8
9
10
python gwas_pipeline.py \
  --vcf genotypes.vcf \
  --pheno phenotypes.txt \
  --output ./gwas_results_resistance \
  --pheno-type binary \
  --model lmm \
  --use-pca \
  --pca-components 5 \
  --correction fdr \
  --sig-threshold 0.05

参数详解

参数 是否必需 描述 默认值
--vcf 是 输入的VCF格式基因型文件。 -
--pheno 是 制表符分隔的表型文件。第一列为样本ID，后续列为性状数据。 -
--output 否 所有输出文件的存放目录。 gwas_output
--maf 否 最小等位基因频率（Minor Allele Frequency）过滤阈值。 0.05
--hwe 否 哈迪-温伯格平衡检验（Hardy-Weinberg Equilibrium）的p值过滤阈值。 0.001
--missing 否 SNP和样本允许的最大缺失率。0.1表示保留缺失率低于10%的SNP和样本。 0.9
--sig-threshold 否 显著性水平阈值。用于FDR的目标q值，或Bonferroni校正的名义alpha水平。 0.05
--use-pca 否 设置此标志以运行PCA并将主成分作为协变量（仅对LMM模型有效）。 False
--pca-components 否 当--use-pca启用时，指定使用的主成分数量。 3
--model 否 选择GWAS分析模型。 lmm
--correction 否 P值校正方法。bonferroni: 严格控制；fdr: 控制假阳性率；none: 不校正。 none
--pheno-type 否 表型的数据类型。binary用于病例/对照研究，quantitative用于连续值性状。 quantitative

工作流程详解

脚本将自动按以下顺序执行分析：

  1. VCF转PLINK: 使用PLINK将输入的VCF文件转换为二进制的bed/bim/fam格式。
  2. 质量控制 (QC): 依次使用PLINK进行三步过滤：
    • 基于MAF (--maf)
    • 基于HWE p值 (--hwe)
    • 基于SNP和样本的缺失率 (--missing)
  3. 表型文件准备:
    • 读取PLINK质控后的样本列表（.fam文件）。
    • 将用户提供的表型文件与样本列表对齐，确保顺序和成员一致。
    • 如果是二分类性状，自动将 0/1 转换为 1/2。
    • 生成可用于后续分析的多表型文件。
  4. 循环执行GWAS: 对表型文件中的每一个性状，执行以下子流程：
    a. 准备数据: 为当前性状创建一个专用的PLINK文件集。
    b. 计算亲缘关系矩阵: 使用GEMMA计算标准化相关亲缘关系矩阵（Kinship Matrix）。
    c. (可选) PCA分析: 如果设置了--use-pca且模型为lmm，使用PLINK计算PCA，并提取指定数量的主成分作为协变量文件。
    d. 运行GWAS: 调用GEMMA执行核心分析。根据--model和--pheno-type选择合适的参数（如 LMM-Wald, LMM-LRT, 或 BSLMM）。
    e. 筛选显著SNP:
    - **LMM模型**: 根据指定的`--correction`方法（`bonferroni`, `fdr`, `none`）和`--sig-threshold`筛选显著关联的SNP。
- **BSLMM模型**: 筛选`gamma`值（后验包含概率）大于`0.01`的SNP。
    f. 保存结果: 将筛选出的显著SNP信息保存到单独的结果文件中。
  5. 生成总结报告: 所有性状分析完成后，生成一份包含所有参数和输入/输出文件路径的总结报告。

输出文件说明

所有输出文件都将保存在您通过--output指定的目录中。主要包括：

  • 日志和报告:

    • gwas_analysis.log: 详细的运行日志，记录了每一步的命令和状态。
    • analysis_summary.txt: 最终的分析总结报告。

  • 对每个性状 [trait_name] 的主要结果:

    • [trait_name]_significant_snps_[correction_method].txt: 最重要的结果文件。包含了通过显著性阈值筛选后的SNP列表及其相关统计信息。
    • output/[trait_name]_gwas.assoc.txt: (LMM模型) GEMMA生成的完整GWAS结果。
    • output/[trait_name]_gwas.param.txt: (BSLMM模型) GEMMA生成的参数估计结果，包含gamma值。
    • output/[trait_name]_gwas.log.txt: GEMMA为该性状分析生成的日志。

  • 关键中间文件:

    • gwas_data_qc.bed/bim/fam: 经过完整质控后的PLINK格式基因型文件。
    • output/kinship_[trait_name].cXX.txt: 为每个性状计算的亲缘关系矩阵。
    • pca_covariates_[trait_name].txt: (如果使用PCA) 为每个性状生成的主成分协变量文件。

重要注意事项

  • PATH环境变量: 请务必确保plink和gemma命令能在任何路径下被系统找到。这是脚本成功运行的前提。
  • 表型文件格式: 严格遵守输入文件格式要求，特别是制表符分隔和表头。对于二分类性状，务必使用0代表对照组，1代表病例组。
  • --missing参数: 此参数直接传递给PLINK的--geno和--mind选项，代表允许的最大缺失率。默认值0.9是一个非常宽松的设置（允许高达90%的缺失），适用于质量较差的数据。对于高质量数据，您可能希望使用更严格的阈值，如--missing 0.1。
  • BSLMM与PCA: BSLMM模型本身可以很好地处理群体结构，因此脚本设计上在bslmm模型下会忽略--use-pca设置，即不使用PCA作为协变量。
  • 计算资源: GWAS分析，特别是亲缘关系矩阵的计算，可能需要较多的内存和计算时间，具体取决于您的样本量和SNP数量。
生物信息学
#生物信息学
Gemma自动化脚本
https://lixiang117423.github.io/article/gemma/
作者
李详【Xiang LI】
发布于
2025年7月11日
许可协议