karyoploteR：怎么用R来展示你的候选区域的基因结构以及BigWig信息，活学活用

首先最重要的参考链接：

其次我就是一个搬运工，希望给有需要的小伙伴，如有老哥看到侵权了，麻烦请私信我删除，谢谢！

番外在和小伙伴分享的时候，小伙伴也发了一个可以完成此功能的链接 trackViewer Vignette
  • 由于我的傻吊服务器装不上包，后面需要linux的基本就是复制粘贴看了一遍过来的。。。 话说服务器装R包我是真的很无语

  • 那么重要的是通过此文我们可以学会绘制出什么样的图呢？见下：


if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("karyoploteR", version = "3.8")
## 哇，按照上面安装得到的是1.8.8版本的，然后会发现后面很多命令都没有。。
# 从新安装
devtools::install_github("bernatgel/karyoploteR") # karyoploteR_1.9.21版本
library(karyoploteR)


TP53.region <- toGRanges("chr17:7,564,422-7,602,719")
> TP53.region 
GRanges object with 1 range and 0 metadata columns:
      seqnames          ranges strand
         <Rle>       <IRanges>  <Rle>
  [1]    chr17 7564422-7602719      *
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths
kp <- plotKaryotype(zoom = TP53.region)
  • 使用kpPlotGenes函数绘制上面展示区域内的基因。
BiocManager::install("TxDb.Hsapiens.UCSC.hg19.knownGene", version = "3.8")
library(TxDb.Hsapiens.UCSC.hg19.knownGene)

> TxDb.Hsapiens.UCSC.hg19.knownGene
TxDb object:
# Db type: TxDb
# Supporting package: GenomicFeatures
# Data source: UCSC
# Genome: hg19
# Organism: Homo sapiens
# Taxonomy ID: 9606
# UCSC Table: knownGene
# Resource URL: http://genome.ucsc.edu/
# Type of Gene ID: Entrez Gene ID
# Full dataset: yes
# miRBase build ID: GRCh37
# transcript_nrow: 82960
# exon_nrow: 289969
# cds_nrow: 237533
# Db created by: GenomicFeatures package from Bioconductor
# Creation time: 2015-10-07 18:11:28 +0000 (Wed, 07 Oct 2015)
# GenomicFeatures version at creation time: 1.21.30
# RSQLite version at creation time: 1.0.0
# DBSCHEMAVERSION: 1.1

genes.data <- makeGenesDataFromTxDb(TxDb.Hsapiens.UCSC.hg19.knownGene,
                                    karyoplot=kp,
                                    plot.transcripts = TRUE, 
                                    plot.transcripts.structure = TRUE)


kp <- plotKaryotype(zoom = TP53.region) # 第一个图层 染色体绘制的区间
kpPlotGenes(kp, data=genes.data) # 第二个图层 区间内各个基因的结构
  • 从上图可以看到在这个区域内的不同转录本的结构图，而不是基因的结构图，接下来将使用mergeTranscripts函数将每个基因对应的所有转录本合并成一个，并且使用addGeneNames加上基因名。使用cex函数增加染色体的名字大小。
genes.data <- addGeneNames(genes.data)
genes.data <- mergeTranscripts(genes.data)

> genes.data
$`genes`
GRanges object with 2 ranges and 2 metadata columns:
        seqnames          ranges strand |     gene_id        name
           <Rle>       <IRanges>  <Rle> | <character> <character>
  55135    chr17 7589389-7606820      + |       55135      WRAP53
   7157    chr17 7565097-7590868      - |        7157        TP53
  -------
  seqinfo: 93 sequences (1 circular) from hg19 genome

$transcripts
$transcripts$`55135`
GRanges object with 1 range and 1 metadata column:
      seqnames          ranges strand |        tx_id
         <Rle>       <IRanges>  <Rle> |  <character>
  [1]    chr17 7589389-7606820      + | 55135.merged
  -------
  seqinfo: 93 sequences (1 circular) from hg19 genome

kp <- plotKaryotype(zoom = TP53.region, cex=2)
kpPlotGenes(kp, data=genes.data)
  • 这才是对于我们展示基因结构的比较好的方式啊。
  • 接下来将使用r0和r1将基因结构图放置在最下面，为上面其他数据元素留出空间。
