updates Jan 2024

eleonore-schneeg · Jan 12, 2024 · 7351c06 · 7351c06
1 parent b84b55c
commit 7351c06
Show file tree

Hide file tree

Showing 51 changed files with 9,692 additions and 74 deletions.
diff --git a/NAMESPACE b/NAMESPACE
@@ -49,6 +49,7 @@ export(multiomics_modules)
 export(multiomics_network)
 export(multiomics_network_matrix)
 export(pathway_analysis_enrichr)
+export(pathway_report)
 export(plot_OpenTarget)
 export(plot_communities)
 export(plot_components)
@@ -61,13 +62,15 @@ export(process_rna)
 export(processed_proteomics)
 export(processed_rna)
 export(protein_DE_analysis)
+export(proteomics_report)
 export(pseudotime_inference)
 export(qc_report)
 export(rna_DE_analysis)
 export(sign_signature)
 export(single_omic_comparisons)
 export(single_omic_report)
 export(vertical_integration)
+export(visualise_pathway)
 export(volcano_interactive)
 export(volcano_interactive_comparison)
 export(volcano_plot_deseq)
@@ -79,6 +82,7 @@ import(circlize)
 import(dplyr)
 import(ggplot2)
 import(igraph)
+import(visNetwork)
 importFrom(ActivePathways,read.GMT)
 importFrom(AnnotationDbi,mapIds)
 importFrom(DESeq2,DESeq)

