软件gemma官方文档

GWAS 分析软件 GEMMA 官方文档笔记。

GEMMA User Manual

简介

缺失数据

基因型缺失

基因型建议先填充。如果不填充，首先会剔除缺失率较高的位点（默认 5%），之后将缺失位点填充为该位点在群体内的基因型均值。

表型缺失

在LMM或BSLMM分析中，会剔除缺失表型的个体，但是所有的个体都会用于构建关系矩阵。但是在构建 G阵的过程中，缺失率和MAF都只会基于分析个体（表型不为缺失的个体）进行计算，如果所有个体的表型均为缺失，此时估计的G阵元素会全是 “nan” 。

输入文件格式

GEMMA 要求4个主要的输入文件，基因型文件、表型文件、亲缘关系文件还有协变量文件。基因型和表型文件可以是plink 格式也可以是BIMBAM格式。

PLINK 的二进制格式

fam文件 GEMMA只读取第二列个体号和第六列表型；通过 -n 改变读取表型的列（-n 1 表示还是读取第六列；-n 2 表示读取第七列，以此类推）

GEMMA将plink里的位点编码为1/0，一般是将bim文件中的第5列编码为1（minor allele），将第6列编码为0（major allele）。（看版本，好像不同版本1/0编码可能是相反的，也就是剂量效应用的等位基因是不同的）。

GEMMA会将表型中的 -9 或者 NA 视为缺失值。如果fam文件中的表型是疾病，推荐的编码方式是control记为0，case记为1。

关系矩阵文件格式

关系矩阵有两种形式，第一种是原始的，第二种关系矩阵的特征值和特征向量。

原始关系矩阵

第一种就是n*n的矩阵形式，行和列的顺序与fam文件相同。

第二种就是三列的形式（个体号-个体号-值，个体号应该是fam文件的第二列），个体号的顺序不需要和fam文件一致，未列出的值均为0。（文档没说，但应该是上三角矩阵）E

plink文件格式的可以通过 -km 1 和 -km 2选择关系矩阵形式（应该1表示原始的矩阵形式，3表示3列的）

协变量文件（可选）

如果没有协变量，GEMMA 默认会有一个截距项（就是均值）；但是如果存在协变量文件，那么GEMMA 就不会自带一个截距项了。

如果有协变量文件，协变量第一列必须是全为1的截距项，因为这个时候 GEMMA 不会自动提供截距项，所以你要自己提供。

GSL ERROR: matrix is singular in lu.c at line 266 errno 1

如果你的某个协变量的值和部分SNP的基因型值完全相同（指乘以一个固定常数后相同），那么就会造成优化算法出错，产生 GSL errors 。为了避免这个问题，你要么先对协变量进行回归分析，然后用残差作为新的表型；要么就剔除那些与协变量相同的snp位点（采用协变量作为表型，基因型数据做为特征）。

It can happen, especially in a small GWAS data set, that some of the covariates will be identical to some of the genotypes (up to a scaling factor). This can cause problems in the optimization algorithm and evoke GSL errors. To avoid this, one can either regress the phenotypes on the covariates and use the residuals as new phenotypes, or use only SNPs that are not identical to any of the covariates for the analysis. The later can be achieved, for example, by performing a standard linear regression in the genotype data, but with covariates as phenotypes.

实际运行发现还有另一个情况也会导致这个问题，就是如果一个协变量完全没有变异，也会导致这个报错。

运行GEMMA

构建关系矩阵

./gemma -bfile [prefix] -gk [num] -o [prefix]

-bfile 二进制plink文件前缀

-gk 表示构建哪种关系矩阵

• -gk 1 计算中心化的亲缘关系矩阵 (the centered relatedness matrix)
• -gk 2 计算标准化的亲缘关系矩阵 (the standardized relatedness matrix)

-o 输出G阵前缀

两种G阵

这里我们定义 \(X\) 为 \(n \times q\) 的基因型矩阵，\(n\) 为样本数，\(q\) 为SNP数目（就是常规的设计矩阵，一行一个样本，一列一个特征）。

\(x_{i}\) 表示第 \(i\) 列SNP的列向量，\(\bar{x_{i}}\) 表示第 \(i\) 列SNP的群体均值，\(v_{x_{i}}\) 表示第 \(i\) 列SNP的群体均值，\(1_{n}\) 表示一个 \(n \times 1\) 的全为1的列向量。

