重演一篇傻瓜式GWAS分析|培养兴趣专用 发布于 2020-05-10 18:58
全基因组关联分析（genome-wideassociationstudy,GWAS）应该算是动物生信分析里最热点的技能之一吧。关于GWAS原理的书籍和文章网上有一大堆，这相对于GWAS跑流程来说只能算一座小山丘。因为很多人花了好几个月研究原理，一旦丢去数据分析，还是Error Error Error....我还是秉承那些晦涩难懂的概念在你跑一次程序之后你会理解70%，因此不要等都准备齐了再去着手实施的观点。
做什么事兴趣最重要，反馈的成就感会让你慢慢喜欢上一件事。一切为了兴趣，特做一期“一站式体验GWAS”，你一定可以拿到属于自己的第一张曼哈顿图。
“GWAS是在某一特定人群中研究遗传突变和表型之间的相关性。GWAS的理论基础是连锁不平衡定律（linkage disequilibrium， LD），既假设观察到的SNP与真正的致病突变（causal variant）之间存在很强的LD。
GWAS是一个逻辑清晰的游戏：我们以身高为例，已知身高的遗传率很高（估计为为0.8），也就是说，A身高1.5，B身高1.8，这30厘米的差异中，有至少24厘米是遗传差异造成的。遗传的基础又是DNA，所以我们应该能够在A和B的基因组上，找到DNA的差异，而这种DNA序列的差异最终贡献了这24厘米的差异，GWAS就是想去找到这种DNA序列的差异。
于是我们找来1000个人，量了每个人的身高（表型），测了每个人的基因组，找出所有的基因组上差异的位点（基因型），对于每个差异的位点都去和表型做一个相关性的分析，给出p-value用来衡量表型和每个位点的相关性。如果和某一个位点非常相关，那我们就找到了能够影响身高的DNA差异！
GWAS分析需要准备的材料：软件、基因型数据和表型数据。
表型数据一般是课题研究者前期实验动物群体饲养获得的数据，网上也有很多共享的数据。
表型数据统计分析：逻辑回归（表型数据为二元）；线性回归（表型数据为连续性变量）；表型数据正态分析（如果不是正态分布，需转换处理为正态分布）；表型数据均值、中值、最大值、最小值；影响因子对表型的影响分析。
基因型数据一般是测序下机数据经SNP calling、基因型填补等后拿到的**.vcf文件，后续根据研究者的课题要求和条件进行质控（Quality control），群体结构分层及亲缘关系/PCA等处理。处理好以上文件后，进行最重要的关联分析(当然以上也很重要)，接下来正式重演这个GWAS分析。（质控、call SNP以及表型数据统计分析无；直接画图 相关可能遇到的问题可在简书或知乎私信我，这篇文章在知乎也有发可以复制代码
数据：测试一组狗全基因组的遗传变异与分类形状（毛皮颜色）之间的关系。
#为GWAS分析创建一个工作目录，并下载数据；
#创建工作目录
mkdir ~/GWAS
cd ~/GWAS
#样本VCF文件和表型数据下载
wget https://de.cyverse.org/dl/d/E0A502CC-F806-4857-9C3A-BAEAA0CCC694/pruned_coatColor_maf_geno.vcf.gzwget https://de.cyverse.org/dl/d/3B5C1853-C092-488C-8C2F-CE6E8526E96B/coatColor.pheno
#pheno为处理好的表型文件 vcf为基因组文件
#解压VCF文件之后 可以看看数据
gunzip pruned_coatColor_maf_geno.vcf.gz
#查看vcf
less pruned_coatColor_maf_geno.vcf
#查看表型数据 前两列也是FID and IID，第三列是表型。
less coatColor.pheno
VCF文件全称为Variant Call Format，表示基因组的变异信息，通常为GATK和Samtools软件处理所得到。
查看文件发现，这个数据涉及53只小狗的476840个SNP,表型：24只黄毛犬 29只深色毛犬
##########################正式开始###########################
#数据有了，接下来安装软件，GWAS分析的软件有很多，这里介绍使用主流软件Plink和vcftools来准备文件，R语言统计画图。
#下载安装Plink 不建议使用 #sudo apt-get install Plink 进行下载
cd /usr/local/bin/
sudo wget http://zzz.bwh.harvard.edu/plink/dist/plink-1.07-x86_64.zip
sudo rm -f plink_linux_x86_64.zip
cd plink-1.07-x86_64/
echo export PATH=$PATH:$(pwd) >> ~/.bashrc
source ~/.bashrc
#下载vcftools，将编译后的文件安装到系统
git clone https://github.com/vcftools/vcftools.git
cd vcftools
./autogen.sh
./configure
make
sudo make install
#安装R和RStudio
#安装R
sudo apt-get update && sudo apt-get install -y r-base r-base-dev
#RStudio
sudo apt-get install gdebi-core
wget https://download2.rstudio.org/server/bionic/amd64/rstudio-server-1.2.5033-amd64.deb
sudo gdebi rstudio-server-1.2.5033-amd64.deb
#qqman包适用于画图 使用R- qqman软件统计画图，安装qqman软件包
Rscript -e "install.packages('qqman', contriburl=contrib.url('http://cran.r-project.org/'))"
#接下来正式使用PLINK 准备文件格式 #将VCF文件转换为Plink可读文件格式（map，ped），然后转换为Plink二进制格式（fam，bed，bim）
cd ~/GWAS
vcftools --vcf pruned_coatColor_maf_geno.vcf --plink --out coatColor
plink --file coatColor --allow-no-sex --dog --make-bed --noweb --out coatColor.binary
#此时应该看到以上格式结尾的5个文件（**.map，**.ped，**.fam，**.bed，**.bim）
#候选等位基因列表创建，awk编辑文本
cat pruned_coatColor_maf_geno.vcf | awk 'BEGIN{FS="\t";OFS="\t";}/#/{next;}{{if($3==".")$3=$1":"$2;}print $3,$5;}' > alt_alleles
#文件准备完毕，运行这个简单的关联分析 #关于模型的选择和协变量的调试还需要按照自己的要求学习修改 #plink各选项含义可查看官方文档；
plink --bfile coatColor.binary --make-pheno coatColor.pheno "yellow" --assoc --reference-allele alt_alleles --allow-no-sex --adjust --dog --noweb --out coatColor
#画曼哈顿图
#统计图
unad_cutoff_sug=$(tail -n+2 coatColor.assoc.adjusted | awk '$10>=0.05' | head -n1 | awk '{print $3}')
unad_cutoff_conf=$(tail -n+2 coatColor.assoc.adjusted | awk '$10>=0.01' | head -n1 | awk '{print $3}')
#绘图
Rscript -e 'args=(commandArgs(TRUE));library(qqman);'\
'data=read.table("coatColor.assoc", header=TRUE); data=data[!is.na(data$P),];'\
'bitmap("coatColor_man.bmp", width=20, height=10);'\
'manhattan(data, p = "P", col = c("blue4", "orange3"),'\
'suggestiveline = 12,'\
'genomewideline = 15,'\
'chrlabs = c(1:38, "X"), annotateTop=TRUE, cex = 1.2);'\
'graphics.off();' $unad_cutoff_sug $unad_cutoff_conf
文件coatColor_man.bmp即为结果图，可传输查看。
查询发现：最相关的突变是已知的控制色素生成的MC1R（c.916C > T）中的无义SNP。
congratulations！ 如果做到这一步，我们已经掌握了简单的Linux命令行执行、大体体验了GWAS分析，是否感觉数据分析原理并没有那么难，快趁热学习一下文中没有涉及的质控内容，这将对接下来学习有很大帮助。
这里需要注意在平常动植物中，不使用plink 进行关联分析，可以在数据过滤处理的时候使用，但是在关联分析的时候不使用。一般是在人类 GWAS 才会使用 plink 进行关联分析。因为它没办法实现复杂模型，就是 MLM 等模型。
在动植物中关联分析 tassel 使用的最多。 上百万标记，几百个样本要几十G 上百G 内存。
gapit 主要是基于 R 软件 。
