R/genelist_specific_profileplot.R
In sethiyap/FungalSporeAnalysis: Genomics data analysis of fungal spores.
Defines functions
genelist_specific_profileplot
Documented in
genelist_specific_profileplot
#' profile plot specific to a genelist
#' @description Profile plot for list of genes to be subsetted from genome
#' feature file.
#' @author pooja sethiya
#' @param feature_txDb A TxDb object of gene feature file, generated by
#' \code{GenomicFeatures::makeTxDbFromGFF()}
#' @param bw_files A character vector containing either name or names of GRanges
#' object of bw or bdg files. GRanges object for bw_files can be generated by: \code{
#' rtracklayer::import()}
#' @param genelist A tibble with two columns \code{1. gene_id and 2.Expression
#' value} without header, to subset genes from the feature file.
#' \code{default:NULL}
#' @param tss Logical, if TRUE plot 1kb up/downstream of start co-ordinate; if
#' FALSE plot 1kb up/downstream of gene-body \code{default: TRUE}
#' @param palette A character vector to choose sequential palette; \itemize{
#' \item white_green \item white_red \item white_blue \item cream_pink \item
#' cream_brown \item cream_green } \code{default: cream_brown}
#' @param max_key Numeric, maximum limit of the color key.\code{default:6}
#' @param min_key Numeric, mimimum limit of the color key.\code{default:0}
#' @param top_line Logical, whether to plot profile with average signal as line
#' plot at the top of plot. \code{default: TRUE}
#' @param ymax Numeric, y-axis maximum limit for lineplot \code{default:6}
#' @param ymin Numeric, y-axis minimum limit for lineplot \code{default:1}
#' @param log2 Logical, whether to plot profile as log-transformed or not.
#' @param rename Logical, whether to remove characters after first \code{- OR
#' _}, suitable for long names of bw objects.\code{default: FALSE}
#' @param output Logical, to return a output of profile plot or not.
#' \code{default: TRUE, with output_name,if FALSE; plot will be returned} **
#' Highly recommended to plot in outfile when .bw files are more than 2.
#' @param output_name A character vector containing name of the output
#' plot.\code{default:Sample}
#'
#' @export
#' @import magrittr
#' @import png
#' @importFrom rtracklayer import.bed
#' @importFrom GenomicFeatures promoters
#' @importFrom broom tidy
#' @importFrom dplyr mutate
#' @importFrom purrr map
#' @import EnrichedHeatmap
#' @import ComplexHeatmap
#' @importFrom grid unit
#' @importFrom grid gpar
#'
#' @return Parallel plots of profiles(log2) for given gene lists, ordered by
#' given expression value (high to low).
#'
#' @examples
#' \dontrun{
#'
#' # Load feature file from system data, to run example data
#' feature_txDb <- AnnotationDbi::loadDb(system.file("extdata/sqllite/an_feature_file_s10_m04_r07.sqlite" , package = "FungalSporeAnalysis"))
#'
#' # Or load feature file as TxDb
#' feature_txDb <- GenomicFeatures::makeTxDbFromGFF("an_feature_file_s10_m04_r07.gff")
#'
#' # Load bw files as GRanges object
#' pol2_veA_wt_spore <- rtracklayer::import.bw("pol2_veA_wt_spore_mix22_CACAGTTGGT_normalized_repeat.bw")
#' TBP_veA_wt_spore <- rtracklayer::import.bw("TBP_veA_wt_spore_mix22_CAGTTGGT_normalized_repeat.bw")
#'
#' # Load GRanges object from system data, to run example data
#' data_files <- system.file("extdata/sysdata.rda" , package = "FungalSporeAnalysis")
#' load(data_files)
#'
#' # Load example gene_list file
#' gene_list <- readr::read_delim(system.file("extdata/genesets/an_spore_pol2.txt" , package = "FungalSporeAnalysis"), delim="\t", col_names = FALSE)
#'
#' #################
#'
#' bw_files <- c("pol2_veA_wt_spore","TBP_veA_wt_spore")
#' genelist_specific_profileplot(feature_txDb=feature_txDb,bw_files = bw_files, genelist=gene_list, output_name="An_Pol2_TBP_TFIIB", ymin=3)
#'
#' }
genelist_specific_profileplot <- function(feature_txDb,bw_files,genelist=NULL,tss=TRUE,palette="cream_brown",max_key=6, min_key=0, top_line=TRUE,ymax=6,ymin=1,log2=TRUE,rename=FALSE,output=TRUE,output_name="Sample"){
if(is.null(genelist)==TRUE){
sub_feature_gr <- GenomicFeatures::genes(feature_txDb)
row_order <- rlang::expr(order(enriched_score(eval(mat1)), decreasing = TRUE))
}
else{
feature_gr <- GenomicFeatures::genes(feature_txDb)
gene_list <- genelist %>% dplyr::arrange(desc(X2))
#--- order by the expression value
sub_feature_gr <- subset(feature_gr,feature_gr$gene_id %in% gene_list$X1)
sub_feature_gr <- sub_feature_gr[match(gene_list$X1,sub_feature_gr$gene_id),]
row_order <- NULL
}
if(tss==TRUE){
tss = GenomicFeatures::promoters(sub_feature_gr, upstream = 0, downstream = 1)
axis_name = c("-1kb","Start","+1kb")
}
else{
tss=sub_feature_gr
axis_name = c("-1kb","Start","Stop","+1kb")
}
## prepare signal data
if(rename==TRUE){
bw_files <- broom::tidy(bw_files) %>%
dplyr::mutate(names=x,x=factor(x, levels=x)) %>%
dplyr::mutate(names=paste(gsub(pattern = "-.*|_.*", replacement="",x=names)))
}
else{
bw_files <- broom::tidy(bw_files) %>%
dplyr::mutate(names=x,x=factor(x, levels=x))
}
print(bw_files)
## generate normalised matrix in tidy way
xx <- bw_files %>%
dplyr::mutate(bw = purrr::map(x, function(ii) {
get(as.character(ii))
})) %>%
dplyr::mutate(norm_matrix = purrr::map(bw, function(ii) {
nn <- EnrichedHeatmap::normalizeToMatrix(ii, tss, value_column = "score",background = 0,
extend = c(1000),w = 50,
smooth = TRUE)
nn[nn<0]=0
return(nn)
}))
if(log2==TRUE){
mat1 <- rlang::expr(log2(x[[1]]+1))
mat2 <- rlang::expr(log2(x[[i]]+1))
}
else{
mat1 <- rlang::expr(x[[1]])
mat2 <- rlang::expr(x[[i]])
}
# for color palette
if(palette=="white_green"){
colors <- c("#f7fcfd","#e5f5f9","#ccece6","#99d8c9", "#41ae76","#00441b", "#004529")
}
if(palette=="white_red"){
colors <- c("#f7fcfd","#fff7f3","#fff7ec","#fee0d2","#ef6548","#d7301f","#b30000", "#7f0000")
}
if(palette=="white_blue"){
colors = c("#eff3ff","#c6dbef","#9ecae1","#6baed6","#2171b5","#08306b","#253494", "#081d58")
}
if(palette=="cream_pink"){
colors <- c("#feebe2","#fcc5c0","#fa9fb5","#f768a1", "#c51b8a","#7a0177", "#49006a")
}
if(palette=="cream_green"){
colors <- c("#f7fcb9","#d9f0a3","#addd8e","#78c679","#41ab5d", "#238443","#006837", "#006837")
}
if(palette=="cream_brown"){
colors <- c("#ffeda0","#fed976","#feb24c","#fd8d3c","#fc4e2a","#e31a1c","#bd0026","#800026")
}
split_factor <- round((max_key/6),1)
breaks = seq(min_key,max_key, by = split_factor)
if(top_line==TRUE){
top_annotation = ComplexHeatmap::HeatmapAnnotation(
lines = EnrichedHeatmap::anno_enriched(gp = grid::gpar(fontsize=12),
ylim=c(ymin, ymax),
# yaxis_side = "right",
# yaxis_facing = "inside",
# yaxis_gp = grid::gpar(fontsize = 10, lwd=1.5)),
axis_param = list(gp=grid::gpar(fontsize = 10, lwd=1.5),
side="right",
facing="inside")))
}
else{
top_annotation = NULL
}
get_enrichment_heatmap_list <- function(x, names, titles, ...) {
ll <- length(x)
## first heatmap
ehml <- EnrichedHeatmap::EnrichedHeatmap(mat = eval(mat1), name = names[[1]], column_title = titles[[1]], show_heatmap_legend = T,
row_order= eval(row_order),
use_raster = TRUE, ...)
## several other heatmaps if length of x > 1.
if (ll > 1) {
for (i in 2:ll) {
print(i)
ehml <- ehml +
EnrichedHeatmap::EnrichedHeatmap(
mat = eval(mat2),
row_order= eval(row_order),
name = ifelse(length(names) >= i, names[i], "NA"),
use_raster = TRUE,
column_title = ifelse(length(titles) >= i, titles[i], "NA"),
show_heatmap_legend = ifelse(length(names) >= i, TRUE, FALSE), ...
) ## legend will be shown only if the name is given for a heatmap.
}
}
return(ehml)
}
ehm_list <- get_enrichment_heatmap_list(x = xx$norm_matrix,names = xx$names,titles = xx$names,
cluster_rows = FALSE,
pos_line = TRUE,
show_row_names = FALSE,
axis_name_rot = 90,
heatmap_legend_param = list(color_bar = "continuous",legend_direction="horizontal"),
axis_name = axis_name,
col = circlize::colorRamp2(breaks = breaks,
colors = as.vector(colors[1:length(breaks)])),
top_annotation = top_annotation
)
if(output==TRUE){
print("plotting")
png(file=paste(output_name, length(sub_feature_gr),"hm.png", sep="_"),width=nrow(bw_files)*1.5,height=5.5,pointsize = 8, res=300,units = "in")
ComplexHeatmap::draw(ehm_list, heatmap_legend_side = "top", gap = grid::unit(1.5, "mm"))
dev.off()
}
else{
return(ComplexHeatmap::draw(ehm_list, heatmap_legend_side = "right", gap = grid::unit(1.5, "mm")))
}
}
sethiyap/FungalSporeAnalysis documentation built on July 31, 2020, 10:26 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions genelist_specific_profileplot
Documented in genelist_specific_profileplot
#' profile plot specific to a genelist
#' @description Profile plot for list of genes to be subsetted from genome
#' feature file.
#' @author pooja sethiya
#' @param feature_txDb A TxDb object of gene feature file, generated by
#' \code{GenomicFeatures::makeTxDbFromGFF()}
#' @param bw_files A character vector containing either name or names of GRanges
#' object of bw or bdg files. GRanges object for bw_files can be generated by: \code{
#' rtracklayer::import()}
#' @param genelist A tibble with two columns \code{1. gene_id and 2.Expression
#' value} without header, to subset genes from the feature file.
#' \code{default:NULL}
#' @param tss Logical, if TRUE plot 1kb up/downstream of start co-ordinate; if
#' FALSE plot 1kb up/downstream of gene-body \code{default: TRUE}
#' @param palette A character vector to choose sequential palette; \itemize{
#' \item white_green \item white_red \item white_blue \item cream_pink \item
#' cream_brown \item cream_green } \code{default: cream_brown}
#' @param max_key Numeric, maximum limit of the color key.\code{default:6}
#' @param min_key Numeric, mimimum limit of the color key.\code{default:0}
#' @param top_line Logical, whether to plot profile with average signal as line
#' plot at the top of plot. \code{default: TRUE}
#' @param ymax Numeric, y-axis maximum limit for lineplot \code{default:6}
#' @param ymin Numeric, y-axis minimum limit for lineplot \code{default:1}
#' @param log2 Logical, whether to plot profile as log-transformed or not.
#' @param rename Logical, whether to remove characters after first \code{- OR
#' _}, suitable for long names of bw objects.\code{default: FALSE}
#' @param output Logical, to return a output of profile plot or not.
#' \code{default: TRUE, with output_name,if FALSE; plot will be returned} **
#' Highly recommended to plot in outfile when .bw files are more than 2.
#' @param output_name A character vector containing name of the output
#' plot.\code{default:Sample}
#'
#' @export
#' @import magrittr
#' @import png
#' @importFrom rtracklayer import.bed
#' @importFrom GenomicFeatures promoters
#' @importFrom broom tidy
#' @importFrom dplyr mutate
#' @importFrom purrr map
#' @import EnrichedHeatmap
#' @import ComplexHeatmap
#' @importFrom grid unit
#' @importFrom grid gpar
#'
#' @return Parallel plots of profiles(log2) for given gene lists, ordered by
#' given expression value (high to low).
#'
#' @examples
#' \dontrun{
#'
#' # Load feature file from system data, to run example data
#' feature_txDb <- AnnotationDbi::loadDb(system.file("extdata/sqllite/an_feature_file_s10_m04_r07.sqlite" , package = "FungalSporeAnalysis"))
#'
#' # Or load feature file as TxDb
#' feature_txDb <- GenomicFeatures::makeTxDbFromGFF("an_feature_file_s10_m04_r07.gff")
#'
#' # Load bw files as GRanges object
#' pol2_veA_wt_spore <- rtracklayer::import.bw("pol2_veA_wt_spore_mix22_CACAGTTGGT_normalized_repeat.bw")
#' TBP_veA_wt_spore <- rtracklayer::import.bw("TBP_veA_wt_spore_mix22_CAGTTGGT_normalized_repeat.bw")
#'
#' # Load GRanges object from system data, to run example data
#' data_files <- system.file("extdata/sysdata.rda" , package = "FungalSporeAnalysis")
#' load(data_files)
#'
#' # Load example gene_list file
#' gene_list <- readr::read_delim(system.file("extdata/genesets/an_spore_pol2.txt" , package = "FungalSporeAnalysis"), delim="\t", col_names = FALSE)
#'
#' #################
#'
#' bw_files <- c("pol2_veA_wt_spore","TBP_veA_wt_spore")
#' genelist_specific_profileplot(feature_txDb=feature_txDb,bw_files = bw_files, genelist=gene_list, output_name="An_Pol2_TBP_TFIIB", ymin=3)
#'
#' }
genelist_specific_profileplot <- function(feature_txDb,bw_files,genelist=NULL,tss=TRUE,palette="cream_brown",max_key=6, min_key=0, top_line=TRUE,ymax=6,ymin=1,log2=TRUE,rename=FALSE,output=TRUE,output_name="Sample"){
if(is.null(genelist)==TRUE){
sub_feature_gr <- GenomicFeatures::genes(feature_txDb)
row_order <- rlang::expr(order(enriched_score(eval(mat1)), decreasing = TRUE))
}
else{
feature_gr <- GenomicFeatures::genes(feature_txDb)
gene_list <- genelist %>% dplyr::arrange(desc(X2))
#--- order by the expression value
sub_feature_gr <- subset(feature_gr,feature_gr$gene_id %in% gene_list$X1)
sub_feature_gr <- sub_feature_gr[match(gene_list$X1,sub_feature_gr$gene_id),]
row_order <- NULL
}
if(tss==TRUE){
tss = GenomicFeatures::promoters(sub_feature_gr, upstream = 0, downstream = 1)
axis_name = c("-1kb","Start","+1kb")
}
else{
tss=sub_feature_gr
axis_name = c("-1kb","Start","Stop","+1kb")
}
## prepare signal data
if(rename==TRUE){
bw_files <- broom::tidy(bw_files) %>%
dplyr::mutate(names=x,x=factor(x, levels=x)) %>%
dplyr::mutate(names=paste(gsub(pattern = "-.*|_.*", replacement="",x=names)))
}
else{
bw_files <- broom::tidy(bw_files) %>%
dplyr::mutate(names=x,x=factor(x, levels=x))
}
print(bw_files)
## generate normalised matrix in tidy way
xx <- bw_files %>%
dplyr::mutate(bw = purrr::map(x, function(ii) {
get(as.character(ii))
})) %>%
dplyr::mutate(norm_matrix = purrr::map(bw, function(ii) {
nn <- EnrichedHeatmap::normalizeToMatrix(ii, tss, value_column = "score",background = 0,
extend = c(1000),w = 50,
smooth = TRUE)
nn[nn<0]=0
return(nn)
}))
if(log2==TRUE){
mat1 <- rlang::expr(log2(x[[1]]+1))
mat2 <- rlang::expr(log2(x[[i]]+1))
}
else{
mat1 <- rlang::expr(x[[1]])
mat2 <- rlang::expr(x[[i]])
}
# for color palette
if(palette=="white_green"){
colors <- c("#f7fcfd","#e5f5f9","#ccece6","#99d8c9", "#41ae76","#00441b", "#004529")
}
if(palette=="white_red"){
colors <- c("#f7fcfd","#fff7f3","#fff7ec","#fee0d2","#ef6548","#d7301f","#b30000", "#7f0000")
}
if(palette=="white_blue"){
colors = c("#eff3ff","#c6dbef","#9ecae1","#6baed6","#2171b5","#08306b","#253494", "#081d58")
}
if(palette=="cream_pink"){
colors <- c("#feebe2","#fcc5c0","#fa9fb5","#f768a1", "#c51b8a","#7a0177", "#49006a")
}
if(palette=="cream_green"){
colors <- c("#f7fcb9","#d9f0a3","#addd8e","#78c679","#41ab5d", "#238443","#006837", "#006837")
}
if(palette=="cream_brown"){
colors <- c("#ffeda0","#fed976","#feb24c","#fd8d3c","#fc4e2a","#e31a1c","#bd0026","#800026")
}
split_factor <- round((max_key/6),1)
breaks = seq(min_key,max_key, by = split_factor)
if(top_line==TRUE){
top_annotation = ComplexHeatmap::HeatmapAnnotation(
lines = EnrichedHeatmap::anno_enriched(gp = grid::gpar(fontsize=12),
ylim=c(ymin, ymax),
# yaxis_side = "right",
# yaxis_facing = "inside",
# yaxis_gp = grid::gpar(fontsize = 10, lwd=1.5)),
axis_param = list(gp=grid::gpar(fontsize = 10, lwd=1.5),
side="right",
facing="inside")))
}
else{
top_annotation = NULL
}
get_enrichment_heatmap_list <- function(x, names, titles, ...) {
ll <- length(x)
## first heatmap
ehml <- EnrichedHeatmap::EnrichedHeatmap(mat = eval(mat1), name = names[[1]], column_title = titles[[1]], show_heatmap_legend = T,
row_order= eval(row_order),
use_raster = TRUE, ...)
## several other heatmaps if length of x > 1.
if (ll > 1) {
for (i in 2:ll) {
print(i)
ehml <- ehml +
EnrichedHeatmap::EnrichedHeatmap(
mat = eval(mat2),
row_order= eval(row_order),
name = ifelse(length(names) >= i, names[i], "NA"),
use_raster = TRUE,
column_title = ifelse(length(titles) >= i, titles[i], "NA"),
show_heatmap_legend = ifelse(length(names) >= i, TRUE, FALSE), ...
) ## legend will be shown only if the name is given for a heatmap.
}
}
return(ehml)
}
ehm_list <- get_enrichment_heatmap_list(x = xx$norm_matrix,names = xx$names,titles = xx$names,
cluster_rows = FALSE,
pos_line = TRUE,
show_row_names = FALSE,
axis_name_rot = 90,
heatmap_legend_param = list(color_bar = "continuous",legend_direction="horizontal"),
axis_name = axis_name,
col = circlize::colorRamp2(breaks = breaks,
colors = as.vector(colors[1:length(breaks)])),
top_annotation = top_annotation
)
if(output==TRUE){
print("plotting")
png(file=paste(output_name, length(sub_feature_gr),"hm.png", sep="_"),width=nrow(bw_files)*1.5,height=5.5,pointsize = 8, res=300,units = "in")
ComplexHeatmap::draw(ehm_list, heatmap_legend_side = "top", gap = grid::unit(1.5, "mm"))
dev.off()
}
else{
return(ComplexHeatmap::draw(ehm_list, heatmap_legend_side = "right", gap = grid::unit(1.5, "mm")))
}
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.