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等

参考:https://www.jianshu.com/p/c0a032a8e7d5

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 194,491评论 5 459
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 81,856评论 2 371
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 141,745评论 0 319
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 52,196评论 1 263
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 61,073评论 4 355
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 46,112评论 1 272
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 36,531评论 3 381
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 35,215评论 0 253
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 39,485评论 1 290
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 34,578评论 2 309
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 36,356评论 1 326
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 32,215评论 3 312
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 37,583评论 3 299
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 28,898评论 0 17
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,174评论 1 250
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 41,497评论 2 341
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 40,697评论 2 335

推荐阅读更多精彩内容