karyoploteR：怎么用R來展示你的候選區(qū)域的基因結(jié)構(gòu)以及BigWig信息止潮，活學(xué)活用

首先最重要的參考鏈接：

其次我就是一個(gè)搬運(yùn)工，希望給有需要的小伙伴政基，如有老哥看到侵權(quán)了拌禾，麻煩請(qǐng)私信我刪除，謝謝闽晦！

番外在和小伙伴分享的時(shí)候，小伙伴也發(fā)了一個(gè)可以完成此功能的鏈接 trackViewer Vignette
  • 由于我的傻吊服務(wù)器裝不上包提岔，后面需要linux的基本就是復(fù)制粘貼看了一遍過來的仙蛉。。碱蒙。 話說服務(wù)器裝R包我是真的很無語

  • 那么重要的是通過此文我們可以學(xué)會(huì)繪制出什么樣的圖呢荠瘪？見下：


if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("karyoploteR", version = "3.8")
## 哇夯巷，按照上面安裝得到的是1.8.8版本的，然后會(huì)發(fā)現(xiàn)后面很多命令都沒有哀墓。趁餐。
# 從新安裝
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函數(shù)繪制上面展示區(qū)域內(nèi)的基因。
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) # 第一個(gè)圖層 染色體繪制的區(qū)間
kpPlotGenes(kp, data=genes.data) # 第二個(gè)圖層 區(qū)間內(nèi)各個(gè)基因的結(jié)構(gòu)
  • 從上圖可以看到在這個(gè)區(qū)域內(nèi)的不同轉(zhuǎn)錄本的結(jié)構(gòu)圖篮绰，而不是基因的結(jié)構(gòu)圖后雷，接下來將使用mergeTranscripts函數(shù)將每個(gè)基因?qū)?yīng)的所有轉(zhuǎn)錄本合并成一個(gè)，并且使用addGeneNames加上基因名吠各。使用cex函數(shù)增加染色體的名字大小臀突。
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)
  • 這才是對(duì)于我們展示基因結(jié)構(gòu)的比較好的方式啊。
  • 接下來將使用r0和r1將基因結(jié)構(gòu)圖放置在最下面贾漏，為上面其他數(shù)據(jù)元素留出空間候学。
