专栏首页实验盒基于RAINBOW的单倍型全基因组关联分析(haplotype-based GWAS)教程

基于RAINBOW的单倍型全基因组关联分析(haplotype-based GWAS)教程

Haplotype-based GWAS(单倍型全基因组关联分析)是基于 haplotype (单倍型)进行的关联分析,在基因组层面寻找与表型相关的变异。

Haplotype 是由一条染色体或线粒体上紧密连锁的多个等位基因组合形成, 每一种组合方式即为一种 haplotype。比如,上图中有 3 个 SNP 和 3 条 染色体,形成了 h1、h2、h3 三种 haplotypes。与单个 SNP 相比,haplotype 提供了更多信息,对定位 causal variant 非常有帮助。

Haplotype-based GWAS 可以用 HAPSTAT、PLINK v1.07、WHAP 等工具进行,不过这些操作通常比较麻烦。另外,haplotypes 数量会随 SNP 数量增大而急剧增大,这会导致检验统计量的自由度非常大,限制了检验的功效。

这里,推荐一个比较新的方法:RAINBOW。它无需 haplotype 的先验信息,把 haplotype 作为 SNP-set (多个 SNPs 组成的集合)处理,分析非常方便。

安装

RAINBOW 这个方法已经集成在 R 包 RAINBOWR 中。RAINBOWR 不仅可以做基于 SNPs 集的GWAS(SNP-set GWAS),也可以做单个 SNP 的 GWAS(Single-SNP GWAS)、分析上位效应(SNP-set x SNP-set interaction)、绘制曼哈顿图和 QQ Plot,功能非常多。