kp <- plotKaryotype(zoom = TP53.region, cex=2)
kpPlotGenes(kp, data=genes.data, r0=0, r1=0.15, gene.name.cex = 2)
  • 下一步将将基因状态HMM结果导入进来，首先使用BiocFileCache 函数下载下来。然后再使用regioneR包的函数toGranges导入R。
BiocManager::install("BiocFileCache", version = "3.8")
library(BiocFileCache)
bfc <- BiocFileCache(ask=FALSE)

> bfc
class: BiocFileCache
bfccache: C:\Users\ql\AppData\Local\BiocFileCache\BiocFileCache\Cache
bfccount: 0
For more information see: bfcinfo() or bfcquery()

K562.hmm.file <- bfcrpath(bfc, "http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeBroadHmm/wgEncodeBroadHmmK562HMM.bed.gz")
K562.hmm <- toGRanges(K562.hmm.file)
K562.hmm

> K562.hmm
GRanges object with 622257 ranges and 4 metadata columns:
           seqnames              ranges strand |              name     score     itemRgb               thick
              <Rle>           <IRanges>  <Rle> |       <character> <numeric> <character>           <IRanges>
       [1]     chr1         10001-10600      * | 15_Repetitive/CNV         0     #F5F5F5         10001-10600
       [2]     chr1         10601-10937      * | 13_Heterochrom/lo         0     #F5F5F5         10601-10937
       [3]     chr1         10938-11937      * |       8_Insulator         0     #0ABEFE         10938-11937
       [4]     chr1         11938-12337      * | 5_Strong_Enhancer         0     #FACA00         11938-12337
       [5]     chr1         12338-13137      * |   7_Weak_Enhancer         0     #FFFC04         12338-13137
       ...      ...                 ...    ... .               ...       ...         ...                 ...
  [622253]     chrX 155256807-155257806      * |       11_Weak_Txn         0     #99FF66 155256807-155257806
  [622254]     chrX 155257807-155259206      * |       8_Insulator         0     #0ABEFE 155257807-155259206
  [622255]     chrX 155259207-155259406      * |   6_Weak_Enhancer         0     #FFFC04 155259207-155259406
  [622256]     chrX 155259407-155259606      * |   7_Weak_Enhancer         0     #FFFC04 155259407-155259606
  [622257]     chrX 155259607-155260406      * | 15_Repetitive/CNV         0     #F5F5F5 155259607-155260406
  -------
  seqinfo: 23 sequences from an unspecified genome; no seqlengths
  • 可以从上面看到在itemRgb列有对应的颜色信息，当我们使用kpPlotRegions函数绘图时候就可以使用这里的颜色来代表区域的颜色。
kp <- plotKaryotype(zoom = TP53.region, cex=2)
kpPlotGenes(kp, data=genes.data, r0=0, r1=0.15, gene.name.cex = 2)
kpPlotRegions(kp, K562.hmm, col=K562.hmm$itemRgb, r0=0.4, r1=0.5)

kp <- plotKaryotype(zoom = TP53.region, cex=2)
kpPlotGenes(kp, data=genes.data, r0=0, r1=0.15, gene.name.cex = 2)
kpPlotRegions(kp, K562.hmm, col=K562.hmm$itemRgb, r0=0.4, r1=0.5)
kpAddLabels(kp, labels = "Chromatin\nState (HMM)", r0=0.4, r1=0.5, cex = 1)
注意控制字体大小
  • 接下来我们可以添加一些表观数据了。使用kpPlotBigwig函数绘制包含BigWig文件的图片。此函数调用的是rtracklayer’s BigWigFile函数来导入BigWig文件。
    注意！由于rtracklayer bigwig管理的限制，kpPlotBigWig在Windows上不起作用。它仅适用于Linux和Mac计算机
  • 首先我们将可视化H3K4me3信息。绘制在0. 
kp <- plotKaryotype(zoom = TP53.region, cex=2)
kpPlotGenes(kp, data=genes.data, r0=0, r1=0.15, gene.name.cex = 2)
kpPlotRegions(kp, K562.hmm, col=K562.hmm$itemRgb, r0=0.4, r1=0.5)
kpAddLabels(kp, labels = "Chromatin\nState (HMM)", r0=0.4, r1=0.5, cex = 1)

bigwig.file <- "http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeBroadHistone/wgEncodeBroadHistoneK562H3k4me3StdSig.bigWig"

