尔云间 一个专门做科研的团队

原创 小果 生信果

前言：与复杂疾病相关的遗传因素的研究方法

基于当前流行病学研究的进展及研究热点，小果在这里为大家简单介绍一下关于全基因组关联研究（Genome-wide association study，GWAS）（以下简称GWAS）。全基因组关联研究发展已有二十多年了，研究人员们发现了大量和人类疾病以及其它表型相关联的基因，GWAS是现代遗传学的重要组成部分，推动了孟德尔随机化和多基因风险评分的发展与应用。

基本概念

GWAS是对多个个体在全基因组范围的遗传变异（标记）多态性进行检测，获得基因型，进而将基因型与可观测的性状，即表型，进行群体水平的统计学分析，根据统计量或显著性 p 值筛选出最有可能影响该性状的遗传变异（标记），挖掘与性状变异相关的基因。

相对于连锁分析的优势

关联定位的相对优势：
o 分辨率高（单碱基水平）
o 研究材料来源广泛，可捕获的变异丰富
o 节省时间

材料选择与群体设计

材料选择的基本原则

1）遗传变异和表型变异丰富
2）群体结构分化不能过于明显（如亚种以上,发生生殖隔离是不能做GWAS的）

样本量
非稀有变异中，对中等变异解释率（10%左右）的位点的检测功效要达到80%以上时，需要的样本量在400左右。
位点的效应越低，需要的样本量越大。

关联分析的三要素

o 测定某一群体的表型数据

o 测量该群体的基因型数据

o 进行关联计算

实操

作为生信分析里最基础的技能之一，有人花了很大的时间和精力都不能够完成一次GWAS。

小果在这里用最简单的数据和代码跑一遍GWAS，希望能够帮助大家更好的理解GWAS。

数据：测试一组狗全基因组的遗传变异与分类形状（毛皮颜色）之间的关系。

操作系统：Linux

1 下载样品vcf文件和表型数据

Wget https://de.cyverse.org/dl/d/E0A502CC-F806-4857-9C3A-BAEAA0CCC694/pruned_coatColor_maf_geno.vcf.gz
wget https://de.cyverse.org/dl/d/3B5C1853-C092-488C-8C2F-CE6E8526E96B/coatColor.phen

2 解压VCF文件之后 查看数据

gunzip pruned_coatColor_maf_geno.vcf.gz

#表型数据 前两列也是FID and IID，第三列是表型。

##查看文件发现，这个数据涉及53只小狗的476840个SNP,表型：24只黄毛犬 29只深色毛犬。

3 安装pink和vcftools

###安装plink
wget http://zzz.bwh.harvard.edu/plink/dist/plink-1.07-x86_64.zip
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

当然也可以直接用conda 装

4 将vcf文件转换成map、ped格式，然后转换为Plink二进制格式（fam，bed，bim）

vcftools --vcf pruned_coatColor_maf_geno.vcf --plink --out coatColor
 plink --file coatColor --allow-no-sex --dog --make-bed --noweb --out coatColor.binary

5 候选等位基因列表创建，awk编辑文本

cat pruned_coatColor_maf_geno.vcf | awk 'BEGIN{FS="\t";OFS="\t";}/#/{next;}{{if($3==".")$3=$1":"$2;}print $3,$5;}'  > alt_alleles

6 关联分析

plink --bfile coatColor.binary --make-pheno coatColor.phen "yellow" --assoc --reference-allele alt_alleles --allow-no-sex --adjust --dog --noweb --out coatColor

7 画图

#数据处理
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}')


#R绘图
data=read.table("coatColor.assoc", header=TRUE); data=data[!is.na(data$P),]
bitmap("coatColor_man.bmp", width=20, height=10)
library(qqman)
png(“man.pdf”)
manhattan(data, p = "P", col = c("blue4", "orange3"),suggestiveline = 12,genomewideline = 15,chrlabs = c(1:38, "X"), annotateTop=TRUE, cex = 1.2);
dev.off();

最后的结果图和生成文件，是不是很简单呢？关注小果，下期将为大家带来更多使用生信技巧。

最后小果祝大家学业有成，事事顺心。

推荐阅读

GO的概念及用R做简单的富集分析

R语言ConsensusClusterPlus包无监督聚类

柱状图-肿瘤某一指标的比较和GSVA结果展示

看小果演示CIBETSORT、xcell两种常用的R语言包

小果教你三分钟看懂多条线共存的ROC图的R语言画法