GWAS分析总体流程
概念
(1)GWAS:
全称“全基因组关联分析”,使用统计模型找到与性状关联的位点,用于分子标记选择(MAS)或者基因定位
(2)GWAS分析的两类性状:
分类性状(阈值性状,质量性状):比如抗病性,颜色等等
质量性状指相对性状的变异呈不连续性,呈现质的中断性变化的性状。由1对或少数几对主基因控制。如鸡羽的芦花斑纹和非芦花斑纹、角的有无、毛色、血型等都属于质量性状。
连续性状(数量性状):比如株高,体重,产量等等(一般是呈现正态分布)
数量性状指相对性状的变异呈连续性,个体之间的差异不明显,很难明确分组。受微效多基因控制,控制数量性状的基因称为数量性状位点(quantitative trait loci, QTLs)。在 QTLs 中, 基因的效应也有大有小。其中, 效应较大的称为主效QTL, 效应较小的称为微效QTL(或微效多基因)。动植物的许多重要经济性状都是数量性状,如作物的产量、成熟期,奶牛的泌乳量,棉花的纤维长度、细度等等。
(3)GWAS的分析方法:
分类性状:logistic回归模型等等
连续性状:GLM,MLM模型等等
分析工具准备
(1) 准备文件
全基因组参考序列
全基因组注释文件
样本重测序文件(双端测序)200个样本左右或以上
(2)各类软件和包
性状分析
R 和 Rstudio
psych R包
lme4 R包
pheatmap R包
reshape2 R包
CMplot R包(绘制 snp 密度图)
处理初始数据
fstqc (质控)
bowtie2 (构建基因组 index ,及序列比对)
samtools (构建 samtools 基因组 index , sam , bam 文件转换)
bcftools (生成 vcf 原始文件)
标记过滤
vcftools linux版
plink linux版
Tassel linux版
分析流程
(1) 表型数据分析处理
将得到的表型数据(一般是数量性状数据)进行分析处理
对数据进行描述性统计,绘制直方图观察数据是否存在不合适的样本数据
检测正态性
剔除不合适的样本
(2) 原始数据处理(识别样本中 snp 位点)
根据全基因组参考序列使用 bowtie2 建立 index 进行比对准备
使用 fastqc 对样本重测序文件进行测序质量检测
使用 fastx_tolltik(或者 trimmomatic ) 去除不良序列
使用 cutadaptor 去掉 read 两端的 adaptor
使用 bowtie2 比对生成 .sam文件
使用 samtools view 将 sam 文件转化为 bam文件
使用 samtools sort 将 .bam 文件进行排序
在 java 环境下使用 picard MarkDuplicates 去除 PCR 冗余 (超声波打断的叫rad需要用Picard处理, 酶切的叫ddrad不用处理可以省略这一步)
使用 bcftools mpileup 获得原始 vcf 文件(这里相当于将所有的序列文件进行合并,里面就含有 snp 的信息)
补充,在做完这些后,可以在 R 中 CMplot 包绘制绘制SNP密度图
(3)基因型过滤(网上教程大多数是从这一步开始,有 vcf 文件,或者将 vcf 文件转换后的 plink 格式的文件,bed,bim,fam 文件)
如果是vfc文件
使用 vcftools 过滤 (关于其参数设置的问题有待研究,一般是过滤掉低于缺失率高于50%的位点,次等位基因深度低于3。在实际中,要更严格,筛除第二等位基因率(次等位基因)频率小于5%(在国际人类基因组单体型图计划(HapMap)中,MAF大于0.05 (5%)的SNP都被作为了调查目标),缺失率大于10%,等位基因的个数要有两个——这个可以调整)
使用 vcftools 将过滤后的文件转换为 plink 格式的文件(或者使用 plink 也可以直接转换),主要得到的是 bed 文件。
plink格式的文件主要有两组五种
ped 和 map 是一组的
bed fam bim 是一组的
两组可通过 plink -recode 参数相互转换
转换后可以使用plink再次过滤(对于计算不同的东西,如群体结构Structure要求位点要少一些筛选的条件就不一样)
(4)群体结构分析
分析之前需要进行第三步的标记筛选,然后根据以下条件去再次筛选:(一般只要按照LD去筛选就可以)
一定的物理距离取一个标记作为代表进行分析
全基因组上抽取一部分标记进行群体结构的分析
按 LD 筛选,LD强度大于一定阈值的标记只保留其中一个用于分析
数据过滤,使用 plink 进行缺失和 maf 筛选
LD 筛选使用 plink 按照 LD 进行筛选
格式转换,然后使用 recode 参数进行转换并得到 str 相关矩阵文件(后续就用该文件进行群体结构分析)(可以根据需求转换成 structure 或者 admixture 格式,structure比较麻烦一些)
利用得出的structure 或者 admixture 格式文件计算出最适 K 值,并在 R 中绘制不同 K 值时最低交叉验证误差的变化
绘制 structure 结构图(有些文章把这个省略掉了,根据文章的侧重不同可选择保留)
PCA分析,计算 PCA 矩阵(还可以通过EIGENSTRAT软件计算主成分,计算各个主成分是否有显著的统计学意义,将P值小于0.05的主成分计算在内),然后绘制 PCA 图
(5)亲缘关系分析
可用于亲缘关系分析的工具有很多,如:GCTA,LDAK,SPAGeDi,EIGENSOFT,TASSEL,现在使用 TASSEL 比较多
GCTA,LDAK 常用于 snp 遗传力估计或者性状遗传力的估计
同样需要前期使用 plink 进行合理筛选
使用相应软件进行亲缘关系矩阵计算(TASSEL 可以使用界面版也可以使用命令行版本)
计算PCA矩阵(还可以通过EIGENSTRAT软件计算主成分,计算各个主成分是否有显著的统计学意义,将P值小于0.05的主成分计算在内),绘制 PCA 图
结果可视化( kinship 值的分布图和 kinship 热图)
(6)关联分析
使用 tassel 进行 GLM/MLM/CMLM 分析(分为界面版和 Linux 版本,界面版操作比较方便,但是用惯了 Linux 系统的肯定不会选界面版)
(这里要根据实验目的,比如体色,生长,低氧等)
界面版可以参考 tassel 使用说明书,下面主要讲 Linux 操作
首先要对 vcf 文件进行排序,这里用到的是一个 perl 脚本,不排序后面操作会报错
转换成 hapmp 格式,也可以不转换直接操作,注意后面的参数设置就行
进行关联分析( GLM 和 MLM )
对 tassel 计算的 GLM 或 MLM 结果进行校正,校正方法 Bonferroni 和 FDR (前者比较绝对,但筛选的结果一定是正确的,后者可能会保留一些有意义的实验结果,看情况使用)
FDR 校正,在 R 中使用 p.adjust 函数进行
Bonferroni 校正比 FDR 严格,得到的有效 SNP 位点会少一些
这里可以参考之前我写的关于两者的区别
使用 CMplot 包进行结果可视化,曼哈顿图,QQplot(QQplot应该是在校正前做出来还是校正后?)
筛选有意义的 SNP 位点(包括潜在有意义的位点)
计算膨胀系数lambda
基因组膨胀因子λ定义为经验观察到的检验统计分布与预期中位数的中值之比,从而量化了因大量膨胀而造成结果的假阳性率。换句话说,λ定义为得到的卡方检验统计量的中值除以卡方分布的预期中值。预期的P值膨胀系数为1,当实际膨胀系数越偏离1,说明存在群体分层的现象越严重,容易有假阳性结果,需要重新矫正群体分层。
(7) 根据 GWAS 结果 找到 QTL区域(这个后面的操作就了解的比较少了,后面学到了会再补充)
利用 R/qtl 软件 MapQTL 软件等定位 QTL
根据 QTL 定位找到相关性状的基因 (这个是用什么软件?)
根据基因的位置和功能来反向验证结果(这里是不是要用富集分析之类的?)
后面还有一连串对 QTL 的分析
(8) 验证实验
验证实验也有很多种,敲除,抑制基因表达,定量PCR等