kpPlotBigWig(kp, data=bigwig.file, r0=0.55, r1=1)
  • 可以看到上面bw文件展示出来的信息峰值太低了，因为我们选取的是默认的ymax, 所以会将y轴的最大值默认为全局的最大值。我们可以通过ymax = "visible.region"来设置。
kp <- plotKaryotype(zoom = TP53.region, cex=2)
kpPlotGenes(kp, data=genes.data, r0=0, r1=0.15, gene.name.cex = 2)
kpPlotRegions(kp, K562.hmm, col=K562.hmm$itemRgb, r0=0.4, r1=0.5)
kpAddLabels(kp, labels = "Chromatin\nState (HMM)", r0=0.4, r1=0.5, cex = 1)

bigwig.file <- "http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeBroadHistone/wgEncodeBroadHistoneK562H3k4me3StdSig.bigWig"

kpPlotBigWig(kp, data=bigwig.file, ymax="visible.region", r0=0.55, r1=1)
  • 接下来我们加入H3K36me3的信息
kp <- plotKaryotype(zoom = TP53.region, cex=2)
kpPlotGenes(kp, data=genes.data, r0=0, r1=0.15, gene.name.cex = 2)
kpPlotRegions(kp, K562.hmm, col=K562.hmm$itemRgb, r0=0.4, r1=0.5)
kpAddLabels(kp, labels = "Chromatin\nState (HMM)", r0=0.4, r1=0.5, cex = 1)

H3K4me3.bw <- "http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeBroadHistone/wgEncodeBroadHistoneK562H3k4me3StdSig.bigWig"

kpPlotBigWig(kp, data = H3K4me3.bw, ymax="visible.region", r0=0.55, r1=1)

H3K36me3.bw <- "http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeBroadHistone/wgEncodeBroadHistoneK562H3k36me3StdSig.bigWig"

kpPlotBigWig(kp, data=H3K36me3.bw, ymax="visible.region", r0 = 1, r1 = 1.3)
  • 我们可以看到不同的peak profiler。但是我们还是缺少一些注释信息，比如每个bw文件对应的修饰信息、峰的高度信息。
kp <- plotKaryotype(zoom = TP53.region, cex=2)
kpPlotGenes(kp, data=genes.data, r0=0, r1=0.15, gene.name.cex = 2)
kpPlotRegions(kp, K562.hmm, col=K562.hmm$itemRgb, r0=0.4, r1=0.5)
kpAddLabels(kp, labels = "Chromatin\nState (HMM)", r0=0.4, r1=0.5, cex = 1)

H3K4me3.bw <- "http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeBroadHistone/wgEncodeBroadHistoneK562H3k4me3StdSig.bigWig"

kpPlotBigWig(kp, data = H3K4me3.bw, ymax="visible.region", r0=0.55, r1=1)

computed.ymax <- kp$latest.plot$computed.values$ymax
kpAxis(kp, ymin=0, ymax=computed.ymax, r0=0.55, r1=1)
kpAddLabels(kp, labels = "H3K4me3", r0=0.55, r1=1, cex=1.6, label.margin = 0.035)


H3K36me3.bw <- "http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeBroadHistone/wgEncodeBroadHistoneK562H3k36me3StdSig.bigWig"

kpPlotBigWig(kp, data=H3K36me3.bw, ymax="visible.region", r0 = 1, r1 = 1.3)

computed.ymax <- kp$latest.plot$computed.values$ymax # 获得所绘制区域y的最大值
kpAxis(kp, ymin=0, ymax=computed.ymax, r0 = 1, r1 = 1.3) # 设置Y轴坐标值
kpAddLabels(kp, labels = "H3K36me3", r0 = 1, r1 = 1.3, cex=1.6, label.margin = 0.035) # 设置Y轴标签
  • 我们可以看到H3K4me3的修饰要高于H3K36me3的修饰的。我们还可以看到，我们已经开始重复代码，并且为此使用循环会更好，我们将使用autotrack函数自动获取r0和r1值。
histone.marks <- c(H3K4me3="wgEncodeBroadHistoneK562H3k4me3StdSig.bigWig",
                   H3K36me3="wgEncodeBroadHistoneK562H3k36me3StdSig.bigWig")

base.url <- "http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeBroadHistone/"

