自己找了一些文章和视频,先总结了一部分,后面再做补充和实操
一. 相关概念理解
(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等