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

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=VCFv4.2
##ALT=<ID=NON_REF,Description=“Represents any possible alternative allele at this location”>
##FILTER=<ID=LowQual,Description=“Low quality”>
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 081
NC_000962.3 1977 . A G 18395.03 . AC=2;AF=1.00;AN=2;DP=470;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=25.36;SOR=0.791 GT:AD:DP:GQ:PL 1/1:0,470:470:99:18409,1410,0
NC_000962.3 4013 . T C 16093.03 . AC=2;AF=1.00;AN=2;DP=400;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=28.73;SOR=0.776 GT:AD:DP:GQ:PL 1/1:0,400:400:99:16107,1202,0

2.2.2*.bed文件

包含基因型信息，以二进制格式存储（因为格式问题难以展示示例）。

2.2.3*.fam文件

包含表型信息，由六列组成的文本文件，从左到右分别为：

示例：
081 081 0 0 0 1
0 278 0 0 0 1
0 279 0 0 0 1
0 283 0 0 0 0
0 289 0 0 0 1
0 291 0 0 0 1
0 299 0 0 0 1

2.2.4*.bim文件

存储基因型和位置等信息，一个包括6列的文件。

示例：
0 . 0 21 CAGT AGGC
0 . 0 29 C A
0 . 0 30 G C
0 . 0 42 T C
0 . 0 43 G A

2.2.5id.txt文件

为了批量处理文件，将样本id分别作为一行，写入文件。
示例（id.txt）：
081
082

2.3实际操作_安装软件

2.3实际操作_流程

2.4.1实际操作_流程_基础文件夹配置

创建每一步需要的文件夹，并把数据文件和参考基因组文件分别拷贝进00_data文件夹，结构如下图
图5.文件设置

2.4.2实际操作_流程_质控和比对

#使用fastp进行质控
for name in `cat 02_id/id.txt` ;do /data/home/lizhiyuan/biotools/fastp/fastp\
-i 00_data/${name}_1.fastq.gz -I 00_data/${name}_2.fastq.gz  \
-o 01_QC/${name}_1.fastq.gz -O 01_QC/${name}_2.fastq.gz ; done
#为参考基因组添加注释
gatk CreateSequenceDictionary -R 09_reference/GCF_000195955.2_ASM19595v2_genomic.fna -O \ 09_reference/GCF_000195955.2_ASM19595v2_genomic.dict && \
bwa index 09_reference/GCF_000195955.2_ASM19595v2_genomic.fna && \
samtools faidx 09_reference/GCF_000195955.2_ASM19595v2_genomic.fna#使用bwa进行比对，而后通过管道“|”将比对后的文件传递给samtools，使用samtools得到样本测序深度
for name in `cat 02_id/id.txt` ;do rm -rf 04_bam/$name && \
mkdir 04_bam/$name && bwa mem  -R '@RG\tID:${name}\tSM:${name}\tLB:lib' \
09_reference/GCF_000195955.2_ASM19595v2_genomic.fna 01_QC/${name}_1.fastq.gz \
01_QC/${name}_2.fastq.gz | samtools view -Sb -o 04_bam/${name}/${name}.bam && \
samtools sort 04_bam/${name}/${name}.bam  >  04_bam/${name}/${name}_sorted.bam && \
samtools index 04_bam/$name/${name}_sorted.bam && \
samtools depth  -aa  04_bam/${name}/${name}_sorted.bam > 04_bam/${name}/${name}.depth;done

你可以在这里找到更多关于bwa软件使用的信息5.

2.4.3实际操作_流程_合并并形成*.vcf文件[6]

