From fe7edee67d4783dcbfb1e7bab80a53f2ac710ea4 Mon Sep 17 00:00:00 2001 From: JoseCorCab Date: Fri, 13 Sep 2024 08:51:57 +0200 Subject: [PATCH] Enrichment for HCPC clusters inplemented --- R/factor_mining.R | 21 +++ R/main_functional_hunter.R | 15 +- inst/templates/functional_report.Rmd | 202 ++++++++++++++------------- inst/templates/pca_enr.Rmd | 17 ++- 4 files changed, 153 insertions(+), 102 deletions(-) diff --git a/R/factor_mining.R b/R/factor_mining.R index 2dab0db..dad63fd 100644 --- a/R/factor_mining.R +++ b/R/factor_mining.R @@ -146,6 +146,27 @@ get_and_parse_pca_eigenvectors <- function(pca_res, cor_pval_cutoff = 0.05, eig } +get_and_parse_clusters <- function(pca_res){ + clusters <- pca_res$res.hcpc$call$X + dims <- colnames(clusters)[grepl("Dim.", colnames(clusters))] + weights <- clusters |> dplyr::group_by(clust) |> dplyr::summarise_at(dims, mean) + dimnames(pca_res$pca_data$svd$V) <- dimnames(pca_res$pca_data$var$cor) + clust_names <- paste("Cluster",weights$clust ) + weights <- as.data.frame(dplyr::select(weights, -clust)) + rownames(weights) <- clust_names + weighted_eigen <- apply(weights,1, get_clust_contributions, + eigenvectors = pca_res$pca_data$svd$V) + colnames(weighted_eigen) <- rownames(weights) + weighted_eigen <- as.data.frame(weighted_eigen) + weighted_eigen <- lapply(weighted_eigen, function(col) setNames(col, rownames(weighted_eigen))) + return(weighted_eigen) +} + +get_clust_contributions <- function(weigths, eigenvectors) { + weighted_eig <- weigths * eigenvectors + rowSums(weighted_eig) +} + filter_eigenvectors <- function(eigenvectors, cor_sig, pval_cutoff) { eigenvec_fil <- lapply(names(eigenvectors), function(dimension) { diff --git a/R/main_functional_hunter.R b/R/main_functional_hunter.R index 4a9c40d..5556e8a 100644 --- a/R/main_functional_hunter.R +++ b/R/main_functional_hunter.R @@ -142,11 +142,13 @@ main_functional_hunter <- function( } if (!is.null(hunter_results$pca_data)) { - # prepare PCA data for KS + # prepare PCA data for topGO + pca_eigenvectors <- lapply(hunter_results$pca_data, get_and_parse_pca_eigenvectors, cor_pval = 1, eig_abs_cutoff = NULL) + pca_clusters <- lapply(hunter_results$pca_data, get_and_parse_clusters) } @@ -215,6 +217,17 @@ main_functional_hunter <- function( scoreOrder = "decreasing", workers = cores, task_size = task_size) + func_results$PCA_clusters_results <- lapply(pca_clusters, + multienricher_topGO, + all_funsys = enrich_dbs, + organism_info = current_organism_info, + p_value_cutoff = pthreshold, + algorithm = "classic", + statistic = "t", + gene_id = "ensembl", + scoreOrder = "decreasing", + workers = cores, + task_size = task_size) } ############################################################ ## ## diff --git a/inst/templates/functional_report.Rmd b/inst/templates/functional_report.Rmd index 3008f19..ee0ec96 100644 --- a/inst/templates/functional_report.Rmd +++ b/inst/templates/functional_report.Rmd @@ -1,101 +1,103 @@ ---- -title: "ExpHunterSuite: Functional Report" -author: "SysBioLab" -output: - html_document: - toc: true - toc_float: true - df_print: paged - fig_width: 12 ---- - - - -```{r config, include=FALSE, message=FALSE} -require(ggplot2) -require(knitr) -require(clusterProfiler) -require(enrichplot) -require(DOSE) -``` - -```{r child="func_initial_details.Rmd"} -``` -```{r child="func_deg_details.Rmd"} -``` -```{r child="func_top_genes.Rmd"} -``` - - -```{r ORA_analysis, echo=FALSE, results='asis', message=FALSE, warning=FALSE} - - -res_ora <- list() -for(funsys in names(flags_ora)) { - if(flags_ora[[funsys]] == TRUE) { - enrich_obj <- func_results$ORA[[funsys]] - exp_res <- knitr::knit_expand("ora_plots_all_objects.Rmd") - res_ora[[funsys]] <- knitr::knit(text=exp_res, quiet=TRUE) - } else { - res_ora[[funsys]] <- paste0("## **No ORA enrichment found for ", funsys, "**\n") - } -} -cat(unlist(res_ora), sep = '\n') - - -res_gsea <- list() -for(funsys in names(flags_gsea)) { - if(flags_gsea[[funsys]] == TRUE) { - enrich_obj <- func_results$GSEA[[funsys]] - exp_res <- knitr::knit_expand("gsea_plots_all_objects.Rmd") - res_gsea[[funsys]] <- knitr::knit(text=exp_res, quiet=TRUE) - } else { - res_gsea[[funsys]] <- paste0("## **No GSEA enrichment found for ", funsys, "**\n") - } -} -cat(unlist(res_gsea), sep = '\n') - -``` - -```{r PCA_enr, echo=FALSE, results='asis', message=FALSE, warning=FALSE} -res_PCA <- list() - -for(funsys in names(flags_ora)) { - if (funsys %in% c("BP", "MF", "CC")){ - pca_enr_all <- func_results$PCA_enrichments$all_genes[[funsys]] - pca_enr_degs <- func_results$PCA_enrichments$DEGs[[funsys]] - pca_enr <- knitr::knit_expand("pca_enr.Rmd") - res_PCA[[funsys]] <- knitr::knit(text=pca_enr, quiet=TRUE) - } -} - -cat(unlist(res_PCA), sep = '\n') -``` - -## **Values of options passed to the Functional Hunter main function** -First column contains the option names; second column contains the given values for each option. -Note that large data objects (e.g. expression results, - organism table and custom annotation files) are not shown. -```{r opt_vals, echo = FALSE, results='asis'} -final_main_params <- func_results$final_main_params[! sapply(func_results$final_main_params, is.list)] -print(knitr::kable(cbind(final_main_params))) +--- +title: "ExpHunterSuite: Functional Report" +author: "SysBioLab" +output: + html_document: + toc: true + toc_float: true + df_print: paged + fig_width: 12 +--- + + + +```{r config, include=FALSE, message=FALSE} +require(ggplot2) +require(knitr) +require(clusterProfiler) +require(enrichplot) +require(DOSE) +``` + +```{r child="func_initial_details.Rmd"} +``` +```{r child="func_deg_details.Rmd"} +``` +```{r child="func_top_genes.Rmd"} +``` + + +```{r ORA_analysis, echo=FALSE, results='asis', message=FALSE, warning=FALSE} + + +res_ora <- list() +for(funsys in names(flags_ora)) { + if(flags_ora[[funsys]] == TRUE) { + enrich_obj <- func_results$ORA[[funsys]] + exp_res <- knitr::knit_expand("ora_plots_all_objects.Rmd") + res_ora[[funsys]] <- knitr::knit(text=exp_res, quiet=TRUE) + } else { + res_ora[[funsys]] <- paste0("## **No ORA enrichment found for ", funsys, "**\n") + } +} +cat(unlist(res_ora), sep = '\n') + + +res_gsea <- list() +for(funsys in names(flags_gsea)) { + if(flags_gsea[[funsys]] == TRUE) { + enrich_obj <- func_results$GSEA[[funsys]] + exp_res <- knitr::knit_expand("gsea_plots_all_objects.Rmd") + res_gsea[[funsys]] <- knitr::knit(text=exp_res, quiet=TRUE) + } else { + res_gsea[[funsys]] <- paste0("## **No GSEA enrichment found for ", funsys, "**\n") + } +} +cat(unlist(res_gsea), sep = '\n') + +``` + +```{r PCA_enr, echo=FALSE, results='asis', message=FALSE, warning=FALSE} +res_PCA <- list() + +for(funsys in names(flags_ora)) { + if (funsys %in% c("BP", "MF", "CC")){ + pca_enr_all <- func_results$PCA_enrichments$all_genes[[funsys]] + pca_enr_degs <- func_results$PCA_enrichments$DEGs[[funsys]] + pca_cl_all <- func_results$PCA_clusters_results$all_genes[[funsys]] + pca_cl_degs <- func_results$PCA_clusters_results$DEGs[[funsys]] + pca_enr <- knitr::knit_expand("pca_enr.Rmd") + res_PCA[[funsys]] <- knitr::knit(text=pca_enr, quiet=TRUE) + } +} + +cat(unlist(res_PCA), sep = '\n') +``` + +## **Values of options passed to the Functional Hunter main function** +First column contains the option names; second column contains the given values for each option. +Note that large data objects (e.g. expression results, + organism table and custom annotation files) are not shown. +```{r opt_vals, echo = FALSE, results='asis'} +final_main_params <- func_results$final_main_params[! sapply(func_results$final_main_params, is.list)] +print(knitr::kable(cbind(final_main_params))) ``` \ No newline at end of file diff --git a/inst/templates/pca_enr.Rmd b/inst/templates/pca_enr.Rmd index 6d5fd9a..f865936 100644 --- a/inst/templates/pca_enr.Rmd +++ b/inst/templates/pca_enr.Rmd @@ -42,10 +42,17 @@ while (last_dim < func_results$pca_data$all_genes$dim_to_keep){ print_go_tables(pca_enr_all) +cat("### **Functional enrichments of HCPC clusters using all genes**\n") +plot(func_results$pca_data$all_genes$res.hcpc, choice="tree", title="") +plot(func_results$pca_data$all_genes$res.hcpc, axes=c(1,2), choice="map", draw.tree= FALSE) + +print_go_tables(pca_cl_all) + + cat("### **Functional enrichments of PCs using only DEGs**\n") last_dim <- 1 -while (last_dim < func_results$pca_data$all_genes$dim_to_keep){ +while (last_dim < func_results$pca_data$DEGs$dim_to_keep){ plot(FactoMineR::plot.PCA(func_results$pca_data$DEGs$pca_data, invisible = "quali", axis = c(last_dim, last_dim + 1), @@ -55,5 +62,13 @@ while (last_dim < func_results$pca_data$all_genes$dim_to_keep){ print_go_tables(pca_enr_degs) +cat("### **Functional enrichments of HCPC clusters using only DEGs**\n") + +plot(func_results$pca_data$DEGs$res.hcpc, choice="tree", title="") +plot(func_results$pca_data$DEGs$res.hcpc, axes=c(1,2), choice="map", draw.tree= FALSE) + +print_go_tables(pca_cl_degs) + + ``` \ No newline at end of file