m6ASeq RNA甲基化测序ChIP-seqATACSeq 开放染色质分析

使用 edgerTMM 算法对 bw 文件的均一化并且根据 bi

2019-08-23  本文已影响40人  热衷组培的二货潜

链接来源

一切版权来源参考链接:此处只是用来记录和不能出去的情况。

涉及脚本

#### A minimal example on how to use EnrichedHeatmap together with rtracklayer:

require(EnrichedHeatmap)
require(rtracklayer)
require(circlize)
require(data.table)

## We start from a BED file with coordinates and load as GRanges:
tmp.targets <- makeGRangesFromDataFrame(
  df = fread("~/your.bed", header = F),
  seqnames.field = "V1", start.field = "V2", end.field = "V3")

## Say we want to take the peak center and extend it by 5kb in each direction:
tmp.extension <- 5000

## Extend center of the peaks by tmp.extension in each direction:
tmp.targets_extended <- resize(tmp.targets, fix = "center", width = tmp.extension*2)

## Now load the content of the bigwig limited to the regions we are interested in.
## This is much quicker than loading the entire bigwig and does not consume so much memory:
tmp.bigwig <- rtracklayer::import("~/your.bigwig" , 
                                  format = "BigWig", 
                                  selection = BigWigSelection(tmp.targets_extended))
    
## create the normalizedMatrix that EnrichedHeatmap accepts as input.
## We use the tmp.targets center (width=1) because from what I understand normalizeMatrix
## does not allow to turn off its extend= option. Therefore we trick it by simply
## providing the peak centers and then let the function extend it by our predefined window size.
normMatrix <- normalizeToMatrix(signal = tmp.bigwig, 
                                target = resize(tmp.targets, fix = "center", width = 1), 
                                background = 0, 
                                keep = c(0, 0.99),      ## minimal value to the 99th percentile
                                target_ratio = 0,
                                mean_mode = "w0",       ## see ?EnrichedHeatmap on other options
                                value_column = "score", ## = the name of the 4th column of the bigwig
                                extend = tmp.extension)

## a color gradient that I personally find visually appealing, which will cover
## the range from the lowest value of normMatrix to the 99th percentile
## (99th perc. avoids extreme values skewing the heatmap):
col_fun = circlize::colorRamp2(quantile(normMatrix, c(0, .99)), c("darkblue", "darkgoldenrod1"))

## heatmap function:
enrHtmp <- EnrichedHeatmap( mat = normMatrix, 
                            pos_line = FALSE, ## no dashed lines around the start
                            border = FALSE,   ## no box around heatmap
                            col = col_fun,    ## color gradients from above
                            column_title = "Nice Heatmap", ## column title 
                            column_title_gp = gpar(fontsize = 15, fontfamily = "sans"),
                            ## these three options produce a high-quality pdf
                            ## while keeping the file size small so that it easily fits
                            ## nto any powerpoint presentation without crashing it
                            use_raster = TRUE, raster_quality = 10, raster_device = "png",
                            ## turn off background colors
                            rect_gp = gpar(col = "transparent"), 
                            ## legend:
                            heatmap_legend_param = list(
                              legend_direction = "horizontal", ## legend horizontal
                              title = "legend_title"),
                            ## options for the profile plot on top
                            top_annotation = HeatmapAnnotation(
                              enriched = anno_enriched(
                                gp = gpar(col = "black", lty = 1, lwd=2),
                                col="black")
                            )
) ## end of EnrichedHeatmap function
    
## Instead of plotting to the Rstudio device save as pdf,
## with width-2 and height-6 I personally find the heatmap visually most appealing,
## it looks good while not being too "fat":
pdf("~/EnrichedHeatmap.pdf", width = 2, height = 6)

## Plot it:
draw(enrHtmp,                            ## plot the heatmap from above 
     heatmap_legend_side = "bottom",     ## we want the legend below the heatmap
     annotation_legend_side = "bottom",  ## 
     padding = unit(c(4, 4, 4, 4), "mm") ## some padding to avoid labels beyond plot borders
    )

dev.off() ## close the pdf
https://twitter.com/ATpoint90/status/1162065802826342407
上一篇下一篇

猜你喜欢

热点阅读