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