本篇笔记主要记录如何在linux命令行中用TASSEL5和PLINK2软件做GWAS分析，以及如何可视化GWAS结果（在R中绘制曼哈顿图和QQ图）。关于GWAS是做什么的，以及基本概念介绍可以看上一篇笔记介绍：群体基因组学——概述和图表详解 - 我的小破站 (shelven.com)
PLINK和TASSEL都是做GWAS分析用的经典软件，两个软件都有linux版本和windows版本，我这里是在linux集群环境中跑的分析流程，所有软件下载linux版本，两款软件可以在下方官网中安装：
PLINK 2.0 (cog-genomics.org)、Tassel
TASSEL官网有tassel pipeline命令行的帮助文档，一些常用的参数可以参考Tassel5PipelineCLI.pdf (bytebucket.org)
TASSEL中文文档可以参考本站Tassel5中文操作手册.pdf
1. 数据格式
plink的数据格式有两套，每套各自的前缀名称相同，一套后缀为.bed、.bin和.fam，另一套后缀为.map和.ped，后面的分析以第一套为例，读取速度比较快。
那么这三个是怎么来的呢？我们做群体重测序后，通过GATK等软件做变异检测最终会生成一个VCF文件，通过软件plink可以将VCF文件转化成上面的三个文件，分别展示一下三个文件的内容和结构：
- fam文件（家族信息文件）有6列，分别为家系编号、个体编号、父系编号（0表示信息缺失）、母系编号（0表示信息缺失）、性别编号（1表示男，2表示女，0表示未知性别）、表型值（1表示对照，2表示病例，0/-9或者其他非数字表示信息缺失）
- bim文件（个体信息文件）有6列，分别为染色体编号、SNP标识、以摩根或厘摩为单位的位置（可以用0）、碱基对坐标、Allele1和Allele2（通常是主效等位基因）
- bed文件为二进制文件
转化成这三个文件主要是为了下一步更快过滤数据，因为bed是二进制文件，计算机处理起来更快。
以上数据来源是贺师兄做的棉花重测序样本的部分变异信息文件。
2. 基因型数据清洗（使用PLINK）
使用plink软件过滤掉最小等位基因频率（Minor Allele Frequency，MAF）小于0.05的变异位点，在关联分析中MAF值太小会造成假阳性。比如说，一个基因座上有两个等位基因A和B，A在群体中的频率为0.6，B在群体中的频率为0.4，那么MAF值就是0.4，可以想到一个基因座上MAF值越小，基因座上的等位基因越单一，MAF太低的位点贡献的信息很少，不与表型做关联分析。
1 | # 过滤MAF小于0.05的SNP，重新生成.bed、.bin和.fam文件 |
可以通过查看过滤前后的bim文件，统计过滤了多少位点（plink软件在屏幕上的标准输出也会显示）：
1 | # plink输出结果 |
需要注意下我这里使用的表型数据都是处理好的，所以只对基因型数据做过滤就行，如果自己有表型数据还要做一下表型数据清洗，比如删除异常值等等。
生成的GWAS_vcf.vcf用于下一步关联分析。
3. 关联分析（使用TASSEL5）
PLINK软件也可以做关联分析，网上也能找到一把教程，但是PLINK只能做GLM模型的关联分析，现在文献中用的比较多的软件是TASSEL5，接下来演示用命令行的Tassel Pipeline做关联分析。
3.1 计算亲缘关系矩阵
1 | # -Xmx10G 控制最大堆（heap size）的大小为10G |
生成的216个样本的亲缘关系矩阵kinship.txt如下：
3.2 主成分分析
1 | # -fork<id> 用于区分不同的pipeline，代表一个pipeline的起始，后面可以是数字或者字符（非空格） |
默认情况下会生成的PCA结果主要展示前5个主成分，且这一步会生成三个txt文件，PCA结果展示在第一个文件pca1.txt：
3.3 可视化亲缘关系矩阵和PCA结果（群体结构分析）
在拿到这两个矩阵之后就可以通过R可视化结果，首先是亲缘关系矩阵热图：
1 | library(data.table) |
用颜色深浅不同表示每个样本之间的亲缘关系远近，颜色越深亲缘关系越近。
PCA结果作图：
1 | library(data.table) |
这里的群体结构差异有但不能很好分类。如果做群体结构分析的时候可以分成明显的几个不同的部分，那就要考虑后续做GWAS分析过程中要考虑群体结构的分层问题了。这里我们只是拿部分数据跑个简单的流程，继续往下。
3.4 基于GLM模型进行GWAS分析
现在分别有了如下的基因型数据（GWAS_vcf.vcf）、表型数据（phenotype.tassel.txt），就可以以PCA结果作为协变量，基于GLM模型进行GWAS分析：
1 | # -excludeLastTrait 移除表型数据的最后一列（不太明白什么意思？官网说可用于删除的最后一列用于MLM或GLM的群体结构） |
主要结果为glm1.txt，结果如下：
我们主要用到的数据是：
- 第一列：表型
- 第二列：SNP名称
- 第三列：染色体编号
- 第四列：处于染色体的位置
- 第六列：p值
3.5 基于MLM模型进行GWAS分析
在tassel中运行MLM模型和GLM模型相似，MLM模型需要亲缘关系矩阵来定义个体之间的关系，以PCA和kinship作为协变量，基于MLM模型进行GWAS分析：
1 | # -mlm -mlmVarCompEst P3D -mlmCompressionLevel Optimum 使用P3D和最优水平压缩 |
主要结果为mlm6.txt（前5个是5个表型的分析数据），结果如下：
我们主要用到的数据是：
- 第一列：表型
- 第二列：SNP名称
- 第三列：染色体编号
- 第四列：处于染色体的位置
- 第七列：p值
3.6 计算核酸多态性值
1 | vcftools --vcf GWAS_vcf.vcf --site-pi --out pi |
用vcftools就可以根据vcf文件统计平均每个位点的π值。
4. GWAS结果可视化
1 | # qqman包绘制曼哈顿图和qq图 |
生成各个表型GWAS分析的曼哈顿图和qq图：
随便查看一个MLM模型的GWAS曼哈顿图和qq图结果：