kp <- plotKaryotype(zoom = TP53.region, cex=2)
kpPlotGenes(kp, data=genes.data, r0=0, r1=0.15, gene.name.cex = 2)
kpPlotRegions(kp, K562.hmm, col=K562.hmm$itemRgb, r0=0.22, r1=0.3)
kpAddLabels(kp, labels = "Chromatin\nState (HMM)", r0=0.22, r1=0.3, cex=2)

for(i in seq_len(length(histone.marks))) {
  bigwig.file <- paste0(base.url, histone.marks[i])
  at <- autotrack(i, length(histone.marks), r0=0.35, r1=1)
  kp <- kpPlotBigWig(kp, data=bigwig.file, ymax="visible.region", 
                     r0=at$r0, r1=at$r1)
  computed.ymax <- ceiling(kp$latest.plot$computed.values$ymax)
  kpAxis(kp, ymin=0, ymax=computed.ymax, numticks = 2, r0=at$r0, r1=at$r1)
  kpAddLabels(kp, labels = names(histone.marks)[i], r0=at$r0, r1=at$r1, 
              cex=1.6, label.margin = 0.035)
}
  • 一旦我们有for循环和autotrack自动跟踪到位，我们可以增加组蛋白标记的数量，一切都将自动调整。
histone.marks <- c(H3K4me3="wgEncodeBroadHistoneK562H3k4me3StdSig.bigWig",
                   H3K36me3="wgEncodeBroadHistoneK562H3k36me3StdSig.bigWig",
                   H3K27ac="wgEncodeBroadHistoneK562H3k27acStdSig.bigWig",
                   H3K9ac="wgEncodeBroadHistoneK562H3k9acStdSig.bigWig",
                   H3K27me3="wgEncodeBroadHistoneK562H3k27me3StdSig.bigWig")

base.url <- "http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeBroadHistone/"

kp <- plotKaryotype(zoom = TP53.region, cex=2)
kpPlotGenes(kp, data=genes.data, r0=0, r1=0.15, gene.name.cex = 2)
kpPlotRegions(kp, K562.hmm, col=K562.hmm$itemRgb, r0=0.22, r1=0.3)
kpAddLabels(kp, labels = "Chromatin\nState (HMM)", r0=0.22, r1=0.3, cex=2)

for(i in seq_len(length(histone.marks))) {
  bigwig.file <- paste0(base.url, histone.marks[i])
  at <- autotrack(i, length(histone.marks), r0=0.35, r1=1, margin = 0.1)
  kp <- kpPlotBigWig(kp, data=bigwig.file, ymax="visible.region", 
                     r0=at$r0, r1=at$r1)
  computed.ymax <- ceiling(kp$latest.plot$computed.values$ymax)
  kpAxis(kp, ymin=0, ymax=computed.ymax, numticks = 2, r0=at$r0, r1=at$r1)
  kpAddLabels(kp, labels = names(histone.marks)[i], r0=at$r0, r1=at$r1, 
              cex=1.6, label.margin = 0.035)
}
  • 我们现在可以调整绘图参数 plotting parameters 来减少边距和表意文字高度并改变颜色以改善plo的整体外观
pp <- getDefaultPlotParams(plot.type=1)
pp$leftmargin <- 0.15
pp$topmargin <- 15
pp$bottommargin <- 15
pp$ideogramheight <- 5
pp$data1inmargin <- 10

kp <- plotKaryotype(zoom = TP53.region, cex=2, plot.params = pp)
kpAddBaseNumbers(kp, tick.dist = 10000, minor.tick.dist = 2000,
                 add.units = TRUE, cex=1.3, digits = 6)
kpPlotGenes(kp, data=genes.data, r0=0, r1=0.1, gene.name.cex = 2)
kpPlotRegions(kp, K562.hmm, col=K562.hmm$itemRgb, r0=0.15, r1=0.18)
kpAddLabels(kp, labels = "Chromatin\nState (HMM)", r0=0.15, r1=0.18, cex=2)