因此两种矩阵的计算方式为 (文章中给的公式)： \[ \begin{aligned} G_{c} &=\frac{1}{p} \sum_{i=1}^{p}\left(\mathrm{x}_{i}-1_{n} \bar{x}_{i}\right)\left(\mathrm{x}_{i}-1_{n} \bar{x}_{i}\right)^{T}, \\ G_{s} &=\frac{1}{p} \sum_{i=1}^{p} \frac{1}{v_{x_{i}}}\left(\mathrm{x}_{i}-1_{n} \overline{x_{i}}\right)\left(\mathrm{x}_{i}-1_{n} \overline{x_{i}}\right)^{T} . \end{aligned} \] 如果换成我们常见的G阵构建方式，如果 \(Z_{1}\) 表示对 \(X\) 每列进行中心化之后的矩阵，\(Z_{2}\) 表示对 \(X\) 每列进行标准化之后的矩阵。上面的式子可以写为下式： \[ \mathbf{G_{c}}=\frac{1}{p} \mathbf{Z_{1}} \mathbf{Z}_{1}^{\top} \]

\[ \mathbf{G_{s}}=\frac{1}{p} \mathbf{Z_{2}} \mathbf{Z}_{2}^{\top} \]

我们发现标准化的G阵就是GCTA计算用的G阵，而中心化的G阵相比于 VanRaden 构建的G阵（见下式）常数项不同， \[ \mathbf{G}=\mathbf{Z Z}^{\prime} / \sum 2 p_{i}\left(1-p_{i}\right) \] 作者解释了一下何时使用这两种矩阵：如果MAF比较低的SNP容易有更大的效应，那么我们倾向于使用标准化的G阵，反之则倾向于中心化的G阵（理由呢？）。但是根据作者的经验，中心化的G阵对于校正群体结构效果更好。原文如下：

结果文件

• prefix.log.txt 日志文件
• prefix.cXX.txt 或者 prefix.sXX.txt 是n*n 的矩阵文件。

一般线性模型

./gemma -bfile [prefix] -lm [num] -o [prefix]

-bfile : plink 二进制文件的前缀

-lm [] 表示用哪个方法进行检验（which frequentist test to use），就是得到P值的不同检验方式。

• -lm 1 表示 wald test
• -lm 2 表示 likelihood ratio test
• -lm 3 表示 score test
• -lm 4 表示 执行上述三种检验

-o 结果文件名称

对于二分类数据，你可以将数据标记为 0 和 1 ，GEMMA 会将二分类数据视为数量性状进行分析（震惊，居然不是用逻辑回归）。这种做法的合理性，可以将这种做法视为广义线性模型的一阶泰勒展开 (?)。

结果文件

结果文件会在当前目录的 output 文件夹中。

prefix.log.txt

日志文件

prefix.assoc.txt

关联结果，这里的 beta 估计值就是这个位点的斜率，最后一列 p_wald 表示 Wald test 计算得到的 P 值。

单性状混合模型GWAS

./gemma -bfile [prefix] -k [filename] -lmm [num] -o [prefix]

-bfile 二进制plink文件前缀

-k 关系矩阵名称

-lm [] 表示用哪个方法进行检验（which frequentist test to use）

• -lm 1 表示 wald test
• -lm 2 表示 likelihood ratio test
• -lm 3 表示 score test
• -lm 4 表示 执行上述三种检验

如果想检测G×E效应，可以增加 -gxe [filename]，gxe 文件中包括一列环境变量。 此时，对于每一个SNP，GEMMA 除了拟合这个 SNP 和环境变量的效应外，还会检验这二者的互作效应。

详细信息

这个算法会通过极大似然估计 (MLE) 或限制极大似然估计 (REML) 方法来估计 \(\lambda\) (应该是加性方差与残差的比值) ，其范围为 \(1 \times 10^{-5}\) (纯粹是环境效应) 到 \(1 \times 10^{5}\) (纯粹是加性效应) 。

对于二分类性状，你可以标记为 0 和 1，GEMMA 会将其视为数量性状进行计算（实锤了，没用逻辑回归）。

结果文件

• prefix.log.txt ：日志文件，其中还包含了 PVE (遗传力) 估计值及其标准误
• prefix.assoc.txt ： GWAS 结果文件，示例如下图。这里新增了倒数第二列，l_remle 为 \(\lambda\) 的限制极大似然估计值（REML）。

多性状混合模型GWAS

./gemma -bfile [prefix] -k [filename] -lmm [num] -n [num1] [num2] [num3] -o [prefix]

这和单性状完全相同，除了 -n 选项，用于指定表型，以空格分隔。（例如 -n 1 2 就是第六列和第七列）。

详细信息

如果表型中有少量缺失值，可以在进行GWAS分析先进行表型值填充（这里没说原理，但是我估计是用 GS 进行预测，我个人不建议这么做，建议剔除缺失表型）。

./gemma -bfile [prefix] -k [filename] -predict -n [num1] [num2] [num3] -o [prefix]

输出文件

• prefix.log.txt ：日志文件，其中还包含了 PVE (遗传力) 和遗传相关。
• prefix.assoc.txt ： GWAS 结果文件，和单性状类似，最后一列是 p 值，p 值前面是几个表型的 beta 估计值，以及这些估计值的方差矩阵 (啥叫 beta 估计值的方差矩阵?)。