R/Limma_export.R

In Rvirgenslane/Hotgenes: For Visualizing Differential Expression Analysis

#' Prepares Limma object for visualization
#' @importFrom limma decideTests
#' @export
#' @param limmafit the fit generated by limma
#' @param design_data data.frame containing Experimental design
#' @param Expression_dat R object containing the expression data
#' used for limma differential expression analysis.
#' @inheritParams limma::decideTests
#' @return list of R objects needed for Shiny_DE_viz function
#' @examples
#' Example_Hotgenes_dir<-system.file("extdata",
#' "Example_Hotgenes.Rdata",
#' package = "Hotgenes", mustWork = TRUE)
#' load(Example_Hotgenes_dir)
#' 
#' library(limma)
#' 
#' exp<-Example_Hotgenes$Normalized_Expression$rld
#' design_m<-Example_Hotgenes$design_data
#' 
#' design_matrix <- model.matrix(~sh*Hrs+Bio_Rep,   
#' data = design_m)
#' 
#' aw <- arrayWeights(exp, design_matrix)
#' 
#' fit <- lmFit(exp, design=design_matrix, weights = aw)
#' fit <- eBayes(fit, robust = TRUE) 
#' 
#' L_out<-Limma_export(Expression_dat = exp, 
#' design_data = design_m, 
#' limmafit = fit)
#' summary(L_out$Output_lists)


Limma_export<-function(limmafit = NULL,
design_data=NULL,
Expression_dat = NULL,
method="global",
adjust.method="BH",
p.value=0.1){

# verify that design and expression data align
if(!all(rownames(design_data) == colnames(Expression_dat))){
stop("design_data does not match Expression_dat")
}

# verify that rownames(limmafit$design) and expression data align 
  if(!all(rownames(limmafit$design) == colnames(Expression_dat))){
    stop("limmafit$design does not match Expression_dat")
  }  

# setting up global variables
res <- decideTests(limmafit, method=method, 
adjust.method=adjust.method,
p.value=p.value)

df_results<-data.frame(res, 
stringsAsFactors = FALSE,
check.names = FALSE)

# Getting adj pvalues
Adj_res <- Limma_Adj_pval(limmafit, method=method, 
adjust.method=adjust.method,
p.value=p.value)

df_Adj_res<-data.frame(Adj_res, 
stringsAsFactors = FALSE,
check.names = FALSE)

# betas
betas_con <-data.frame(limmafit[["coefficients"]], check.names = FALSE)

# Output lists
all_contrasts<-colnames(limmafit)[-1]

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

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

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


# Getting output
for (i in all_contrasts) {

sig_genes<-rownames(df_results[df_results[i] != 0 &
!is.na(df_results[i]),][i])

Output_lists[[i]]<-sig_genes
Output_DE[[i]]<-data.frame(Genes=sig_genes,
log2FoldChange=limmafit[["coefficients"]][sig_genes,i],  
stat=limmafit[["t"]][sig_genes,i],   
pvalue=limmafit[["p.value"]][sig_genes,i], 
padj=df_Adj_res[sig_genes,i],
stringsAsFactors = FALSE)



Output_DE[[i]]<-Output_DE[[i]][order(Output_DE[[i]]$padj, decreasing = FALSE),]
rownames(Output_DE[[i]])<-Output_DE[[i]]$Genes
Enriched_by[[i]]<-Output_DE[[i]][Output_DE[[i]]$log2FoldChange > 0,]$Genes
Depleted_by[[i]]<-Output_DE[[i]][Output_DE[[i]]$log2FoldChange < 0,]$Genes
Output_DE[[i]]$Genes<-NULL
}

# fixing names
names(Output_DE)<-make.names(names(Output_DE), unique = FALSE, allow_ = TRUE)
names(Output_lists)<-make.names(names(Output_lists), unique = FALSE, allow_ = TRUE)
names(Enriched_by)<-make.names(names(Enriched_by), unique = FALSE, allow_ = TRUE)
names(Depleted_by)<-make.names(names(Depleted_by), unique = FALSE, allow_ = TRUE)
names(betas_con)<-make.names(names(betas_con), unique = FALSE, allow_ = TRUE)


# expression matrix -------------------------------------------------------
Normalized_data<-data.frame(Expression_dat, check.names = FALSE)
Normalized_Expression<-list(Normalized_data=Normalized_data)

# transformation_plots

#dev.off() ## clean up device

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


Exported_plots<-list(transformation_plots=NULL,
boxplot_Normalized_data=boxplot_Normalized_data)


# Getting expression data

Temp_data<-data.frame(t(Normalized_data))


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

# Final list
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)
}