for(i in seq_len(length(histone.marks))) {
  bigwig.file <- paste0(base.url, histone.marks[i])
  at <- autotrack(i, length(histone.marks), r0=0.23, r1=1, margin = 0.1)
  kp <- kpPlotBigWig(kp, data=bigwig.file, ymax="visible.region",
                     r0=at$r0, r1=at$r1, col = "cadetblue2")
  computed.ymax <- ceiling(kp$latest.plot$computed.values$ymax)
  kpAxis(kp, ymin=0, ymax=computed.ymax, numticks = 2, r0=at$r0, r1=at$r1)
  kpAddLabels(kp, labels = names(histone.marks)[i], r0=at$r0, r1=at$r1, 
              cex=1.6, label.margin = 0.035)
}
  • 我们甚至可以添加其他实验峰值并使用嵌套自动跟踪autotrack来定位它们
base.url <- "http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeBroadHistone/"
histone.marks <- c(H3K4me3="wgEncodeBroadHistoneK562H3k4me3StdSig.bigWig",
                   H3K36me3="wgEncodeBroadHistoneK562H3k36me3StdSig.bigWig",
                   H3K27ac="wgEncodeBroadHistoneK562H3k27acStdSig.bigWig",
                   H3K9ac="wgEncodeBroadHistoneK562H3k9acStdSig.bigWig",
                   H3K27me3="wgEncodeBroadHistoneK562H3k27me3StdSig.bigWig")

DNA.binding <- c(CTCF="wgEncodeBroadHistoneK562CtcfStdSig.bigWig",
                 EZH2="wgEncodeBroadHistoneK562Ezh239875StdSig.bigWig",
                 POL2="wgEncodeBroadHistoneK562Pol2bStdSig.bigWig",
                 P300="wgEncodeBroadHistoneK562P300StdSig.bigWig",
                 HDAC1="wgEncodeBroadHistoneK562Hdac1sc6298StdSig.bigWig",
                 HDAC2="wgEncodeBroadHistoneK562Hdac2a300705aStdSig.bigWig")


pp <- getDefaultPlotParams(plot.type=1)
pp$leftmargin <- 0.15
pp$topmargin <- 15
pp$bottommargin <- 15
pp$ideogramheight <- 5
pp$data1inmargin <- 10

kp <- plotKaryotype(zoom = TP53.region, cex=2, plot.params = pp)
kpAddBaseNumbers(kp, tick.dist = 10000, minor.tick.dist = 2000,
                 add.units = TRUE, cex=1.3, digits = 6)
kpPlotGenes(kp, data=genes.data, r0=0, r1=0.1, gene.name.cex = 2)
kpPlotRegions(kp, K562.hmm, col=K562.hmm$itemRgb, r0=0.15, r1=0.18)
kpAddLabels(kp, labels = "Chromatin\nState (HMM)", r0=0.15, r1=0.18, cex=2)

#Histone marks
total.tracks <- length(histone.marks)+length(DNA.binding)
out.at <- autotrack(1:length(histone.marks), total.tracks, margin = 0.3, r0=0.23)

for(i in seq_len(length(histone.marks))) {
  bigwig.file <- paste0(base.url, histone.marks[i])
  at <- autotrack(i, length(histone.marks), r0=out.at$r0, r1=out.at$r1, margin = 0.1)
  kp <- kpPlotBigWig(kp, data=bigwig.file, ymax="visible.region",
                     r0=at$r0, r1=at$r1, col = "cadetblue2")
  computed.ymax <- ceiling(kp$latest.plot$computed.values$ymax)
  kpAxis(kp, ymin=0, ymax=computed.ymax, numticks = 2, r0=at$r0, r1=at$r1)
  kpAddLabels(kp, labels = names(histone.marks)[i], r0=at$r0, r1=at$r1, 
              cex=1.6, label.margin = 0.035)
}

#DNA binding proteins
out.at <- autotrack((length(histone.marks)+1):total.tracks, total.tracks, margin = 0.3, r0=0.23)

for(i in seq_len(length(DNA.binding))) {
  bigwig.file <- paste0(base.url, DNA.binding[i])
  at <- autotrack(i, length(DNA.binding), r0=out.at$r0, r1=out.at$r1, margin = 0.1)
  kp <- kpPlotBigWig(kp, data=bigwig.file, ymax="visible.region",
                     r0=at$r0, r1=at$r1, col = "darkolivegreen1")
  computed.ymax <- ceiling(kp$latest.plot$computed.values$ymax)
  kpAxis(kp, ymin=0, ymax=computed.ymax, numticks = 2, r0=at$r0, r1=at$r1)
  kpAddLabels(kp, labels = names(DNA.binding)[i], r0=at$r0, r1=at$r1, 
              cex=1.6, label.margin = 0.035)
}
  • 并添加一个主标题，几个额外的标签，并调整一些参数（文字大小等...），以获得更好的最终图像
