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)
欢迎大家学习交流！