#将*.bam格式的文件转换成*.g.vcf格式的中间文件
for name in `cat 02_id/id.txt` ;do gatk HaplotypeCaller \-R 09_reference/GCF_000195955.2_ASM19595v2_genomic.fna \--emit-ref-confidence GVCF \-I 04_bam/${name}/${name}_sorted.bam \-O 05_gvcf/${name}.g.vcf;done#将*.g.vcf文件合并
gatk CombineGVCFs -R 09_reference/GCF_000195955.2_ASM19595v2_genomic.fna \
--variant 05_gvcf/081.g.vcf \
--variant 05_gvcf/082.g.vcf \
-O 07_vcf/merge.g.vcf#将合并后的merge.g.vcf文件转化成*.vcf文件
for name in `cat 02_id/id.txt` ;do gatk GenotypeGVCFs \-R 09_reference/GCF_000195955.2_ASM19595v2_genomic.fna \-V 07_vcf/merge.g.vcf \-O 07_vcf/merge.vcf ;done

2.4.4实际操作_流程_转换格式

#将vcf格式转换为plink格式
cd 06_plink && vcftools --vcf ../07_vcf/merge.vcf --plink --out merge_plink && plink --file merge_plink --make-bed --out merge_plink

2.4.5实际操作_流程_单个基因和表型的关联