base.url <- "http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeBroadHistone/"
histone.marks <- c(H3K4me3="wgEncodeBroadHistoneK562H3k4me3StdSig.bigWig",
                   H3K36me3="wgEncodeBroadHistoneK562H3k36me3StdSig.bigWig",
                   H3K27ac="wgEncodeBroadHistoneK562H3k27acStdSig.bigWig",
                   H3K9ac="wgEncodeBroadHistoneK562H3k9acStdSig.bigWig",
                   H3K27me3="wgEncodeBroadHistoneK562H3k27me3StdSig.bigWig")

DNA.binding <- c(CTCF="wgEncodeBroadHistoneK562CtcfStdSig.bigWig",
                 EZH2="wgEncodeBroadHistoneK562Ezh239875StdSig.bigWig",
                 POL2="wgEncodeBroadHistoneK562Pol2bStdSig.bigWig",
                 P300="wgEncodeBroadHistoneK562P300StdSig.bigWig",
                 HDAC1="wgEncodeBroadHistoneK562Hdac1sc6298StdSig.bigWig",
                 HDAC2="wgEncodeBroadHistoneK562Hdac2a300705aStdSig.bigWig")


pp <- getDefaultPlotParams(plot.type=1)
pp$leftmargin <- 0.15
pp$topmargin <- 15
pp$bottommargin <- 15
pp$ideogramheight <- 5
pp$data1inmargin <- 10
pp$data1outmargin <- 0

kp <- plotKaryotype(zoom = TP53.region, cex=3, plot.params = pp)
kpAddBaseNumbers(kp, tick.dist = 10000, minor.tick.dist = 2000,
                 add.units = TRUE, cex=2, tick.len = 3)
kpAddMainTitle(kp, "Epigenetic Regulation in K562", cex=4)
kpPlotGenes(kp, data=genes.data, r0=0, r1=0.1, gene.name.cex = 2.5)
kpPlotRegions(kp, K562.hmm, col=K562.hmm$itemRgb, r0=0.15, r1=0.18)
kpAddLabels(kp, labels = "Chromatin\nState (HMM)", r0=0.15, r1=0.18, cex=2.5)

#Histone marks
total.tracks <- length(histone.marks)+length(DNA.binding)
out.at <- autotrack(1:length(histone.marks), total.tracks, margin = 0.3, r0=0.23)
kpAddLabels(kp, labels = "Histone marks", r0 = out.at$r0, r1=out.at$r1, cex=3.5,
            srt=90, pos=1, label.margin = 0.14)

for(i in seq_len(length(histone.marks))) {
  bigwig.file <- paste0(base.url, histone.marks[i])
  at <- autotrack(i, length(histone.marks), r0=out.at$r0, r1=out.at$r1, margin = 0.1)
  kp <- kpPlotBigWig(kp, data=bigwig.file, ymax="visible.region",
                     r0=at$r0, r1=at$r1, col = "cadetblue2")
  computed.ymax <- ceiling(kp$latest.plot$computed.values$ymax)
  kpAxis(kp, ymin=0, ymax=computed.ymax, tick.pos = computed.ymax, 
         r0=at$r0, r1=at$r1, cex=1.6)
  kpAddLabels(kp, labels = names(histone.marks)[i], r0=at$r0, r1=at$r1, 
              cex=2.2, label.margin = 0.035)
}

#DNA binding proteins
out.at <- autotrack((length(histone.marks)+1):total.tracks, total.tracks, margin = 0.3, r0=0.23)

kpAddLabels(kp, labels = "DNA-binding proteins", r0 = out.at$r0, r1=out.at$r1,
             cex=3.5, srt=90, pos=1, label.margin = 0.14)
