一文学会基因型和表型关联分析算法GEMMA

备注：本文同时发布在公众号“基因的生物信息学分析”，链接在这里

1.介绍

1.1 介绍_简介

GEMMA全称为：Genome-wide Efficient Mixed Model Association algorithm，即基于全基因组混合模型关联算法的工具[1]。GEMMA是一款基于混合线性模型的GWAS分析软件。可精确和快速地进行单SNP的GWAS、多SNP的GWAS和多性状的GWAS分析。

1.2 介绍_优点

1.2.1介绍_优点_排除了连锁不平衡的干扰3级标题

人类和细菌群落中本身的SNP差异造成了突变的连锁不平衡（即基因组本身序列有差异），造成了组内差异大于组间差异，进而导致通过GWAS分析找到的显著表达的基因只是由不同群体基因组的差异引起的，产生了大量的假阳性(图1)[2]。
图1.群里结构或遗传相关性带来的假阳性
而GEMMA以混合线性模型关联表型和基因型，通过计算群体中个体基因型的相似性矩阵，在一定程度上排除了群体基因组连锁不平衡的干扰，进而去除了这些假阳性（图2.）。GEMMA的算法计算式子如下：

表型=平均数+基因型效应+群体结构+环境

Y = W * α + X * β + u + ɛ

其中u~ MVN_n(0,λτ^-1K)， ɛ ~ MVN_n(0,λτ^-1In)
即
零假设H0：每一个SNP的β=0（即基因型与表型无关）， 备择假设H1：存在β≠0，存在某个SNP的β≠0（即基因型与表型相关）
图2.计算相似性矩阵并结合多元线性回归模型

1.2.2介绍_优点_速度快

使用固定的处理器配置对数据进行分析在比较中，GEMMA在精确算法中速度最快，在近似算法中也仅仅慢于GRAMMAR；GEMMA准确度和EMMA更接近(图3)[3]。但图3中只有一部分算法，新的算法如phylogenetic convergence test还在不断的被开发出来。
图3.方法学比较

处理器：
n Intel Xeon L5420 2.50 GHz CPU
数据：
HDL-C:m=99, n=681 and p=1,885,197
Crohn‘s disease：m=n=4686 and p=442,001

2.实际操作

2.1实际操作_分析流程概述

图4.流程

2.2实际操作_输入文件格式

GEMMA支持BIMBAM格式和plink二进制格式，主要包括以.bed、.fam、.bim。这3种文件需要由标准格式的.vcf文件转化而来。限于篇幅，只选取部分文件讲解格式。

2.2.1*.vcf文件

vcf（Variant Call Format）文件是存储变异位点的标准格式，可以用来表示单核苷酸多态性、插入缺失、结构变异、拷贝数变异等。因为vcf格式很复杂，本文中未提及的请查看https://samtools.github.io/hts-specs/或http://www.cog-genomics.org/plink/2.0/input#vcf。
前几行以“#”开头的行的含义：

简称	含义
##fileformat	VCF格式版本号
##FILTER	显示这个文件已经进行了过滤
##contig	参考基因组contig信息
##INFO	INFO列中各简写的含义。ID、Number、Type、Description主要有几个tag标记：AD、DP、GQ、GT、PL

不以“#”开头的信息列的含义：

简称	含义
CHROM	表示变异位点位于哪个染色体
POS	变异位点相对于参考基因组所在的位置，若为删除或位移第一个碱基所在的位置。
ID	变异位点名称（对应dbSNP数据库中的ID；若没有，则默认用）
ALT	该位点突变后的碱基类型类型，若有多个，则用逗号分隔。
REF	该位点参考基因组的碱基类型。
QUAL	可以理解为所call出来的变异位点的质量值Q。Q=-10lgP，P表示这个位点发生错误的概率。因此，当Q=20时，错误率P=0.01。
FILTER	针对质量值等变量过滤之后，在FILTER一栏都会留下过滤记录，通过过滤的位点FILTER一栏显示“PASS”；没有通过过滤的位点的FILTER一栏为其它信息。若FILTER一栏为“.”，说明没有进行过滤。
INFO	有关该位点的额外信息
FORMAT	变异位点格式，其中字母简写的含义在文件中以“#”开头的行中
SMAPLE	使用的样本名称，由bam文件中@RG的SM标签决定。对应的数据必定有GT（genotype，样本基因型）信息。
GT:AD:DP	GT信息中，两个数字之间以斜线分隔则表示二倍体样本位于两条染色体上的基因型。0表示该位点的与参考基因