diff --git a/R/TF_enrichment.R b/R/TF_enrichment.R
@@ -19,12 +19,19 @@
 Transcription_Factor_enrichment <- function(communities,
                                             weights,
                                             TF_gmt = "Enrichr_Queries.gmt",
-                                            threshold = 0.2,
+                                            threshold = 1.5,
                                             direction = "negative") {
   TF_genelist <- ActivePathways::read.GMT(TF_gmt)
   TF_genelist_list <- lapply(TF_genelist, function(x) {
     x$genes
   })
+
+  mean_weight <- mean(weights$weights_df$protein$Weights)
+  sd_weight <- sd(weights$weights_df$protein$Weights)
+
+  # Identify rows where Weight is over 2 SD above the mean
+  threshold <- mean_weight + threshold * sd_weight
+
 
   if (direction == "negative") {
     TF_targets <- names(TF_genelist_list)[names(TF_genelist_list) %in% weights$weights_df$protein$Feature[which(weights$weights_df$protein$Weights <= -threshold)]]

diff --git a/R/extract_weights.R b/R/extract_weights.R
@@ -4,7 +4,7 @@
 #' @param model MOFA model from integration in
 #' `multiassay@metadata$integration$MOFA`
 #' @param factor Factor to extract weight from
-#' @param threshold Absolute threshold to filter weights
+#' @param threshold Number of standard deviation away from the mean weight
 #' @param sense_check_variable default to `NULL`. High weights should
 #' coincide with stronger correlation if the sense_check_variable is an
 #' important driver of variation in the designated factor. Will be used to
@@ -26,7 +26,7 @@
 #'
 extract_weigths <- function(model,
                             factor = 1,
-                            threshold = 0.3,
+                            threshold = 1.5,
                             sense_check_variable = NULL) {
   weights_rna_1 <- MOFA2::get_weights(model,
     view = "mRNA",
@@ -36,7 +36,6 @@ extract_weigths <- function(model,
     as.data.frame = FALSE
   )
 
-
   weights_prot_1 <- MOFA2::get_weights(model,
     view = "proteins",
     factor = factor,
@@ -56,6 +55,12 @@ extract_weigths <- function(model,
 
   rna_1$model_feature <- rownames(rna_1)
   protein_1$model_feature <- rownames(protein_1)
+
+  mean_weight <- mean(rna_1$Weights)
+  sd_weight <- sd(rna_1$Weights)
+
+  # Identify rows where Weight is over 2 SD above the mean
+  threshold <- mean_weight + threshold * sd_weight
 
   protein_positive <- protein_1 %>%
     dplyr::arrange(dplyr::desc(Weights), dplyr::desc(Feature)) %>%

diff --git a/R/integrative_results_unsupervised.R b/R/integrative_results_unsupervised.R
@@ -1,7 +1,7 @@
 #' Downstream analysis for integrated multi-omics data
 #'
 #' @param multiassay Multiassay experiment object generated by Omix
-#' @param integration_model Possible unsupervised integration methods are `MOFA`,
+#' @param integration_model Possible unsupervised integration methods are `MOFA`
 #' @param dependent Dependent variable for groups.
 #' @param correlation_threshold Absolute correlation threshold to draw edges
 #' in the network
@@ -11,7 +11,7 @@
 #' @param sense_check_variable sense check
 #' @param covariates cov
 #' @param time tim
-#' @param weights_threshold weights threshold
+#' @param weights_threshold weights threshold (+/- SD from the mean)
 #' @param community_detection community detection method
 #' @param pseudotime_analysis TRUE/ FALSE
 #' @param TF_fp file path GMT file with Transcription Factors and target genes. Check `https://maayanlab.cloud/chea3/`
@@ -33,7 +33,7 @@ integrative_results_unsupervised <- function(multiassay,
                                              dependent = "diagnosis",
                                              sense_check_variable = "PHF1",
                                              covariates = c("PHF1", "amyloid", "AT8"),
-                                             weights_threshold = 0.2,
+                                             weights_threshold = 1.5,
                                              correlation_threshold = 0.4,
                                              community_detection = "leading_eigen",
                                              pseudotime_analysis = TRUE,

diff --git a/R/pathway_report.R b/R/pathway_report.R
@@ -0,0 +1,62 @@
+#' Generate a report for single omic analysis
+#'
+#' @param pathways list of pathways 
+#' @param report_folder_path report_folder_path folder path to save the report.
+#' @param report_file  filename for report (without an extension).
+#' @param database Enrichment database `GO_Molecular_Function_2021`,`GO_Cellular_Component_2021`,
+#'  `GO_Biological_Process_2021`, `Reactome_2016` , `KEGG_2021_Human` , `MSigDB_Hallmark_2020`
+#' @family Report
+#' @return Single omic analyses report html
+#' @export
+#'
+pathway_report <- function(pathways,
+                              report_folder_path = getwd(),
+                              report_file = "single_omic_report_Omix",
+                              database='Reactome_2016',
+                              num_path=20){
+  uniomic=list()
+  uniomic$param$pathways=pathways
+  uniomic$param$database=database
+  uniomic$param$num_path=num_path
+
+
+  report_file <- tools::file_path_sans_ext(report_file)
+
+  cli::cli_h2("Generating report for pathways analyses")
+
+
+  metadata_tmp_path <- file.path(tempdir(), "metadata.qs")
+
+  cli::cli_text("Writing temp files for report...")
+  qs::qsave(
+    uniomic,
+    metadata_tmp_path
+  )
+
+  krd <- file.path(tempdir(), "krdqc")
+  intd <- file.path(tempdir(), "idqc")
+  dir.create(krd, showWarnings = FALSE)
+  dir.create(intd, showWarnings = FALSE)
+
+  cli::cli_text("Generating Pathway report...")
+  rmarkdown::render(
+    system.file(
+      "rmarkdown/templates/pathways/skeleton.Rmd",
+      package = "Omix"
+    ),
+    params = list(
+      metadata_path = metadata_tmp_path
+    ),
+    output_dir = report_folder_path,
+    output_file = report_file,
+    knit_root_dir = krd,
+    intermediates_dir = intd,
+    quiet = TRUE
+  )
+
+  report_file_name <- paste(report_file, ".html", sep = "")
+
+
+
+}
+
diff --git a/R/process_protein.R b/R/process_protein.R
@@ -228,7 +228,7 @@ process_protein <- function(multiassay,
       outliers
 
     processing_outputs[["remove_sample_outliers"]]=list(sample_outliers= dim5 - dim6,
-                                                         percentage_outliers= (dim5 - dim6)/dim4,
+                                                         percentage_outliers= (dim5 - dim6)/dim5,
                                                          detected_sample_outliers= outliers)
   }
 

diff --git a/R/proteomics_report.R b/R/proteomics_report.R
@@ -0,0 +1,70 @@
+#' Generate a report for single omic analysis
+#'
+#' @param multiassay Multiassay experiment object generated by Omix
+#' @param report_folder_path report_folder_path folder path to save the report.
+#' @param report_file  filename for report (without an extension).
+#' @param slot Slot of interest
+#' @param database Enrichment database `GO_Molecular_Function_2021`,`GO_Cellular_Component_2021`,
+#'  `GO_Biological_Process_2021`, `Reactome_2016` , `KEGG_2021_Human` , `MSigDB_Hallmark_2020`
+#' @param clinical_covariates clinical covariates for table1
+#' @family Report
+#' @return Single omic analyses report html
+#' @export
+#'
+proteomics_report <- function(multiassay,
+                               report_folder_path = getwd(),
+                               report_file = "single_omic_report_Omix",
+                               slot='ADvsControl',
+                               database='Reactome_2016',
+                               num_path=20,
+                               clinical_covariates=c('PHF1','amyloid','PMD')){
+  uniomic=list()
+  uniomic=multiassay@metadata[c('DEP')]
+  uniomic$param$slot=slot
+  uniomic$param$database=database
+  uniomic$param$num_path=num_path
+  uniomic$metadata=multiassay@colData
+  uniomic$parameters_analysis_protein = multiassay@metadata$parameters$single_omic$de$protein
+  uniomic$parameters_processing_protein = multiassay@metadata$parameters$processing$protein
+  uniomic$clinical_covariates= clinical_covariates
+
+  report_file <- tools::file_path_sans_ext(report_file)
+
+  cli::cli_h2("Generating report for single omic analyses")
+
+
+  metadata_tmp_path <- file.path(tempdir(), "metadata.qs")
+
+  cli::cli_text("Writing temp files for report...")
+  qs::qsave(
+    uniomic,
+    metadata_tmp_path
+  )
+
+  krd <- file.path(tempdir(), "krdqc")
+  intd <- file.path(tempdir(), "idqc")
+  dir.create(krd, showWarnings = FALSE)
+  dir.create(intd, showWarnings = FALSE)
+
+  cli::cli_text("Generating single omic analyses report...")
+  rmarkdown::render(
+    system.file(
+      "rmarkdown/templates/proteomics/skeleton.Rmd",
+      package = "Omix"
+    ),
+    params = list(
+      metadata_path = metadata_tmp_path
+    ),
+    output_dir = report_folder_path,
+    output_file = report_file,
+    knit_root_dir = krd,
+    intermediates_dir = intd,
+    quiet = TRUE
+  )
+
+  report_file_name <- paste(report_file, ".html", sep = "")
+
+
+
+}
+