GWAS 总体流程理解版
2020-05-30 本文已影响0人
奔跑的Forrest
自己找了一些文章和视频，先总结了一部分，后面再做补充和实操
一. 相关概念理解
（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 分析 （这里要根据实验目的，比如体色，生长，低氧等）
- 对结果进行校正，校正方法 Bonferroni 和 FDR （前者比较绝对，但筛选的结果一定是正确的，后者可能会保留一些有意义的实验结果，看情况使用）
- 结果可视化，曼哈顿图，QQplot（QQplot应该是在校正前做出来还是校正后？）
- 筛选有意义的 SNP 位点（包括潜在有意义的位点）
（这里面有个关于 λ 的计算，其结果和意义还要再找一些相关的东西来看看）
（7） 根据 GWAS 结果 找到 QTL区域（这个后面的操作就了解的比较少了）
- 利用 R/qtl 软件 MapQTL 软件等定位 QTL
- 根据 QTL 定位找到相关性状的基因 （这个是用什么软件？）
- 根据基因的位置和功能来反向验证结果（这里是不是要用富集分析之类的？）
- 后面还有一连串对 QTL 的分析
（8） 验证实验
验证实验也有很多种，敲除，抑制基因表达，定量PCR等
