2020.11.9【WGS/GWAS】丨全基因组分析（关联分析）全流程(下)

2021.07.30 已更新数据分析部分内容：点我查看

一. 摘要

经过为期半个月的东拼西凑 研发测试，作者终于整理出了一个从VCF开始的GWAS后期分析流程。当然要感谢很多大佬提供的代码 帮助，在文章中也附上参考链接。对GWAS还不够熟悉的朋友，可以看一下我之前整理的一份PPT学习笔记《遗传进化与GWAS研究》。

二. SNP 检测与注释

使用软件：snpEff
安装方式：
miniconda3: conda install snpeff
SnpEff and SnpSift (pcingola.github.io)
使用流程：
下载参考基因组及注释文件（注意vcf染色体为数字，使用RefSeq或者GenBank编号需要统一染色体序列编号）

java -jar /home/yangxin/miniconda3/pkgs/snpeff-5.0-0/share/snpeff-5.0-0/snpEff.jar build -c /home/yangxin/miniconda3/pkgs/snpeff-5.0-0/share/snpeff-5.0-0/snpEff.config -gtf22 -v GCF_002163495.1_Omyk_1.0_genomic
#参数说明
#java -jar: Java环境下运行程序
#-c snpEff.config配置文件路径
#-gff3 设置输入基因组注释信息是gff3格式，如果是gtf文件请改为-gtf22
#-v 设置在程序运行过程中输出的日志信息

最后的-v参数 设置输入的基因组版本信息，和snpEff.config配置文件中添加的信息一致

java -jar /home/yangxin/miniconda3/pkgs/snpeff-5.0-0/share/snpeff-5.0-0/snpEff.jar build -c /home/yangxin/miniconda3/pkgs/snpeff-5.0-0/share/snpeff-5.0-0/snpEff.config -gff3 -v GCF_002163495.1_Omyk_1.0_genomic #搭建基因组数据库

处理染色体编号（可选）
由于我下载的基因组使用的RefSeq序列编号，于是，我需要对染色体号进行转换。

awk '{$1="";print $0}' OFS="\t" AxiomGT1.calls.vcf > col_2.vcf 
如果你的染色体与数据库不一致，那么可以手动生成染色体序列文件col_chr2.vcf
paste -d '' col_chr2.vcf col_2.vcf >AxiomGT1.vcf #合并序列
sed -n l AxiomGT1.vcf #检查vcf文件格式是否一致，注意分隔符为“\t”

注意：cat、head、less看不出来分隔符号的问题

java -jar ~/miniconda3/pkgs/snpeff-5.0-0/share/snpeff-5.0-0/snpEff.jar eff -c ~/miniconda3/pkgs/snpeff-5.0-0/share/snpeff-5.0-0/snpEff.config GCF_002163495.1_Omyk_1.0_genomic AxiomGT1.vcf > AxiomGT1.snp.eff.vcf -csvStats AxiomGT1_positive.csv -stats AxiomGT1_positive.html #对snp进行注释
awk -F "|" '{print $1,$2,$4,$5,$8,$9}' OFS="\t" AxiomGT1.snp.eff.vcf > SNV.anno.xls #生成注释文件
#参数-F "|"用于处理表头信息，如果表头正常可以不用这个参数

结果展示：注释文件表格
流程参考：SnpEff使用方法

三. 群体SNP 过滤

使用软件：vcftools
安装方式：
miniconda3: conda install vcftools
使用流程：
基本过滤参数 Calling rate < 95%，次要等位基因计数 < 3。

vcftools --vcf AxiomGT1.calls.vcf --max-missing 0.95 --mac 3 --recode --recode-INFO-all --out AxiomGT1.calls.filter.vcf
#--max-missing 0.95，即过滤calling rate低于95％的基因型
#--mac 3过滤次要等位基因计数小于3的SNP

vcftools --vcf AxiomGT1.calls.filter.recode.vcf --missing-indv
#检测样本数据丢失率,丢失率过高需要调整参数。
结果展示
过滤前：
过滤后：

四. 群体分层分析

群体分层是指群体内存在亚群的现象，亚群内部个体间的相互关系大于整个群体内部个体间的 平均亲缘关系。由于不同的亚群间，某些稳定的等位基因频率不同，当将两个亚群混合进行关联分 析时， 就会导致假阳性结果的产生。所以，进行关联分析之前，一定要先进行群体分层分析。

群体分层分析有：群体进化树分析和主成分分析，两者的结果可以进行相互验证。

4.1 群体进化树分析