RAINBOWR 稳定版可直接通过 CRAN 安装,最新版则需要通过 GitHub(https://github.com/KosukeHamazaki/RAINBOWR):

#### 稳定版 RAINBOW ####
install.packages("RAINBOWR")  

#### 最新版 RAINBOW ####
install.packages("devtools")  
devtools::install_github("KosukeHamazaki/RAINBOW")

RAINBOWR 依赖其他的一些包,如果安装过程报错,检查一下这些包是否有安装:Rcpp (win 用户则需 Rtools), rgl, tcltk, Matrix, cluster, MASS, pbmcapply, optimx, methods, ape, stringr, pegas, ggplot2, ggtree, scatterpie, phylobase, haplotypes, ggimage

RAINBOWR 安装失败通常是因为个别包自动安装失败,需要按提示手动安装缺失的包。其中, ggtree 需从 Bioconductor 安装:BiocManager::install("ggtree")

数据格式

分析需要三个文件,分别是记录每个个体基因型的文件(geno_score)、基因型位置信息文件(geno_map) 以及表型文件(pheno)。注意,下面的演示例子中,第一行为 header,第一列是行名。

基因型文件

基因型文件 geno_score 需要将每个基因型编码为 -1、0、1 的形式,如果按 additive model 计算的话, -1 代表祖先纯合子,0 代表杂合子,1 代表突变纯合子。每一列为一个个体,每一行为一个 SNP 在不同个体中的基因型:

snp                L1        L2        L3        L4        L5        L6
id1000223         1         1        -1        -1        -1        -1
id1000556        -1        -1         1         1        -1         1
id1000673        -1         1        -1         1        -1         1
id1000830        -1         1         1         1         1         1
id1000955        -1         1         1         1        -1         1
id1001073        -1        -1        -1         1         1        -1

如果不知道怎么转换出基因型文件,可用 plink 的 --recode A 编码成 0、1、2 的方式,用 R 读取结果文件后让每个数字减去 1,再转置一下数据框。

基因型位置信息文件

基因型位置信息文件 geno_map 包含每一个 SNP 的名字、染色体和物理位置:

snp             marker       chr       pos
id1000223 id1000223         1    420422
id1000556 id1000556         1    655693
id1000673 id1000673         1    740153
id1000830 id1000830         1    913806
id1000955 id1000955         1   1041748
id1001073 id1001073         1   1172387

表型文件

表型文件 pheno 每一行为一个人,每一列为一个表型,可以包含多个表型:

id      Flowering.time.at.Arkansas Flowering.time.at.Faridpur 
L1                   75.08333333                         64 
L3                          89.5                         66 
L4                          94.5                         67 
L5                          87.5                         70 
L6                   89.08333333                         73 
L7                           105                       <NA> 

分析

以 RAINBOWR 的实例数据为例:

加载包:

require(RAINBOWR)

读取实例数据:

data("Rice_Zhao_etal")
Rice_geno_score <- Rice_Zhao_etal$genoScore
Rice_geno_map <- Rice_Zhao_etal$genoMap
Rice_pheno <- Rice_Zhao_etal$pheno

过滤 MAF < 0.05 的 SNPs:

x.0 <- t(Rice_geno_score)
MAF.cut.res <- MAF.cut(x.0 = x.0, map.0 = Rice_geno_map)
x <- MAF.cut.res$x
map <- MAF.cut.res$map

选择分析的表型:

trait.name <- "Flowering.time.at.Arkansas"
y <- Rice_pheno[, trait.name, drop = FALSE]

整合转换数据:

modify.data.res <- modify.data(pheno.mat = y, geno.mat = x, map = map,
                               return.ZETA = TRUE, return.GWAS.format = TRUE)
pheno.GWAS <- modify.data.res$pheno.GWAS
geno.GWAS <- modify.data.res$geno.GWAS
ZETA <- modify.data.res$ZETA

进行单个 SNP 的 GWAS,以前 4 个 主成分和亲缘关系 ZETA 进行校正,并绘制曼哈顿图和 QQ 图:

normal.res <- RGWAS.normal(pheno = pheno.GWAS, geno = geno.GWAS,
                           plot.qq = TRUE plot.Manhattan = TRUE
                           ZETA = ZETA, n.PC = 4, P3D = TRUE, count = FALSE)

进行 SNP-Set GWAS,以前 4 个 主成分和亲缘关系 ZETA 进行校正,并绘制曼哈顿图和 QQ 图:

SNP_set.res <- RGWAS.multisnp(pheno = pheno.GWAS, geno = geno.GWAS, ZETA = ZETA, plot.qq = TRUE, plot.Manhattan = TRUE, count = FALSE, n.PC = 4, test.method = "LR", kernel.method = "linear", gene.set = NULL, test.effect = "additive", window.size.half = 5, window.slide = 11)

对于一个给定的 SNP,它对应的 SNP set 大小是 2 * window.size.half + 1,也就是说 SNP set 是以一个 SNP 为中心,上下游各自拓展 window.size.half 个SNPs。参数 window.slide 控制了中心 SNP 滑动的步长,如果要 SNPs 分成 bin,可以按 window.slide = 2 * window.size.half + 1 设定。

在上面的这个例子中,SNP-set大小是 11,滑动步长也是 11 个 SNPs。

查看 SNP-Set GWAS结果,结果第 4 列为 负对数转换后的值 -log10(p):··· See(SNP_set.res$D) ···

SNP-set 也可以通过通过 gene.set 参数自行指定:

gene (or haplotype block)  marker
gene_1  id1000556
gene_1  id1000673
gene_2  id1000830
gene_2  id1000955
gene_2  id1001516
…  …

FAQ

提示 configure: error: gdal-config not found or not executable 报错 需要安装 GDAL(https://gdal.org/index.html),可使用conda安装:conda install -c conda-forge gdal

安装 ggtree 报错 package ‘ggtree’ is not available for this version of R,是因为 ggtree 需要通过 Bioconductor 安装:

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("ggtree")

ggimage 安装过程中报错 ERROR: compilation failed for package ‘magick’,需安装 Magick++。可直接用 conda 安装:conda install -c conda-forge imagemagick

参考

  1. Multi-SNP haplotype analysis methods for association analysis (DOI: 10.1007/978-1-4939-7274-6_24)
  2. https://cran.r-project.org/web/packages/RAINBOWR
  3. https://github.com/KosukeHamazaki/RAINBOWR

本文分享自微信公众号 - 实验盒(gh_8a85afc0b064),作者:实验盒

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2021-01-24

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • Haplotype Reference Consortium:最大规模的单倍型数据库

    在进行基因型填充时,reference panel的选择对填充结果的影响非常大,HapMap包含了3百多万个SNP位点,420个单倍型,1000G包含了8千多万...

    生信修炼手册
  • 使用IMPUTE2进行基因型填充

    需要两个基本元素,第一个是检测样本的分型结果,即图中所示的study genotypes, 第二个元素称之为reference panel, 对应图中的refe...

    生信修炼手册
  • snpQT-又一个人基因组SNP填充和GWAS流程

    发现搜索引擎是个神奇的东西,偶然想起的关键词一搜索,获得的就是意想不到的结果,我以imputation+qc搜索,就找到了snpQT(发音Snip Cute)这...

    用户1075469
  • GWAS中的genotype imputation简介

    GWAS用于寻找与疾病或者特定性状相关联的SNP位点,为了更加有效的挖掘信息,GWAS需要大样本量和高密度的SNP分型结果,最佳的分型方案当然是全基因组测序,然...

    生信修炼手册
  • 采用plink挑选tagSNPs

    tagSNPs叫做标签SNP, 用来代表一组高度连锁不平衡的SNP位点。对于一组高度连锁不平衡的SNP位点而言,在遗传时这些位点往往同时遗传,其包含的信息是冗余...

    生信修炼手册
  • GWAS做完了,下一步做什么?

    通过GWAS分析可以找出与疾病相关联的SNP位点,然而我们的根本目的是找出可能导致疾病发生的SNP位点,这些位点位于GWAS分析结果中,完成了GWAS分析之后,...

    生信修炼手册
  • 使用Minimac进行基因型填充

    Minimac是一款经典的基因型填充软件,该软件也是以内存消耗小,运行速度快而著称,历经了MaCH, minimac, minimac2, minmac3多个版...

    生信修炼手册
  • GWAS和GS的结合:SSGWAS的应用

    满血复活。注意:这个blupf90的新功能,貌似有点问题,好几个人测试显示SSGWAS结果P-value不显示。还未坐实,待我测试后公布。

    邓飞
  • GWAS分析新软件 | GMATs:解析复杂性状和复杂遗传机制的高效工具

    飞哥介绍:终于向作者要来了PPT,在动物大会上听了这个做GWAS的软件,一直想学习,今天作者回复了PPT内容,先分享一下。个人认为这款软件的特色:

    邓飞
  • ANFD-HLA在不同人群中的频率数据库

    在研究SNP时,我们有类似1000G,HapMap, Exac 等数据库,提供了不同人群中的频率信息。对于HLA的研究而言,也有存储频率信息的数据库-ANFD。

    生信修炼手册
  • 笔记 | GWAS 操作流程1:下载数据

    这里,总结一下GWAS的学习笔记,GWAS全称“全基因组关联分析”,使用统计模型找到与性状关联的位点,用于分子标记选择(MAS)或者基因定位,这次学习的教程是p...

    邓飞
  • GWAS ATLAS:最全面的GWAS数据库

    GWAS ATLAS数据库收录了来自4756个人类不同表型的GWAS结果,并进行了不同表型间的遗传相关性分析,对应的文献发表在nature genetics上,...

    生信修炼手册
  • ANNOVAR 软件用法还可以更复杂

    这次耗费15个小时系统性的回顾了该软件,希望可以做到教学上的最佳教程。虽然其它杂七杂八中文教程没有看的必要性,但是其英文文档是需要反复读的。

    生信技能树
  • GWAS大家都耳熟能详, TWAS又是何方神圣

    GWAS称之为全基因组关联分析,是研究复杂疾病遗传易感性的一种方法,已经广泛应用于各种复杂疾病中,识别到了许多与疾病相关的SNP位点,然而GWAS识别到的很多S...

    生信修炼手册
  • 仅2张图分析如何发到顶刊PNAS?

    大家好,今天和大家分享的是2020年3月发表在Proceedings of the National Academy of Sciences o...

    科研菌
  • 想要进行gene prioritization分析,请看这里!

    通过GWAS分析可以识别到与性状关联的SNP位点,然而从生物学角度出发,我们更想了解的是哪些基因或者通路导致了这些位点与性状的关联现象。为了解决这一问题,科学家...

    生信修炼手册
  • 整理了一些自己可能会用到的R包

    用户7010445
  • BAIR最新RL算法超越谷歌Dreamer,性能提升2.8倍

    此次研究的本质在于回答一个问题—使用图像作为观测值(pixel-based)的 RL 是否能够和以坐标状态作为观测值的 RL 一样有效?传统意义上,大家普遍认为...

    机器之心
  • 【直播】我的基因组64:用gwas来预测健康风险

    既然我的基因组已经确定无误了,那么我可以放心的探索啦!正好这几个月收集了朋友圈的一些GWAS新闻,就是某个科研团体研究了几千甚至上万人的基因型,找到了某几个位点...

    生信技能树

扫码关注云+社区

领取腾讯云代金券