for(i in seq_len(length(DNA.binding))) {
  bigwig.file <- paste0(base.url, DNA.binding[i])
  at <- autotrack(i, length(DNA.binding), r0=out.at$r0, r1=out.at$r1, margin = 0.1)
  kp <- kpPlotBigWig(kp, data=bigwig.file, ymax="visible.region",
                     r0=at$r0, r1=at$r1, col = "darkolivegreen1")
  computed.ymax <- ceiling(kp$latest.plot$computed.values$ymax)
  kpAxis(kp, ymin=0, ymax=computed.ymax, tick.pos = computed.ymax, 
         r0=at$r0, r1=at$r1, cex=1.6)
  kpAddLabels(kp, labels = names(DNA.binding)[i], r0=at$r0, r1=at$r1, 
              cex=2.2, label.margin = 0.035)
}
  • 与karyoploteR一样，我们可以更改绘图区域（在这种情况下放大）以绘制基因组的任何部分，例如，重叠区域的详细视图。
TP53.promoter.region <- toGRanges("chr17:7586000-7596000")
kp <- plotKaryotype(zoom = TP53.promoter.region, cex=3, plot.params = pp)
kpAddBaseNumbers(kp, tick.dist = 10000, minor.tick.dist = 2000,
                 add.units = TRUE, cex=2, tick.len = 3)
kpAddMainTitle(kp, "Epigenetic Regulation in K562", cex=4)
kpPlotGenes(kp, data=genes.data, r0=0, r1=0.1, gene.name.cex = 2.5)
kpPlotRegions(kp, K562.hmm, col=K562.hmm$itemRgb, r0=0.15, r1=0.18)
kpAddLabels(kp, labels = "Chromatin\nState (HMM)", r0=0.15, r1=0.18, cex=2.5)

#Histone marks
total.tracks <- length(histone.marks)+length(DNA.binding)
out.at <- autotrack(1:length(histone.marks), total.tracks, margin = 0.3, r0=0.23)
kpAddLabels(kp, labels = "Histone marks", r0 = out.at$r0, r1=out.at$r1, cex=3.5,
            srt=90, pos=1, label.margin = 0.14)

for(i in seq_len(length(histone.marks))) {
  bigwig.file <- paste0(base.url, histone.marks[i])
  at <- autotrack(i, length(histone.marks), r0=out.at$r0, r1=out.at$r1, margin = 0.1)
  kp <- kpPlotBigWig(kp, data=bigwig.file, ymax="visible.region",
                     r0=at$r0, r1=at$r1, col = "cadetblue2")
  computed.ymax <- ceiling(kp$latest.plot$computed.values$ymax)
  kpAxis(kp, ymin=0, ymax=computed.ymax, tick.pos = computed.ymax, 
         r0=at$r0, r1=at$r1, cex=1.6)
  kpAddLabels(kp, labels = names(histone.marks)[i], r0=at$r0, r1=at$r1, 
              cex=2.2, label.margin = 0.035)
}

#DNA binding proteins
out.at <- autotrack((length(histone.marks)+1):total.tracks, total.tracks, margin = 0.3, r0=0.23)

kpAddLabels(kp, labels = "DNA-binding proteins", r0 = out.at$r0, r1=out.at$r1,
             cex=3.5, srt=90, pos=1, label.margin = 0.14)
for(i in seq_len(length(DNA.binding))) {
  bigwig.file <- paste0(base.url, DNA.binding[i])
  at <- autotrack(i, length(DNA.binding), r0=out.at$r0, r1=out.at$r1, margin = 0.1)
  kp <- kpPlotBigWig(kp, data=bigwig.file, ymax="visible.region",
                     r0=at$r0, r1=at$r1, col = "darkolivegreen1")
  computed.ymax <- ceiling(kp$latest.plot$computed.values$ymax)
  kpAxis(kp, ymin=0, ymax=computed.ymax, tick.pos = computed.ymax, 
         r0=at$r0, r1=at$r1, cex=1.6)
  kpAddLabels(kp, labels = names(DNA.binding)[i], r0=at$r0, r1=at$r1, 
              cex=2.2, label.margin = 0.035)
}
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 211,265评论 6赞 490
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 90,078评论 2赞 385
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 156,852评论 0赞 347
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 56,408评论 1赞 283
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 65,445评论 5赞 384
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 49,772评论 1赞 290
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 38,921评论 3赞 406
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 37,688评论 0赞 266
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 44,130评论 1赞 303
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,467评论 2赞 325
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,617评论 1赞 340
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 34,276评论 4赞 329
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 39,882评论 3赞 312
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,740评论 0赞 21
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,967评论 1赞 265
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 46,315评论 2赞 360
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 43,486评论 2赞 348

推荐阅读更多精彩内容