设置

首先，我们需要从本地存储库中获取一些R包(在家里我们可以直接使用biocLite安装它们)。

源(“http://192.168.0.9/setRepos.R”)安装。packages(c("h5vc", "h5vcData")) install.packages("BSgenome.Hsapiens.NCBI.GRCh38") # We need this genome to define our reference, you should use a genome that corresponds to the bam files you want to use

通常，HDF5计数文件将存储核苷酸计数数据和随附的样本数据以及用于绘图的函数h5vc包。如果不想为数据创建HDF5文件表示，可以定义一个新函数，以确保数据和示例数据以正确的格式呈现。

首先我们装载所需的包裹，h5vc用于绘图和BSgenome.Hsapiens.NCBI.GRCh38提供参考基因组序列。我们使用BSgenome.Hsapiens.NCBI.GRCh38因为这个例子本演讲我将在这里使用的文件与人类参考基因组的修订是一致的。你应该在这里使用一个合适的基因组对象，例如，如果你想在拟南芥你可能会用BSgenome.Athaliana.TAIR.TAIR9包。

需要(h5vc)要求(BSgenome.Hsapiens.NCBI.GRCh38)

接下来，我们找到示例数据本演讲文件，并创建sampleData对象，该对象只是data.frame描述我们想要处理的样本。

文件<- list。文件(系统。file("extdata"， package = "h5vcData")， "Pt.*bam$") bamFiles <-文件。路径(系统。file("extdata"， package = "h5vcData")， files) sampleData <- data.frame(SampleFiles = files, Sample = sapply(strsplit(files, split = "\\.")， function(x) x[1])， Type = rep(c("Case"， "Control")， length(files)/2)， Column = seq(length(files))， Patient = substr(files, 1,4)， stringsAsFactors = FALSE #傻，但重要的是要得到正确的跟踪标签
样本类型列患者## 1 Pt10Cancer。bam Pt10Cancer病例1 Pt10 ## 2 Pt10Control。bam Pt10Control Control 2 Pt10 ## 3 Pt17Cancerbam Pt17Cancer病例3 Pt17 ## 4 Pt17Control。bam Pt17Control Control 4 Pt17 ## 5 Pt18Cancer。bam Pt18Cancer病例5 Pt18 ## 6 Pt18Control。bam Pt18Control Control 6 Pt18 ## 7 Pt20Cancerbam Pt20Cancer病例7 Pt20 ## 8 Pt20Control。bam Pt20Control Control 8 Pt20 ## 9 Pt23Cancer。bam Pt23Cancer病例9 Pt23 ## 10 Pt23Control。控制10 Pt23 ## 11 Pt25Cancer。bam Pt25Cancer病例11 Pt25 ## 12 Pt25Control。bam Pt25Control Control 12 Pt25

定义绘图函数

最后我们定义plotBAMs函数值为a农庄对象描述我们要绘制的基因组间隔的列表本演讲文件之后,sampleData对象和引用对象(即BSgenome.Hsapiens.NCBI.GRCh38在我们这里)。

plotams <- function(ranges, bamFiles, sampleData, reference){require(h5vc) theData <- tallyRanges(bamFiles, ranges = ranges, reference = reference) for(idx in seq(length(ranges))){#这是需要的，因为我们不使用HDF5来存储数据theData[[idx]]$h5dapplyInfo <- list(Blockstart = start(ranges[idx])， Blockend = end(ranges[idx]))} mismatchPlot(theData, sampleData)}

为DNMT3A基因创建示例数据图

在这里，我们定义了我们想要绘制的基因组范围，在我们的例子中，它们对应于DNMT3A基因的外显子。的农庄当我们想要可视化变量调用时，应该来自一个具有变量调用的VCF文件，但为了这个例子，DNMT3A外显子就足够了。

require(GenomicRanges) dnmt3a <- read.table(system. range)file("extdata"， "dnmt3a.txt"， package = "h5vcData")， header=TRUE, stringsAsFactors = FALSE) dnmt3a <- with(dnmt3a, GRanges(seqname, ranges = IRanges(start = start, end = end))) dnmt3a <-减少(dnmt3a) dnmt3a
##具有29个范围和0个元数据列的GRanges对象:## seqnames ranges strand ##    ## [1] 2 [25227855,25234420] * ## [2] 2 [25235706,25235825] * ## [3] 2 [25236930,25239215] * ## [5] 2 [25239315,25239513] * ## ... ... ... ...## [25] 2 [25313133,25314161] * ## [27] [25328633, 25328744] * ## [27] [25341826,25341885] * ## [29] 2 [25342430,25342590] * ## ------- ## seqinfo:来自未知基因组的1个序列;没有seqlengths

让我们在前4个样本(对象)中绘制一些DNMT3A外显子(2,3,4)Hsapiens定义在BSgenome.Hsapiens.NCBI.GRCh38包装并包含基因组参考序列2)

p <- plotams (ranges = dnmt3a[2:4]， bamFiles = bamFiles[1:4]， sampleData = sampleData[1:4，]， reference = Hsapiens) print(p)

p是一个ggplot对象，因此我们可以通过添加主题命令等来修改它的外观:

打印(p + theme(text = element_text(color = "hotpink")))

嗯，也许把文本变成粉红色不是最好的主意，让我们试试其他的方法:

Print (p + theme(text = element_text(size = 16)， axis.text.x = element_text(angle = 45, size = 12, vjust = 1, hjust = 1)， panel.margin.x = unit(1, units = "cm")))

练习

  1. 应用你的知识来创建一个浏览器应用。这个应用程序应该让你选择哪些样本和外显子(DNMT3A)应该绘制，然后显示选定的区域。
  2. 运用你的知识VariantAnnotation包来实现一个闪亮的应用程序或函数集，加载一个VCF文件，并围绕变量绘制范围。你可能需要一套合适的本演讲文件来执行此操作，即。本演讲包含在VCF文件中引起变量调用的数据的文件。