专栏首页R语言交流中心R语言实现基因组浏览器可视化功能

R语言实现基因组浏览器可视化功能

做生物信息的同仁们应该对基因组浏览器(IGV)都很熟悉,今天给大家介绍下在R语言中如何实现基因组的浏览。首先我们需要用到R包Gviz。其安装主要是用到bioconductor源,如果R版本>3.5可以运行下面的代码:

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

如果R版本<3.5需要运行另一段代码:

source("https://bioconductor.org/biocLite.R")
biocLite("Gviz")

当然你执行晚上面的安装后,可能会报下面错误:

就还需要我们再安装下错误中的包GenomeInfoDb。还是需要bioconductor安装流程。

接下来我们先导入我们需要的样例数据:

library(Gviz)
library(GenomicRanges)
data(cpgIslands)
data(geneModels)

chr <-as.character(unique(seqnames(cpgIslands)))#获取染色体名称
gen <- genome(cpgIslands)#获取参考序列名称

以上就是数据的信息获取,接下来就是如何绘制我们想要的可视化图像:

首先是基础的获取track信息,所用的函数是AnnotationTrack,他可以灵活的去做任何的定位,类似UCSC的定位方式输入的可以是data.frame,IRanges ,GRangesList Irange格式。

atrack <- AnnotationTrack(cpgIslands,name = "CpG")

plotTracks(atrack)

接下来我们增加范围坐标,所用到的函数是GenomeAxisTrack,绘制范围坐标轴。调用实例:

gtrack <- GenomeAxisTrack()
plotTracks(list(gtrack, atrack))

然后还可以在染色体G带图标注染色体的位置,需要用到函数IdeogramTrack,其中主要的参数genome(染色体名称),chromosome(参考序列):

itrack <- IdeogramTrack(genome = gen,chromosome = chr)
plotTracks(list(itrack, gtrack, atrack))

然后就是更加详细的信息的展示,我们需要用到GeneRegionTrack:

我们看下实例:

grtrack <- GeneRegionTrack(geneModels,genome = gen,chromosome = chr, name = "Gene Model")
 plotTracks(list(itrack,gtrack, atrack, grtrack))

如果想注释每个model,那么需要加一个参数transcriptAnnotation = "symbol":

grtrack <- GeneRegionTrack(geneModels,genome = gen,chromosome = chr, name = "Gene Model",transcriptAnnotation = "symbol")

当然,我们也可以只显示一定范围内的models:

plotTracks(list(itrack, gtrack, atrack,grtrack),from = 26700000, to = 26750000)

接下来我们看下如何显示基因参考序列,我们需要用到函数SequenceTrack:

这个函数的使用实例用到了另一个包BSgenome.Hsapiens.UCSC.hg19(源自bioconductor):

library(BSgenome.Hsapiens.UCSC.hg19)
strack <- SequenceTrack(Hsapiens,chromosome = chr)
plotTracks(list(itrack, gtrack, atrack,grtrack, strack), from = 26591822, to = 26591852, cex = 0.8)

最后我们可以把对应的数据添加到我们的图中,比如测序深度等数据。

set.seed(255)
lim <- c(26700000, 26750000)
coords <- sort(c(lim[1], sample(seq(from= lim[1],to = lim[2]), 99), lim[2]))
dat <- runif(100, min = -10, max = 10)
dtrack <- DataTrack(data = dat, start =coords[-length(coords)],end = coords[-1], chromosome = chr, genome = gen,name ="Uniform")


plotTracks(list(itrack, gtrack, atrack,grtrack, dtrack), from = lim[1], to = lim[2])

当然,我们还可以通过type函数进行图形样式的选择默认是散点图,也可以柱状图等具体的选择可以参考下面红框中的参数:

plotTracks(list(itrack, gtrack, atrack,grtrack, dtrack), from = lim[1], to = lim[2], type = "h")

以上都是单组数据的绘图,那么多组数据也是可以进行展示的需要用到参数groups:

data(twoGroups)
dTrack <- DataTrack(twoGroups, name ="uniform")
plotTracks(dTrack, groups =rep(c("control", "treated"),  each = 3), type = c("a","p", "confint"))

当然还有更高级的,那就是多样本的多行展示:

data(dtHoriz)
dtHoriz <- dtHoriz[1:6, ]
plotTracks(dtHoriz, type ="horiz", groups = rownames(values(dtHoriz)),showSampleNames = TRUE,cex.sampleNames = 0.6,separator = 1)

我们还可以发现在IGV中可以在顶部显示测序的峰值,那么如何在此包中显示峰值,我们直接看下实例:

afrom <- 2960000
ato <- 3160000
 alTrack <-AlignmentsTrack(system.file(package = "Gviz","extdata","gapped.bam"), isPaired = TRUE)
 bmt<- BiomartGeneRegionTrack(genome = "hg19", chromosome ="chr12",start = afrom, end = ato, filter = list(with_ox_refseq_mrna =TRUE),stacking = "dense")
 plotTracks(c(bmt, alTrack), from = afrom,to = ato,  chromosome ="chr12")

当然我们也可以不显示下面所有的tracks,通过type来设置:

plotTracks(c(alTrack, bmt), from = afrom,to = ato, chromosome = "chr12", type = "coverage")

它还有个强大的功能那就是给出剪切事件的可视化展示:

plotTracks(c(bmt, alTrack), from = afrom +12700,to = afrom + 15200, chromosome = "chr12", type =c("coverage", "sashimi"))

不仅可以可视化剪切事件,同时还能对指定范围相关的事件进行筛选,通过参数sashimiFilter, sashimiFilterTolerance 。

introns <- GRanges("chr12",IRanges(start = c(2973662,2973919), end = c(2973848, 2974520)))
plotTracks(c(bmt, alTrack), from = afrom +12700,to = afrom + 15200, chromosome = "chr12", type =c("coverage", "sashimi"), sashimiFilter = introns,sashimiFilterTolerance = 5L)

欢迎大家学习交流!

本文分享自微信公众号 - R语言交流中心(R_statistics)

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

原始发表时间:2019-06-26

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

我来说两句

0 条评论
登录 后参与评论

相关文章

  • R(二)近期记录

    这个功能很简单也很常用,但是不加注意还是容易写错,比如只对每一行的前两个元素求和:

    一只羊
  • 我的5年Python7年R,述说她们的差异在哪里?

    首次接触R语言是在2012年读研的时候,有一门课程是统计分析与R语言,清晰地记得期末考试时,由于把答案给同学抄,最终落了个重考的后果(重考92分)。那个时候真的...

    1480
  • 学习Julia与弯道超车

    看一下Julia官网上的Benchmark,Julia综合速度,是R语言的42倍,是Python的15倍,是Java的3倍,是Fortran的1倍,和C语言速度...

    邓飞
  • 使用bowtie2和samblaster一步到位的干净比对

    运行速度很慢,现在有高效工具啦,比如sambamba主要有filter,merge,slice和duplicate等七个功能来处理sam/bam文件,几乎可以替...

    生信技能树
  • R语言作图(一)violin plot

    即便小仙同学决定学习R语言来提升自己作图的“逼格”的时候,心中还有有些疑虑的(嘿嘿,我这么懒,可不愿意做无用功了?)。仔细想了想,貌似又找到了两个学习的理由。

    一只羊
  • SQL注入语法

    字段数为8,我为了好辨识用的数字,我建议用null,来代表数字进行测试,因为数字的兼容性不高,容易出现异常

    天钧
  • 什么是 XLNet ? 为什么它的性能优于 BERT?

    XLNet:NLP领域中一个新的预训练方法,相比BERT可以显著提高20个任务的准确率。

    AI研习社
  • 手把手教你比较两个模型的预测能力

    各位科研芝士的朋友,大家好。最近学习到用NRI进行模型比较,起初当听到NRI这个词的时候,我的表情可能是这样的。

    百味科研芝士
  • R语言 常见函数知识点梳理与解析 | 精选分析

    R语言 控制流:for、while、ifelse和自定义函数function|第5讲

    1480
  • C++ OpenCV之鼠标响应事件

    在OpenCV中也存在鼠标的操作,今天我们先介绍一下鼠标中的操作事件,用于为之后的GrabCut分割来做个前提。

    Vaccae

扫码关注云+社区

领取腾讯云代金券