绘制进化树的软件比较多，但是测试几个之后还是觉得tassel可视化界面绘图最方便，下面给大家介绍一下
使用软件：Tassel
软件安装：miniconda3: conda install tassel5.0
官网下载：Tassel (bitbucket.io)
使用流程（可视化画图）：
step1：打开Tassel
step2：导入并打开vcf文件

step3：创建进化树文档

step4：绘制进化树

step5：调整图像

功能简介：
左侧可以针对输入文件的参数进行增删；
左下角X,Y±用来调整进化树的长宽；
左上角File可以输出各种格式的文件；
在工具栏Type可以选择进化树类型。（最后一个圈图比较好看）
结果展示

注意：因为VCF文件在Tassel进行关联分析的时候有最大字符限制，只截取前面10个字符，但是绘制进化树不会处理字符。所以担心绘制进化树观看不方便的话可以提前处理。

Tassel使用参考：https://bitbucket.org/tasseladmin/tassel-5-source/wiki/UserManual
使用手册中文版(2014)：https://wenku.baidu.com/view/4c0cabbe9b6648d7c0c746b2.html

4.2 群体主成分分析

使用软件：Plink
软件安装：miniconda3: conda install Plink
使用流程：

vcftools --vcf AxiomGT1.calls.vcf --plink --out GT1
#生成plink需要的ped/map格式
plink --file GT1 --make-bed --out GT1 --chr-set 30
#转换为bed格式,非人类染色体，使用--chr-set参数设置染色体数目
Plink --threads 8 --bfile snp --pca 10 --out pca
#--threads 线程数 --bfile 输入.bed文件 --pca 主成分的成分数 --out输出的文件名

###########R-PCA###########

mydat<-read.table("pca.eigenvec",as.is = T,header = T,stringsAsFactors = F)
png("PCA.png",width = 7000,height = 7000,pointsize = 160)

结果展示

Tassel，GCTA 或者 EIGENSOFT也可以进行分析，可以多尝试一下。

4.3 群体遗传结构分析

这一部分我在之前有专门的文章进行使用报错介绍，可以直接跳转参考：Admixture群体遗传结构分析
使用软件：admixture
软件安装：miniconda3: conda install admixture
使用流程：
前面与主成分分析相同，需要使用plink将VCF文件转换成bed/bim/fam几个文件

vcftools --vcf AxiomGT1.calls.vcf --plink --out GT1
#生成plink需要的ped/map格式
plink --file GT1 --make-bed --out GT1 --chr-set 30
#转换为bed格式,非人类染色体，使用--chr-set参数设置染色体数目

Admixture工具处理

for K in 2 3 4 5 6 7 8 9 10; \
do admixture --cv GT1_test.bed $K | tee log${K}.out; done
#2 3 4 5 6 7 8 9 10分成的群体结构数 hapmap3.bed 输入文件
grep -h CV log*.out

#查看最佳K值（最小值）， 输出最佳K值文件或展示所有K值结果
#R绘图

for (k in 2:10){
file_name = paste("GT1_test.",k,".Q",sep="")
result_name = paste(file_name,".PDF",sep="")
pdf(result_name,width=15,height=3)
tbl <- read.table(file_name,header = T)
barplot(t(as.matrix(tbl)), col=rainbow(k),xlab="Individual", ylab="Ancestry", border=NA)
dev.off()
}

结果展示：
以k=3，k=8为例进行展示

五. 连锁不平衡分析

2018 年 10 月 15 日，Bioinformatics 杂志上在线发表了 PopLDdecay 软件，用于分析连锁不平衡的衰减（Linkage disequilibrium (LD) decay）。PopLDdecay 可以直接读取 VCF 文件，相比于计算 LD 常用的 PLINK 软件，这一特点简化了格式转换的繁琐。

使用软件：PopLDdecay

源码下载：PopLDdecay-3.41

使用流程：
如果想要分表型进行绘制，我们需要先建立一个list，分表型进行计算，list格式为一列该表型的样品名称

PopLDdecay -InVCF AxiomGT1.calls.vcf -OutStat RESISTANT -SubPop RESISTANT.list #分型计算
gunzip RESiSTANT.stat.gz #其他表型操作相同

根据不同表型建立分组list，单行格式为stat路径，分隔符，stat文件前缀

perl ~/soft/PopLDdecay-3.41/bin/Plot_MultiPop.pl -inList file.list -output Fig

结果展示

参考链接：PopLDdecay: 基于 VCF 文件的快速、高效计算连锁不平衡的工具

六. 选择消除分析

