Writting eigenvectors table to analyze their functions

seoanezonjic · Jul 17, 2024 · e93da28 · e93da28
1 parent 25e3c0b
commit e93da28
Show file tree

Hide file tree

Showing 3 changed files with 41 additions and 16 deletions.
diff --git a/R/main_degenes_Hunter.R b/R/main_degenes_Hunter.R
@@ -172,7 +172,7 @@ main_degenes_Hunter <- function(
     ############################################################
     ##             PERFORM EXPRESION ANALYSIS                 ##
     ############################################################
-    dir.create(output_files)
+    check_and_create_dir(output_files)
 
 
     exp_results <- perform_expression_analysis(modules, replicatesC, 
@@ -217,7 +217,7 @@ main_degenes_Hunter <- function(
     if(grepl("X", modules)) { # CASE X: diffcoexp
       cat('Correlation analysis is performed with diffcoexp\n')
       path <- file.path(output_files, "Results_diffcoexp")
-      dir.create(path)
+      check_and_create_dir(path)
 
       results_diffcoexp <- analysis_diffcoexp(data = raw_filter,
                                            path = path,
@@ -315,6 +315,8 @@ main_degenes_Hunter <- function(
     final_results[["numeric_factors"]] <- numeric_factors
     final_results[["string_factors"]] <- string_factors
     final_results[["PCA_res"]] <- PCA_res
+    final_results[["library_sizes"]] <- library_sizes
+
     if(!is.null(combinations_WGCNA)){
       final_results <- c(final_results, combinations_WGCNA)
     }

diff --git a/R/main_functional_hunter.R b/R/main_functional_hunter.R
@@ -141,6 +141,14 @@ main_functional_hunter <- function(
         }
     }
 
+
+    # prepare PCA data for KS
+    pca_eigenvectors <- list(pca_all_genes  = hunter_results$pca_all_genes,
+                             pca_degs = hunter_results$pca_degs)
+
+    pca_eigenvectors <- lapply(pca_eigenvectors, name_all_columns)
+
+    pca_eigenvectors <- lapply(pca_eigenvectors, parse_eigenvectors)
     ############################################################
     ##                                                        ##
     ##                     PERFORM ENRICHMENTS                ## 
@@ -159,6 +167,7 @@ main_functional_hunter <- function(
     }
 
     if("ORA" %in% enrich_methods){
+
         deg_enr_ora  <- multienricher_ora(all_funsys = enrich_dbs, 
             genes_list = prev_genes, organism_info = current_organism_info, 
             pvalueCutoff = pthreshold, qvalueCutoff = qthreshold, 
@@ -184,6 +193,24 @@ main_functional_hunter <- function(
         }
     }
 
+    if ("topGO" %in% enrich_methods){
+        func_results$topGO <- multienricher_topGO(all_funsys=enrich_dbs, 
+                                                   universe=universe,
+                                                   genes_list=prev_genes,
+                                                   organism_info=current_organism_info,
+                                                   p_value_cutoff = pthreshold)
+    }
+
+    func_results$PCA_enrichments <- lapply(pca_eigenvectors, 
+                                           multienricher_topGO, 
+                                           all_funsys = enrich_dbs,
+                                           organism_info = current_organism_info,
+                                           p_value_cutoff = pthreshold,
+                                           algorithm = "classic", 
+                                           statistic = "ks",
+                                           gene_id = "ensembl", 
+                                           scoreOrder = "decreasing")
+
     ############################################################
     ##                                                        ##
     ##                     EXPORT DATA                        ## 

diff --git a/R/write_report.R b/R/write_report.R
@@ -52,6 +52,7 @@ write_expression_report <- function(exp_results,
     numeric_factors <- exp_results[["numeric_factors"]]
     string_factors <- exp_results[["string_factors"]] 
     PCA_res <- exp_results[["PCA_res"]]  
+    library_sizes <- exp_results[["library_sizes"]]
     outf <- file.path(normalizePath(output_files),"DEG_report.html")
     rmarkdown::render(file.path(template_folder, 'main_report.Rmd'),
                       output_file = outf, intermediates_dir = output_files)
@@ -78,24 +79,17 @@ write_expression_data <- function(final_results, output_files){
 }
 
 write_pca_data <- function(PCA_res, output_files){
+
     pca_output <- file.path(output_files, "PCA_results")
     dir.create(pca_output)
-    all_genes_pca <- PCA_res$all_genes$dim_data_merged
+    all_genes_pca <- PCA_res$all_genes
 
-    all_genes_merged_metrics <- merge_dim_table_metrics(all_genes_pca)
+    write_general_pca(all_genes_pca, pca_output, "all_genes_")
+
 
-    write.table(all_genes_merged_metrics, file = file.path(pca_output, "all_genes_dim_metrics.txt"),sep = "\t", quote = FALSE, row.names=FALSE)
-    hcpc_table <- PCA_res$all_genes$res.hcpc$call$X
-    hcpc_table$samples <- rownames(hcpc_table)
-    write.table(hcpc_table, file = file.path(pca_output, "all_genes_hcpc.txt"),sep = "\t", quote = FALSE, row.names=FALSE)
-    prevalent_pca <- PCA_res$DEGs$dim_data_merged
+    prevalent_pca <- PCA_res$DEGs
     if (!is.null(prevalent_pca)){
-        hcpc_table_deg <- PCA_res$DEGs$res.hcpc$call$X
-        hcpc_table_deg$samples <- rownames(hcpc_table_deg)
-        write.table(hcpc_table_deg, file = file.path(pca_output, "prevalent_hcpc.txt"),sep = "\t", quote = FALSE, row.names=FALSE)
-        prevalent_merged_metrics <- merge_dim_table_metrics(prevalent_pca)
-        write.table(prevalent_merged_metrics, file = file.path(pca_output, "prevalent_dim_metrics.txt"),sep = "\t", quote = FALSE, row.names=FALSE)
-
+       write_general_pca(prevalent_pca, pca_output, "prevalent_")
     }
 
 }
@@ -104,7 +98,9 @@ write_general_pca <- function(pca_data, output_files, tag = ""){
   pca_res <- pca_data$dim_data_merged
 
   merged_metrics <- merge_dim_table_metrics(pca_res)
-
+
+  dimnames(pca_data$pca_data$svd$V) <- dimnames(pca_data$pca_data$var$cor)
+  write.table(pca_data$pca_data$svd$V, quote = FALSE, file = file.path(output_files, paste0(tag, "eigenvectors.txt")), sep = "\t")
   write.table(merged_metrics, file = file.path(output_files, paste0(tag, "dim_metrics.txt")),sep = "\t", quote = FALSE, row.names=FALSE)
   hcpc_table <- pca_data$res.hcpc$call$X
   hcpc_table$samples <- rownames(hcpc_table)