GWAS 总体流程理解版

自己找了一些文章和视频，先总结了一部分，后面再做补充和实操

一. 相关概念理解

（1）GWAS：

全称“全基因组关联分析”，使用统计模型找到与性状关联的位点，用于分子标记选择（MAS）或者基因定位

（2）GWAS分析的两类性状：

• 分类性状（阈值性状，质量性状）：比如抗病性，颜色等等

质量性状指相对性状的变异呈不连续性，呈现质的中断性变化的性状。由1对或少数几对主基因控制。如鸡羽的芦花斑纹和非芦花斑纹、角的有无、毛色、血型等都属于质量性状。

• 连续性状（数量性状）：比如株高，体重，产量等等（一般是呈现正态分布）

数量性状指相对性状的变异呈连续性，个体之间的差异不明显，很难明确分组。受微效多基因控制，控制数量性状的基因称为数量性状位点(quantitative trait loci, QTLs)。在 QTLs 中, 基因的效应也有大有小。其中, 效应较大的称为主效QTL, 效应较小的称为微效QTL(或微效多基因)。动植物的许多重要经济性状都是数量性状，如作物的产量、成熟期，奶牛的泌乳量，棉花的纤维长度、细度等等。

（3）GWAS的分析方法：

• 分类性状：logistic回归模型等等
• 连续性状：GLM，MLM模型等等

二、准备工作

（1） 准备文件

• 全基因组参考序列
• 全基因组注释文件
• 样本重测序文件（双端测序）200个样本左右或以上

（2）各类软件和包

1. 性状分析
   R 和 Rstudio
   psych R包
   lme4 R包
   pheatmap R包
   reshape2 R包
   CMplot R包（绘制 snp 密度图）

2. 处理初始数据
   fstqc （质控）
   bowtie2 （构建基因组 index ，及序列比对）
   samtools （构建 samtools 基因组 index ， sam , bam 文件转换）
   bcftools （生成 vcf 原始文件）

3. 标记过滤
   vcftools linux版
   plink linux版
   Tassel linux版

三、实验流程

（1） 表型数据分析处理

将得到的表型数据（一般是数量性状数据）进行分析处理

• 对数据进行描述性统计，绘制直方图观察数据是否存在不合适的样本数据
• 检测正态性
• 剔除不合适的样本

（2） 原始数据处理（识别样本中 snp 位点）

1. 根据全基因组参考序列使用 bowtie2 建立 index 进行比对准备
2. 使用 fastqc 对样本重测序文件进行测序质量检测
3. 使用 fastx_tolltik（或者 trimmomatic ） 去除不良序列
4. 使用 cutadaptor 去掉 read 两端的 adaptor
5. 使用 bowtie2 比对生成 .sam文件
6. 使用 samtools view 将 sam 文件转化为 bam文件
7. 使用 samtools sort 将 .bam 文件进行排序
8. 在 java 环境下使用 picard MarkDuplicates 去除 PCR 冗余 （超声波打断的叫rad需要用Picard处理， 酶切的叫ddrad不用处理可以省略这一步）
9. 使用 bcftools mpileup 获得原始 vcf 文件（这里相当于将所有的序列文件进行合并，里面就含有 snp 的信息）

补充，在做完这些后，可以在 R 中 CMplot 包绘制绘制SNP密度图

（3）基因型过滤（网上教程大多数是从这一步开始，有 vcf 文件，或者将 vcf 文件转换后的 plink 格式的文件，bed,bim,fam 文件）

如果是vfc文件

1. 使用 vcftools 过滤 （关于其参数设置的问题有待研究，一般是过滤掉低于缺失率高于50%的位点，次等位基因深度低于3。在实际中，要更严格，筛除第二等位基因率（次等位基因）频率小于5%（在国际人类基因组单体型图计划（HapMap）中，MAF大于0.05 （5%）的SNP都被作为了调查目标），缺失率大于10%,等位基因的个数要有两个——这个可以调整）
2. 使用 vcftools 将过滤后的文件转换为 plink 格式的文件（或者使用 plink 也可以直接转换），主要得到的是 bed 文件。

plink格式的文件主要有两组五种
ped 和 map 是一组的
bed fam bim 是一组的
两组可通过 plink -recode 参数相互转换

3. 转换后可以使用plink再次过滤（对于计算不同的东西，如群体结构Structure要求位点要少一些筛选的条件就不一样）

（4）群体结构分析

分析之前需要进行第三步的标记筛选，然后根据以下条件去再次筛选：（一般只要按照LD去筛选就可以）

• 一定的物理距离取一个标记作为代表进行分析
• 全基因组上抽取一部分标记进行群体结构的分析
• 按 LD 筛选，LD强度大于一定阈值的标记只保留其中一个用于分析

1. 数据过滤，使用 plink 进行缺失和 maf 筛选
2. LD 筛选使用 plink 按照 LD 进行筛选
3. 格式转换，然后使用 recode 参数进行转换并得到 str 相关矩阵文件（后续就用该文件进行群体结构分析）（可以根据需求转换成 structure 或者 admixture 格式，structure比较麻烦一些）
4. 利用得出的structure 或者 admixture 格式文件计算出最适 K 值，并在 R 中绘制不同 K 值时最低交叉验证误差的变化
5. 绘制 structure 结构图（有些文章把这个省略掉了，根据文章的侧重不同可选择保留）
6. PCA分析，计算 PCA 矩阵（还可以通过EIGENSTRAT软件计算主成分，计算各个主成分是否有显著的统计学意义，将P值小于0.05的主成分计算在内），然后绘制 PCA 图

（5）亲缘关系分析

可用于亲缘关系分析的工具有很多，如：GCTA，LDAK，SPAGeDi，EIGENSOFT，TASSEL，现在使用 TASSEL 比较多

GCTA，LDAK 常用于 snp 遗传力估计或者性状遗传力的估计

1. 同样需要前期使用 plink 进行合理筛选
2. 使用相应软件进行亲缘关系矩阵计算（TASSEL 可以使用界面版也可以使用命令行版本）
3. 计算PCA矩阵（还可以通过EIGENSTRAT软件计算主成分，计算各个主成分是否有显著的统计学意义，将P值小于0.05的主成分计算在内），绘制 PCA 图
4. 结果可视化（ kinship 值的分布图和 kinship 热图）

（6）关联分析

1. 使用 tassel 进行 GLM/MLM/CMLM 分析 （这里要根据实验目的，比如体色，生长，低氧等）
2. 对结果进行校正，校正方法 Bonferroni 和 FDR （前者比较绝对，但筛选的结果一定是正确的，后者可能会保留一些有意义的实验结果，看情况使用）
3. 结果可视化，曼哈顿图，QQplot（QQplot应该是在校正前做出来还是校正后？）
4. 筛选有意义的 SNP 位点（包括潜在有意义的位点）
   （这里面有个关于 λ 的计算，其结果和意义还要再找一些相关的东西来看看）

（7） 根据 GWAS 结果 找到 QTL区域（这个后面的操作就了解的比较少了）

1. 利用 R/qtl 软件 MapQTL 软件等定位 QTL
2. 根据 QTL 定位找到相关性状的基因 （这个是用什么软件？）
3. 根据基因的位置和功能来反向验证结果（这里是不是要用富集分析之类的？）
4. 后面还有一连串对 QTL 的分析

（8） 验证实验

验证实验也有很多种，敲除，抑制基因表达，定量PCR等

关于关联分析 QTL 和实验验证后面会再做补充。