kp <- plotKaryotype(zoom = TP53.region, cex=2)
kpPlotGenes(kp, data=genes.data, r0=0, r1=0.15, gene.name.cex = 2)
  • 下一步將將基因狀態(tài)HMM結(jié)果導(dǎo)入進(jìn)來，首先使用BiocFileCache 函數(shù)下載下來纵散。然后再使用regioneR包的函數(shù)toGranges導(dǎo)入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列有對(duì)應(yīng)的顏色信息，當(dāng)我們使用kpPlotRegions函數(shù)繪圖時(shí)候就可以使用這里的顏色來代表區(qū)域的顏色伍掀。
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)
注意控制字體大小
  • 接下來我們可以添加一些表觀數(shù)據(jù)了。使用kpPlotBigwig函數(shù)繪制包含BigWig文件的圖片硕盹。此函數(shù)調(diào)用的是rtracklayer’s BigWigFile函數(shù)來導(dǎo)入BigWig文件符匾。
    注意！由于rtracklayer bigwig管理的限制瘩例，kpPlotBigWig在Windows上不起作用啊胶。它僅適用于Linux和Mac計(jì)算機(jī)
  • 首先我們將可視化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文件展示出來的信息峰值太低了垛贤，因?yàn)槲覀冞x取的是默認(rèn)的ymax, 所以會(huì)將y軸的最大值默認(rèn)為全局的最大值焰坪。我們可以通過ymax = "visible.region"來設(shè)置。
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聘惦。但是我們還是缺少一些注釋信息某饰，比如每個(gè)bw文件對(duì)應(yīng)的修飾信息、峰的高度信息善绎。
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 # 獲得所繪制區(qū)域y的最大值
kpAxis(kp, ymin=0, ymax=computed.ymax, r0 = 1, r1 = 1.3) # 設(shè)置Y軸坐標(biāo)值
kpAddLabels(kp, labels = "H3K36me3", r0 = 1, r1 = 1.3, cex=1.6, label.margin = 0.035) # 設(shè)置Y軸標(biāo)簽
  • 我們可以看到H3K4me3的修飾要高于H3K36me3的修飾的黔漂。我們還可以看到，我們已經(jīng)開始重復(fù)代碼禀酱，并且為此使用循環(huán)會(huì)更好炬守，我們將使用autotrack函數(shù)自動(dòng)獲取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循環(huán)和autotrack自動(dòng)跟蹤到位剂跟，我們可以增加組蛋白標(biāo)記的數(shù)量减途，一切都將自動(dòng)調(diào)整酣藻。
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)
}
  • 我們現(xiàn)在可以調(diào)整繪圖參數(shù) 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)
}
  • 我們甚至可以添加其他實(shí)驗(yàn)峰值并使用嵌套自動(dòng)跟蹤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)
}
  • 并添加一個(gè)主標(biāo)題，幾個(gè)額外的標(biāo)簽鳍置，并調(diào)整一些參數(shù)（文字大小等...）辽剧，以獲得更好的最終圖像
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一樣，我們可以更改繪圖區(qū)域（在這種情況下放大）以繪制基因組的任何部分税产，例如怕轿，重疊區(qū)域的詳細(xì)視圖。
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)
}
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
  • 人面猴
    序言：七十年代末砖第，一起剝皮案震驚了整個(gè)濱河市，隨后出現(xiàn)的幾起案子环凿，更是在濱河造成了極大的恐慌梧兼，老刑警劉巖，帶你破解...
    沈念sama閱讀 206,126評(píng)論 6贊 481
  • 死咒
    序言：濱河連續(xù)發(fā)生了三起死亡事件智听，死亡現(xiàn)場(chǎng)離奇詭異羽杰，居然都是意外死亡，警方通過查閱死者的電腦和手機(jī)到推，發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 88,254評(píng)論 2贊 382
  • 救了他兩次的神仙讓他今天三更去死
    文/潘曉璐 我一進(jìn)店門考赛，熙熙樓的掌柜王于貴愁眉苦臉地迎上來，“玉大人莉测，你說我怎么就攤上這事颜骤。” “怎么了捣卤？”我有些...
    開封第一講書人閱讀 152,445評(píng)論 0贊 341
  • 道士緝兇錄：失蹤的賣姜人
    文/不壞的土叔 我叫張陵忍抽，是天一觀的道長。 經(jīng)常有香客問我董朝，道長鸠项，這世上最難降的妖魔是什么？ 我笑而不...
    開封第一講書人閱讀 55,185評(píng)論 1贊 278
  • ?港島之戀（遺憾婚禮）
    正文 為了忘掉前任子姜，我火速辦了婚禮祟绊，結(jié)果婚禮上，老公的妹妹穿的比我還像新娘哥捕。我一直安慰自己牧抽，他們只是感情好，可當(dāng)我...
    茶點(diǎn)故事閱讀 64,178評(píng)論 5贊 371
  • 惡毒庶女頂嫁案：這布局不是一般人想出來的
    文/花漫 我一把揭開白布遥赚。 她就那樣靜靜地躺著阎姥，像睡著了一般。 火紅的嫁衣襯著肌膚如雪鸽捻。 梳的紋絲不亂的頭發(fā)上呼巴，一...
    開封第一講書人閱讀 48,970評(píng)論 1贊 284
  • 城市分裂傳說
    那天泽腮，我揣著相機(jī)與錄音，去河邊找鬼衣赶。 笑死诊赊，一個(gè)胖子當(dāng)著我的面吹牛，可吹牛的內(nèi)容都是我干的府瞄。 我是一名探鬼主播碧磅，決...
    沈念sama閱讀 38,276評(píng)論 3贊 399
  • 雙鴛鴦連環(huán)套：你想象不到人心有多黑
    文/蒼蘭香墨 我猛地睜開眼，長吁一口氣：“原來是場(chǎng)噩夢(mèng)啊……” “哼遵馆！你這毒婦竟也來了鲸郊？” 一聲冷哼從身側(cè)響起，我...
    開封第一講書人閱讀 36,927評(píng)論 0贊 259
  • 萬榮殺人案實(shí)錄
    序言：老撾萬榮一對(duì)情侶失蹤货邓，失蹤者是張志新（化名）和其女友劉穎秆撮，沒想到半個(gè)月后，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體换况，經(jīng)...
    沈念sama閱讀 43,400評(píng)論 1贊 300
  • ?護(hù)林員之死
    正文 獨(dú)居荒郊野嶺守林人離奇死亡职辨，尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 35,883評(píng)論 2贊 323
  • ?白月光啟示錄
    正文 我和宋清朗相戀三年，在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了戈二。 大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片舒裤。...
    茶點(diǎn)故事閱讀 37,997評(píng)論 1贊 333
  • 活死人
    序言：一個(gè)原本活蹦亂跳的男人離奇死亡，死狀恐怖觉吭，靈堂內(nèi)的尸體忽然破棺而出腾供，到底是詐尸還是另有隱情，我是刑警寧澤鲜滩，帶...
    沈念sama閱讀 33,646評(píng)論 4贊 322
  • ?日本核電站爆炸內(nèi)幕
    正文 年R本政府宣布台腥，位于F島的核電站，受9級(jí)特大地震影響绒北，放射性物質(zhì)發(fā)生泄漏黎侈。R本人自食惡果不足惜，卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 39,213評(píng)論 3贊 307
  • 男人毒藥：我在死后第九天來索命
    文/蒙蒙 一闷游、第九天 我趴在偏房一處隱蔽的房頂上張望峻汉。 院中可真熱鬧皮仁，春花似錦坟漱、人聲如沸追驴。這莊子的主人今日做“春日...
    開封第一講書人閱讀 30,204評(píng)論 0贊 19
  • 一樁弒父案，背后竟有這般陰謀
    文/蒼蘭香墨 我抬頭看了看天上的太陽赫悄。三九已至寺酪，卻和暖如春裸准，著一層夾襖步出監(jiān)牢的瞬間梅尤，已是汗流浹背柜思。 一陣腳步聲響...
    開封第一講書人閱讀 31,423評(píng)論 1贊 260
  • 情欲美人皮
    我被黑心中介騙來泰國打工岩调， 沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留，地道東北人赡盘。 一個(gè)月前我還...
    沈念sama閱讀 45,423評(píng)論 2贊 352
  • 代替公主和親
    正文 我出身青樓号枕，卻偏偏與公主長得像，于是被迫代替她去往敵國和親陨享。 傳聞我的和親對(duì)象是個(gè)殘疾皇子葱淳，可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 42,722評(píng)論 2贊 345

推薦閱讀更多精彩內(nèi)容