群体进化选择消除分析简单的说就是基因组某区域由于受到了选择而消除多态性。选择消除分析是正向选择在物种基因组上留下的印迹。与野生祖先相比，栽培或驯化的物种发生选择消除的区域遗传多样性显著降低，这是驯化区域的典型特征。

使用软件：vcftools
软件安装：略

6.1 Fst 分析

群体的固定系数F反映了群体等位基因杂合性水平。固定系数F是F统计量（Fst）的一个特例。Fst分析表示群体的分化程度，值越大，群体分化程度越高，受选择程度越高。

##对每一个SNP变异位点进行计算
vcftools --vcf AxiomGT1.calls.vcf --weir-fst-pop 1_population.txt --weir-fst-pop 2_population.txt  --out Fst_AximoGT1

##按照步距来计算
vcftools --vcf test.vcf --weir-fst-pop 1_population.txt --weir-fst-pop 2_population.txt  --out p_1_2_bin --fst-window-size 500000 --fst-window-step 50000
#test.vcf是SNP calling 过滤后生成的vcf 文件；
#p_1_2_3 生成结果的prefix
#1_population.txt是一个文件,包含群体一中所有个体，一般每行一个个体。个体名字要和vcf的名字对应。
#2_population.txt 包含了群体二中所有个体。
#计算的窗口是500kb，而步长是50kb （根据你的需其可以作出调整）。我们也可以只计算每个点的Fst，去掉参数（--fst-window-size 500000 --fst-window-step 50000）即可。

#R-manhattan绘图

library(qq)
data<-read.table("Fst_AxiomGT1.weir.fst",header=T)
pdf(file="Fst_manhattan.pdf",width=20,height=8)
png(file="Fst_manhattan.png",width=2400,height=960)
#sc3 = subset(data,CHROM=="1") 选择任意染色体进行分析
manhattan(data,chr="CHROM",bp="POS",p="WEIR_AND_COCKERHAM_FST",logp=F,main="Manhattan Plot",col = c("#8DA0CB","#E78AC3","#A6D854","#FFD92F","#E5C494","#66C2A5","#FC8D62"),annotatePval = 0.01)
dev.off()

结果展示

6.2 θπ分析

π用来分析碱基多态性，多态性越低，受选择程度越高。

vcftools --vcf AxiomGT1.calls.vcf --keep RESISTANT.list --recode --out AxiomGT1.RESISTANT

vcftools --vcf AxiomGT1.calls.vcf --keep SUSCEPTIBLE.list --recode --out AxiomGT1.SUSCEPTIBLE_500K
#对两个表型进行分类

vcftools --vcf AxiomGT1.RESISTANT.recode.vcf --window-pi 500000 --window-pi-step 500000 --out AxiomGT1.RESISTANT_500K
#以500K为步长分别计算pi值

paste -d '\t' RESISTANT_typename.list AxiomGT1.RESISTANT_500K.windowed.pi > type_RESISTANT_500K.windowed.pi
#单独建立表型标签list，给windowed.pi文件打上表型标签

#R绘图

library(ggplot2)
data<-read.table("allsample_500K.windowed.pi",header=T)
for (i in 1:29){
#result_name = paste("chr",i,".pdf",sep="")
#pdf(result_name,width=15,height=3)
result_name = paste("chr",i,".png",sep="")
png(result_name,width=1500,height=300)
chrom = subset(data,CHROM==i)
xlab = paste("Chromosome",i,"(MB)",sep="")
p <- ggplot(chrom,aes(x=BIN_END/1000000,y=PI,group=Phenotype,color=Phenotype,shape=Phenotype)) + geom_line(size=0.5) + xlab(xlab)+ ylab("Pi") + theme_bw()
print(p)
dev.off()
}

结果展示

七. 全基因组关联分析

全基因组关联分析（Genome-wide association study，GWAS）是一种对全基因组范围内的常见遗传变异（单核苷酸多态性和拷贝数）基因总体关联分析的方法，该方法以自然群体为研究对象， 以长期重组后保留下来的基因（位点）间连锁不平衡（linkage disequilibrium, LD）为基础，将目标性状表型的多样性与基因（或标记位点）的多态性结合起来分析，可直接鉴定出与表型变异密切相 关且具有特定功能的基因位点或标记位点。采用 GWAS 技术在全基因组范围内进行研究，能够一次性对多个性状进行定位，适用于定位性状关联区间、功能基因研究、开发性状选育标记等方面的研究。
可用于分析的软件较多，如GEMMA、Plink或者Tassel。分析方式也分一般线性模式和混合线性模式。这里以GEMMA分析为例。

7.1 性状关联分析

