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 密度图)
  1. 处理初始数据
    fstqc (质控)
    bowtie2 (构建基因组 index ,及序列比对)
    samtools (构建 samtools 基因组 index , sam , bam 文件转换)
    bcftools (生成 vcf 原始文件)
  1. 标记过滤
    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 参数相互转换

  1. 转换后可以使用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 分析(分为界面版和 Linux 版本,界面版操作比较方便,但是用惯了 Linux 系统的肯定不会选界面版)
    (这里要根据实验目的,比如体色,生长,低氧等)

界面版可以参考 tassel 使用说明书,下面主要讲 Linux 操作

  • 首先要对 vcf 文件进行排序,这里用到的是一个 perl 脚本,不排序后面操作会报错
  • 转换成 hapmp 格式,也可以不转换直接操作,注意后面的参数设置就行
  • 进行关联分析( GLM 和 MLM )
  1. 对 tassel 计算的 GLM 或 MLM 结果进行校正,校正方法 Bonferroni 和 FDR (前者比较绝对,但筛选的结果一定是正确的,后者可能会保留一些有意义的实验结果,看情况使用)

FDR 校正,在 R 中使用 p.adjust 函数进行
Bonferroni 校正比 FDR 严格,得到的有效 SNP 位点会少一些
这里可以参考之前我写的关于两者的区别

  1. 使用 CMplot 包进行结果可视化,曼哈顿图,QQplot(QQplot应该是在校正前做出来还是校正后?)
  2. 筛选有意义的 SNP 位点(包括潜在有意义的位点)
  3. 计算膨胀系数lambda

基因组膨胀因子λ定义为经验观察到的检验统计分布与预期中位数的中值之比,从而量化了因大量膨胀而造成结果的假阳性率。换句话说,λ定义为得到的卡方检验统计量的中值除以卡方分布的预期中值。预期的P值膨胀系数为1,当实际膨胀系数越偏离1,说明存在群体分层的现象越严重,容易有假阳性结果,需要重新矫正群体分层。

(7) 根据 GWAS 结果 找到 QTL区域(这个后面的操作就了解的比较少了,后面学到了会再补充)

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

(8) 验证实验

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

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