#将表型信息输入到fam文件（GEMMA开发者建议使用0作为对照组，使用1作为实验组），使用GEMMA计算表示遗传距离远近的矩阵，用于去除连锁不平衡带来的假阳性
#创建文件用于储存矩阵
mkdir 08_results/00_matrix
#-bfile表示plink文件的前缀;-gk表示计算的矩阵类型，1表示把数据归一化，2表示计算数据的标准差;-o表示产出文件的前缀是什么
gemma -bfile 06_plink/merge_plink -gk 1 -o merge_plink_matrix
#将产生的矩阵移动到
mv output/* 08_results/00_matrix && rm -rf output
#分析基因型和表型的关联,-bfile plink二进制文件前缀; -k 指定遗传距离矩阵的位置; -lmm 选择计算显著性的方法，1表示使用wald检验计算显著性，2表示使用likelihood ratio test计算显著性，3表示使用score test计算显著性
gemma -bfile 06_plink/merge_plink -k 08_results/00_matrix/merge_plink_matrix.cXX.txt -lmm 1 -o merge_plink_assoc

2.4.6实际操作_流程_基因和多个表型的关联

#-n表示在表型文件中哪几个表型用于分析
gemma -bfile 06_plink/merge_plink -k 08_results/00_matrix/merge_plink_matrix.cXX.txt -lmm 1 -n 6  -o

3.结果分析

3.1结果分析_文件

产生的结果文件merge_plink_assoc.assoc.txt中，从左到右分别为：

示例：
chr ps n_mis n_obs allele1 allele0 af beta se p_wald
0 4187485 5 4271 TG CG 0 -4.88E-01 7.64E-02 1.86E-10
0 3829770 12 4264 T C 0 -4.89E-01 7.64E-02 1.82E-10

3.2结果分析_整合多个结果文件

可以使用python中的pandas[7]包高效整合不同*assoc.txt文件

# -*- coding: utf-8 -*-
import pandas as pd
import glob
import os
os.chdir(r"D:/linux_lizhiyuan/MTB_ann_nohead_noPEPPE_assoc/outputs_matrix_standard_assoc")
i = 0
#设置i为0，即False
for filename in glob.glob("*.assoc.txt"):
#批量打开文件名尾为.assoc.txt的文件print(filename)feature = pd.read_csv(filename, delimiter="\t", header=0)feature['file'] = filename#添加一列，储存文件来源feature ['Software'] = "GEMMA"#添加一列，储存使用的软件feature=pd.DataFrame(feature)#将数据设定为数据框格式if i == False:feature[1:].to_csv('assoc_gemma_merge.csv', mode='+a', index=None)#将第一个文件写入新表格时保留表头if i == True:feature[1:].to_csv('assoc_gemma_merge.csv', mode='+a', header=None, index=None)# 将第二个或更多文件写入新表格时省略表头i = True#给i赋值，进而使i成为正值（Ture）

合成的文件格式如下：

3.3结果分析_绘图

根据得到的P值，可以使用R包qqman [8] 和ggplot2[9]分别绘制QQ图曼哈顿图，QQ图中点和线越靠近结果越好，曼哈顿图中位点的数据点位置越高越显著，有关二图图注含义的解释可以查看参考资料[10]。

library(qqman)
library(ggplot2)
setwd("08_results/01_assoc")gwasResults1 <- read.delim('merge_plink_assoc.assoc.txt', sep = '\t')
#读入文件，以制表符分隔gwasResults1$SNP1 <- seq(1, nrow(gwasResults1), 1)
#为了绘制X轴，设定gwasResults的新列名称为SNP1，内容为从1到行数目的数列
gwasResults1$ps=as.numeric(levels(gwasResults1$ps))subdata <- subset(gwasResults1,P < 0.00001)
#可以根据P值等条件筛选数据，进而只对部分数据注释文本#ggplot2 绘制曼哈顿图
p_ggplot2 <- ggplot(gwasResults1, aes(ps, -log(P, 10)),color = ps) +#以SNP1为x轴，以-log（P，10）为Y轴，以作为绘图条件作图theme(panel.grid = element_blank(), axis.line = element_line(colour = 'black'), panel.background = element_rect(fill = 'transparent')) +#绘制横轴和纵轴geom_point(size=0.9) +#绘制点，并且不显示各个点对应的染色体的图注labs(x = 'ps', y = expression(''~-log[10]~'(P)')) +#绘制横纵轴的坐标#取每一个染色体的样本编号的中位数作为分界，以染色体的值作为标签geom_hline(yintercept = c(5, 10), color = c('blue', 'red'), size = 0.5)  +#在-log(P, 10)为5和10的地方分别绘制绿线和红线geom_text(data=gwasResults1,aes(label=ps,color="blue"),nudge_y = 0, nudge_x = 80,check_overlap = TRUE, show.legend = FALSE,size = 1.5)
#这里如果data选择subdata，即可只对部分数据注释文本
ggsave("merge_plink_assoc.assoc.txt_man.png", p_ggplot2, width = 10, height = 4)
#保存图片（图4）#qqman绘制qq图
qq <- qq(gwasResults1$P, main = "Q-Q plot of GWAS p-values", xlim = c(0, 7), ylim = c(0, 12), pch = 18, col = "blue4", cex = 1.5, las = 1)
ggsave(' merge_plink_assoc.assoc.txt_qqplot.png', qq, width = 10, height = 4)
#绘制并保存QQ图

图6.曼哈顿图（横轴为碱基位置，假设基因组长度在4500左右）
图7.qq图。

4.参考资料

[1]:X Zhou,M Stephens. Efficient multivariate linear mixed model algorithms for genome-wide association studies[J]. Nat Methods, 2014, 11(4): 407-9.
[2]:JH Sul, LS Martin,E Eskin. Population structure in genetic studies: Confounding factors and mixed models[J]. PLoS Genet, 2018, 14(12): e1007309.
[3]:X Zhou,M Stephens. Genome-wide efficient mixed-model analysis for association studies[J]. Nature Genetics, 44(7): 821-4.
[4]:https://www.jianshu.com/p/957efb50108f.
[5]:https://github.com/lh3/bwa
[6]:X Zhou. GEMMA User Manual[J] May 18, 2016.
[7]:https://pandas.pydata.org/
[8]:https://cran.r-project.org/web/packages/qqman/vignettes/qqman.html.
[9]:http://ggplot2.tidyverse.org, https://github.com/tidyverse/ggplot2.
[10]:https://www.jianshu.com/p/987859ae503c.

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

1. 对基因型和表型做t.test

   t.test参考部分 参考链接: https://zhuanlan.zhihu.com/p/126351774 https://zhuanlan.zhihu.com/p/123907459 因为只有数 ...

2. Python --深入浅出Apriori关联分析算法（二） Apriori关联规则实战

   上一篇我们讲了关联分析的几个概念,支持度,置信度,提升度.以及如何利用Apriori算法高效地根据物品的支持度找出所有物品的频繁项集. Python --深入浅出Apriori关联分析算法(一) 这次 ...

3. FP-Growth关联分析算法在网络监控领域的应用

   关联分析算法在网络监控领域的应用: 在现今网络规模大,涉及专业多,告警总数大的现在,迫切需要提高对海量告警的分析能力,实现对告警数据的挖掘,提高对有价值告警的提取,简化监控人员的工作,提高排障效率.常 ...

4. pythonapriori算法特点_Python --深入浅出Apriori关联分析算法（一）

   在美国有这样一家奇怪的超市, 它将啤酒与尿布这样两个奇怪的东西放在一起进行销售,并且最终让啤酒与尿布这两个看起来没有关联的东西的销量双双增加 .这家超市的名字叫做沃尔玛. 你会不会觉得有些不可思议?虽 ...

5. 一个实例带你搞懂Apriori关联分析算法

   关联分析 Apriori算法 优点:易编码实现. 缺点:在大数据集上可能较慢. 适用数据类型:数值型或者标称型数据. 关联分析是一种在大规模数据集中寻找有趣关系的任务.这些关系可以有两种形式:频繁项集 ...

6. python 关联分析算法的包_Python 极简关联分析（购物篮分析）

   关联分析,也称购物篮分析,本文目的: 基于订单表,用最少的python代码完成数据整合及关联分析 文中所用数据下载地址: 使用Python Anaconda集成数据分析环境,下载mlxtend机器学习 ...

7. Apriori关联分析算法 -尿布与啤酒的故事

   ✴关联分析概述 选择物品间的关联规则也就是要寻找物品之间的潜在关系.要寻找这种关系,有两步,以超市为例 找出频繁一起出现的物品集的集合,我们称之为频繁项集.比如一个超市的频繁项集可能有{{啤酒,尿布} ...

8. python --深入浅出Apriori关联分析算法Apriori关联...

    一.基础知识 上次我们介绍了几个关联分析的概念,支持度,置信度,提升度.这次我们重点回顾一下置信度和提升度: 置信度(Confidence):置信度是指如果购买物品A,有较大可能购买物品B.计算方式 ...

9. python关联规则apriori算法_Python --深入浅出Apriori关联分析算法（二） Apriori关联规则实战...

   上一篇我们讲了关联分析的几个概念,支持度,置信度,提升度.以及如何利用Apriori算法高效地根据物品的支持度找出所有物品的频繁项集. 这次呢,我们会在上次的基础上,讲讲如何分析物品的关联规则得出关联 ...

最新文章

1. GitHub日收12000星，微软新命令行工具引爆程序员圈！
2. 【 Vivado 】工程模式下运用Tcl脚本示范
3. C++标准库与STL简介
4. HDU 5828 Rikka with Sequence (线段树+剪枝优化)
5. 329. 矩阵中的最长递增路径
6. bootstrapV5+（资源篇）
7. Cesium场景导出为图片
8. 获得执行计划方法三-sql_trace
9. 计算机网络之JSONP跨域
10. 8Manage轻松解决采购过程“脏乱差”问题
11. 检查Mysql引擎的方法
12. Android基于代理的插件化思路分析
13. 查找所引用的文献在某种期刊下的引用格式（引用风格）
14. DiskGenius 强行拆分黑苹果HFS硬盘分区以给Windows扩容
15. 读书郎上市背后隐忧：业绩下滑明显，市场地位较靠后，竞争力存疑
16. Unity 3D作业七：人物模型
17. RationalDMIS 2020 网络报表/网络编程连接设置
18. 免费下论文及查重投稿的10来个方法
19. word详细使用方法（①）
20. 微信公众号小程序如何做流媒体视频直播？

热门文章

1. uni-app移动端-H5-微信小程序在线预览pdf，图片，视频
2. 部署KVM虚拟化平台
3. 60 个程序员才懂的梗！太形象了！(笑死了！！！！！！！！！)
4. java 多线程并行_Fork and Join: Java Can Excel at Painless Parallel Programming Too! | Oracle 中国...
5. 考研冲刺阶段注意事项？看完一定能考上
6. 你真的了解HR问你的问题么？
7. 微信小程序收集formId
8. 图书分类怎么写用php,PHP开发简单图书后台管理系统实现图书统计
9. C# 模仿360安全卫士玻璃按钮 修正版（源码）
10. 如何在 Nervos CKB 上开发智能合约