使用软件：GEMMA
源码下载：Releases · genetics-statistics/GEMMA · GitHub
使用流程：

vcftools --vcf AxiomGT1.calls.vcf --plink --out GT1
#生成plink需要的ped/map格式
plink --file GT1 --make-bed --out GT1 --chr-set 30
#转换为bed格式,非人类染色体，使用--chr-set参数设置染色体数目
~/soft/gemma-0.98.1-linux-static -bfile GT1 -gk 2 -p p.txt
~/soft/gemma-0.98.1-linux-static -bfile GT1 -k output/result.sXX.txt -lmm 1 -p p.txt
#-bfile 读取plink的二进制文件（bed\bim\fam）
#-gk 2 标准化的方法计算G矩阵
#-p 读取表型数据（手动生成，不能省掉）

下面使用名为qqman的R包绘制曼哈顿图和qq图，需提前下载安装

library(qqman) # 载入包
data <- read.table("result.assoc.txt",header = TRUE) #读取数据
data1 <- data[,c(1,2,3,12)] #按照规则截取列
data2 <- na.omit(data1) # 删除含有NA的整行
par(cex=0.8) #设置点的大小
#color_set <- rainbow(9)

前面的分析步骤是通用，后面针对分析所有染色体和单个染色体进行了分拆。

所有染色体manhattan图

#pdf(file="GWAS_manhattan.pdf", width=20, height=8)
png(file="GWAS_manhattan.png",width=2000, height=800)
manhattan(data2,chr="chr",snp="rs",bp="ps",p="p_wald",main="Manhattan Plot",col = c("#8DA0CB","#E78AC3","#A6D854","#FFD92F","#E5C494","#66C2A5","#FC8D62"),annotatePval = 0.01) #suggestiveline = FALSE 更加显著
dev.off() # 保存图片

结果展示

单个染色体manhattan图

data <- read.table("result.assoc.txt",header = TRUE) #读取数据
data1 <- data[,c(1,2,3,12)] #按照规则截取列
data2 <- na.omit(data1) # 删除含有NA的整行
par(cex=0.8) #设置点的大小
for (i in 1:29){
result_name = paste("chr",i,".pdf",sep="")
pdf(result_name,width=15,height=6)
#result_name = paste("chr",i,".png",sep="")
#png(result_name,width=900,height=450)
chrom = subset(data2,chr==i)
main_title = paste("Chromosome",i,"(MB)",sep="")
manhattan(chrom,chr="chr",snp="rs",bp="ps",p="p_wald",main=main_title,col = c("#8DA0CB","#E78AC3","#A6D854","#FFD92F","#E5C494","#66C2A5","#FC8D62"),suggestiveline = FALSE, annotatePval = 0.01)
#suggestiveline = FALSE 去除差异显著线
dev.off()
}

结果展示

绘制qq图

pdf(file=“qqplot.pdf”, width=8, height=8)
png(file=“qqplot.png”)
qq(data2$p_wald) # 画qq图
dev.off()
结果展示

7.2 多重假设检验校正

对通过GEMMA得到的result.assoc.txt文件进行切割处理

awk '{print $2,$1,$3,$12}' result.assoc.txt >GWAS.adjust.xls

在excel里根据P值求出校正值，具体操作为对p值升序排序，校正值公式为，FDR=p*m/k ,m为snp总数，k为当前snp排名。
结果展示

7.3 目标性状相关区域基因功能注释

通过gwas4d进行在线注释（之后单独出一篇文章讲解）
结果展示

参考链接：根红苗正的GWAS分析软件：GEMMA_邓飞----育种数据分析之放飞自我-CSDN博客

八. 结语

比起半个月收藏参考的流程，这篇文章又进行了一些补充和完善。然而，从各方面来讲都还有很大的优化空间，所以，这只是一个1.0版本，如果对流程有问题，或者提供更好的建议，欢迎大家加VX：bbplayer2021 （木青）进群交流，备注 申请加入生信交流群。

九. 参考链接

SnpEff使用方法：https://www.jianshu.com/p/a6e46d0c07ee
Tassel使用参考：https://bitbucket.org/tasseladmin/tassel-5-source/wiki/UserManual
Tassel使用手册中文版(2014)：https://wenku.baidu.com/view/4c0cabbe9b6648d7c0c746b2.html
PopLDdecay: 基于 VCF 文件的快速、高效计算连锁不平衡的工具:http://wap.sciencenet.cn/blog-656335-1168505.html
根红苗正的GWAS分析软件：GEMMA_邓飞----育种数据分析之放飞自我-CSDN博客https://blog.csdn.net/yijiaobani/article/details/106215223