R/DEseq2_export.R

In Rvirgenslane/Hotgenes: For Visualizing Differential Expression Analysis

#' Prepares DESeq2 object for visualization
#' @importFrom grDevices dev.off recordPlot
#' @import dplyr
#' @importFrom grDevices colorRampPalette pdf
#' @importFrom graphics boxplot
#' @importFrom stats coef cutree setNames
#' @importFrom utils head
#' @importFrom rlang .data
#' @export
#' @import DESeq2
#' @import ggplot2
#' @import SummarizedExperiment
#' @param DEseq2_object the output object generated by DESeq2
#' @param ... additional parameters for lfcShrink function
#' @param lfcShrink_type "normal" (Default) is the original DESeq2 shrinkage
#' estimator. Other options is "apeglm".
#' @param resultsSel values from DESeq2 resultsNames function. If NULL
#' (default) DE for all results will be reported.
#' @param padj_cut the adjusted pvalue cut off for differential
#' expression.  Default = 0.1
#' @param rlog_export Logical value, if TRUE (default), data
#' normalized by rlog will be exported. If FALSE, only 
#' varianceStabilizingTransformation normalization will be used.
#' @seealso \code{\link[DESeq2]{lfcShrink}}
#' @return list of R objects needed for Shiny_DE_viz function
#' @examples
#' dds_con_dir<-system.file("extdata",
#' "dds_con.Rdata", package = "Hotgenes", mustWork = TRUE)
#' load(dds_con_dir)
#' library(DESeq2)
#' resultsNames(dds_con)
#' Example_Hotgenes<-DEseq2_export(DEseq2_object = dds_con, 
#' resultsSel = contrast_selec<-c("shEWS","Hrs2","Hrs6",
#' "shEWS.Hrs2","shEWS.Hrs6"), padj_cut = 0.1)
#' Example_Hotgenes$Exported_plots$transformation_plots
#' Example_Hotgenes$Exported_plots$boxplot_rlog
#' Example_Hotgenes$Exported_plots$boxplot_vsd

DEseq2_export<-function(DEseq2_object = NULL,
lfcShrink_type ="normal",
resultsSel=NULL,
rlog_export=TRUE,
padj_cut=0.1,...){

# setting up global variables

# extracting design data
design_data<-data.frame((colData(DEseq2_object)))
design_data$sizeFactor<-NULL

# betas
betas_con <- data.frame(coef(DEseq2_object), check.names = FALSE)

# Output lists

if(!is.null(resultsSel)){
    all_contrasts<-resultsNames(DEseq2_object)[-c(1)]
    
    all_contrasts<-all_contrasts[all_contrasts %in% resultsSel]
}else if(is.null(resultsSel)){
    all_contrasts<-resultsNames(DEseq2_object)[-c(1)]
}



Output_lists <- setNames(vector(length(all_contrasts),
mode="list"), all_contrasts)

Output_DE <- Output_lists

# Enrichment/Depletion lists
Enriched_by <- Output_lists
Depleted_by <- Output_lists


for (i in all_contrasts) {
res <- lfcShrink(DEseq2_object, coef = i, type = lfcShrink_type, quiet = FALSE,...)

# Filtering readout based on adjusted pvalue
cutL<-sum(res$padj < padj_cut, na.rm=TRUE)

# Prepares a table for the significant genes
res_con<-data.frame(res)
# head(res_con)
Sig_DE<-res_con[order(res_con$padj),]
Sig_DE<-Sig_DE[seq(0,cutL),]

if(lfcShrink_type == "apeglm"){
    Temp_Sig_DE<-Sig_DE
    Temp_Sig_DE$logP<--log10(Sig_DE$pvalue)
    Sig_DE$stat <- Temp_Sig_DE$logP/sign(Temp_Sig_DE$log2FoldChange)
}


if(cutL!=0){
Output_DE[[i]]<-Sig_DE

Enriched_by[[i]]<-rownames(Sig_DE[Sig_DE$log2FoldChange > 0,])
Depleted_by[[i]]<-rownames(Sig_DE[Sig_DE$log2FoldChange < 0,])
Output_lists[[i]]<-rownames(Sig_DE)
}
}

# stabilizing variance for plots and expression matrix --------------------

vsd <- varianceStabilizingTransformation(DEseq2_object, blind = FALSE)
df_vsd<-data.frame(assay(vsd), check.names = FALSE)

# plot selection ----------------------------------------------------------
log2_plus_1_norm<-log2(counts(DEseq2_object, normalized=TRUE)[, seq(1,2)]+1)

# boxplot(df_rld)
boxplot(df_vsd)
boxplot_vsd <- recordPlot()
dev.off() ## clean up device


# rlog_export
if(rlog_export==TRUE){
rld <- rlog(DEseq2_object, blind = FALSE)
df_rld<-data.frame(assay(rld), check.names = FALSE)

Normalized_Expression<-list(vsd=df_vsd,
rld=df_rld)

# boxplot(df_rld)
boxplot(df_rld)
boxplot_rlog <- recordPlot()
dev.off() ## clean up device




df <- bind_rows(
    as_tibble(log2_plus_1_norm) %>% mutate(transformation = "log2(x + 1)"),
    as_tibble(df_rld[, seq(1,2)]) %>% mutate(transformation = "rlog"),
    as_tibble(df_vsd[, seq(1,2)]) %>% mutate(transformation = "vst"))
colnames(df)[seq(1,2)] <- c("x", "y")

# transformation_plots
transformation_plots <- ggplot(df, aes(x = .data$x, y = .data$y)) +
    geom_hex(bins = 80) +
    coord_fixed() + facet_grid( . ~ transformation)

Exported_plots<-list(transformation_plots=transformation_plots,
                     boxplot_rlog=boxplot_rlog,
                     boxplot_vsd=boxplot_vsd)

}else if(rlog_export!=TRUE){
    Normalized_Expression<-list(vsd=df_vsd)
    
    df <- bind_rows(
        as_tibble(log2_plus_1_norm) %>% mutate(transformation = "log2(x + 1)"),
        as_tibble(df_vsd[, seq(1,2)]) %>% mutate(transformation = "vst"))
    colnames(df)[seq(1,2)] <- c("x", "y")
    
    # transformation_plots
    transformation_plots <- ggplot(df, aes(x = .data$x, y = .data$y)) +
        geom_hex(bins = 80) +
        coord_fixed() + facet_grid( . ~ transformation)
    
    
    Exported_plots<-list(transformation_plots=transformation_plots,
                         boxplot_vsd=boxplot_vsd)
}



# making query_dataset ----------------------------------------------------
# rlog_export
if(rlog_export==TRUE){
Temp_data<-data.frame(t(df_rld))

}else if(rlog_export!= TRUE){
Temp_data<-data.frame(t(df_vsd))
}

# Merging
Final_DF<-merge(design_data,
Temp_data,
by = "row.names")
rownames(Final_DF)<-Final_DF$Row.names
Final_DF$Row.names<-NULL

#
Export<-list(Output_DE=Output_DE,
Output_lists=Output_lists,
Enriched_by=Enriched_by,
Depleted_by=Depleted_by,
betas_con=betas_con,
Exported_plots=Exported_plots,
Normalized_Expression=Normalized_Expression,
design_data=design_data,
Query_set=Final_DF)
return(Export)


}