Merge pull request #268 from mikelove/master

Add LFC threshold tests for DESeq2
stemangiola · Mar 3, 2023 · c2b019d · c2b019d
2 parents 279b53c + 4babb46
commit c2b019d
Show file tree

Hide file tree

Showing 6 changed files with 107 additions and 54 deletions.
diff --git a/R/functions.R b/R/functions.R
@@ -873,19 +873,20 @@ get_differential_transcript_abundance_bulk_voom <- function(.data,
 #' @param omit_contrast_in_colnames If just one contrast is specified you can choose to omit the contrast label in the colnames.
 #' @param ... Additional arguments for DESeq2
 #'
-#' @return A tibble with edgeR results
+#' @return A tibble with DESeq2 results
 #'
 get_differential_transcript_abundance_deseq2 <- function(.data,
-																											 .formula,
-																											 .sample = NULL,
-																											 .transcript = NULL,
-																											 .abundance = NULL,
-																											 .contrasts = NULL,
-																											 method = "edgeR_quasi_likelihood",
-																											 scaling_method = "TMM",
-																											 omit_contrast_in_colnames = FALSE,
-																											 prefix = "",
-																											 ...) {
+                                                         .formula,
+                                                         .sample = NULL,
+                                                         .transcript = NULL,
+                                                         .abundance = NULL,
+                                                         .contrasts = NULL,
+                                                         method = "deseq2",
+                                                         test_above_log2_fold_change = NULL,
+                                                         scaling_method = "TMM",
+                                                         omit_contrast_in_colnames = FALSE,
+                                                         prefix = "",
+                                                         ...) {
 
 	# Check if contrasts are of the same form
 	if(
@@ -905,11 +906,6 @@ get_differential_transcript_abundance_deseq2 <- function(.data,
 		omit_contrast_in_colnames = FALSE
 	}
 
-	if (find.package("acepack", quiet = TRUE) %>% length %>% equals(0)) {
-		message("tidybulk says: Installing acepack needed for analyses")
-		install.packages("acepack", repos = "https://cloud.r-project.org")
-	}
-
 	# Check if package is installed, otherwise install
 	if (find.package("DESeq2", quiet = TRUE) %>% length %>% equals(0)) {
 		message("tidybulk says: Installing DESeq2 needed for differential transcript abundance analyses")
@@ -918,6 +914,10 @@ get_differential_transcript_abundance_deseq2 <- function(.data,
 		BiocManager::install("DESeq2", ask = FALSE)
 	}
 
+        if (is.null(test_above_log2_fold_change)) {
+          test_above_log2_fold_change <- 0
+        }
+
 	# # Print the design column names in case I want contrasts
 	# message(
 	# 	sprintf(
@@ -952,7 +952,7 @@ get_differential_transcript_abundance_deseq2 <- function(.data,
 		tidybulk_to_SummarizedExperiment(!!.sample, !!.transcript, counts) %>%
 
 		# DESeq2
-		DESeq2::DESeqDataSet( design = .formula) %>%
+		DESeq2::DESeqDataSet(design = .formula) %>%
 		DESeq2::DESeq(...)
 
 	# Read ft object
@@ -966,7 +966,7 @@ get_differential_transcript_abundance_deseq2 <- function(.data,
 				(deseq2_object@colData[,parse_formula(.formula)[1]] %>%
 				 	class %in% c("numeric", "integer", "double")) 	~
 				(.) %>%
-				DESeq2::results() %>%
+				DESeq2::results(lfcThreshold=test_above_log2_fold_change) %>%
 				as_tibble(rownames = quo_name(.transcript)),
 
 			# Simple comparison discrete
@@ -976,13 +976,13 @@ get_differential_transcript_abundance_deseq2 <- function(.data,
 					parse_formula(.formula)[1],
 					deseq2_object@colData[,parse_formula(.formula)[1]] %>% levels %>% .[2],
 					deseq2_object@colData[,parse_formula(.formula)[1]] %>% levels %>% .[1]
-				)) %>%
+				), lfcThreshold=test_above_log2_fold_change) %>%
 				as_tibble(rownames = quo_name(.transcript)),
 
 			# Simple comparison discrete
 			my_contrasts %>% is.null %>% not() & omit_contrast_in_colnames	~
 				(.) %>%
-				DESeq2::results(contrast = my_contrasts[[1]])%>%
+				DESeq2::results(contrast = my_contrasts[[1]], lfcThreshold=test_above_log2_fold_change)%>%
 				as_tibble(rownames = quo_name(.transcript)),
 
 			# Multiple comparisons NOT USED AT THE MOMENT
@@ -994,7 +994,7 @@ get_differential_transcript_abundance_deseq2 <- function(.data,
 						~ 	deseq2_obj %>%
 
 							# select method
-							DESeq2::results(contrast = my_contrasts[[.x]])	%>%
+							DESeq2::results(contrast = my_contrasts[[.x]], lfcThreshold=test_above_log2_fold_change)	%>%
 
 							# Convert to tibble
 							as_tibble(rownames = quo_name(.transcript)) %>%

diff --git a/R/functions_SE.R b/R/functions_SE.R
@@ -1052,16 +1052,19 @@ get_differential_transcript_abundance_bulk_voom_SE <- function(.data,
 #' @param omit_contrast_in_colnames If just one contrast is specified you can choose to omit the contrast label in the colnames.
 #' @param ... Additional arguments for DESeq2
 #'
-#' @return A tibble with edgeR results
+#' @return A tibble with DESeq2 results
 #'
 get_differential_transcript_abundance_deseq2_SE <- function(.data,
-																												 .formula,
-																												 .contrasts = NULL,
-																												 method = "edgeR_quasi_likelihood",
-																												 scaling_method = "TMM",
-																												 omit_contrast_in_colnames = FALSE,
-																												 prefix = "",
-																												 ...) {
+                                                            .formula,
+                                                            .contrasts = NULL,
+                                                            method = "deseq2",
+
+                                                            test_above_log2_fold_change = NULL,
+
+                                                            scaling_method = "TMM",
+                                                            omit_contrast_in_colnames = FALSE,
+                                                            prefix = "",
+                                                            ...) {
 
 
 	# Check if contrasts are of the same form
@@ -1077,11 +1080,6 @@ get_differential_transcript_abundance_deseq2_SE <- function(.data,
 		omit_contrast_in_colnames = FALSE
 	}
 
-	if (find.package("acepack", quiet = TRUE) %>% length %>% equals(0)) {
-		message("Installing acepack needed for analyses")
-		install.packages("acepack", repos = "https://cloud.r-project.org")
-	}
-
 	# Check if package is installed, otherwise install
 	if (find.package("DESeq2", quiet = TRUE) %>% length %>% equals(0)) {
 		message("Installing DESeq2 needed for differential transcript abundance analyses")
@@ -1090,14 +1088,18 @@ get_differential_transcript_abundance_deseq2_SE <- function(.data,
 		BiocManager::install("DESeq2", ask = FALSE)
 	}
 
+        if (is.null(test_above_log2_fold_change)) {
+          test_above_log2_fold_change <- 0
+        }
+
 	my_contrasts = .contrasts
 
 	deseq2_object =
 		.data %>%
 
 		# DESeq2
 		DESeq2::DESeqDataSet( design = .formula) %>%
-		DESeq2::DESeq()
+		DESeq2::DESeq(...)
 
 	# Return
 	list(
@@ -1114,7 +1116,7 @@ get_differential_transcript_abundance_deseq2_SE <- function(.data,
 					(deseq2_object@colData[,parse_formula(.formula)[1]] %>%
 					 	class %in% c("numeric", "integer", "double")) 	~
 					(.) %>%
-					DESeq2::results() %>%
+					DESeq2::results(lfcThreshold=test_above_log2_fold_change) %>%
 					as_tibble(rownames = "transcript"),
 
 				# Simple comparison discrete
@@ -1124,13 +1126,13 @@ get_differential_transcript_abundance_deseq2_SE <- function(.data,
 						parse_formula(.formula)[1],
 						deseq2_object@colData[,parse_formula(.formula)[1]] %>% as.factor() %>% levels %>% .[2],
 						deseq2_object@colData[,parse_formula(.formula)[1]] %>% as.factor() %>% levels %>% .[1]
-					)) %>%
+					), lfcThreshold=test_above_log2_fold_change) %>%
 					as_tibble(rownames = "transcript"),
 
 				# Simple comparison discrete
 				my_contrasts %>% is.null %>% not() & omit_contrast_in_colnames	~
 					(.) %>%
-					DESeq2::results(contrast = my_contrasts[[1]])%>%
+					DESeq2::results(contrast = my_contrasts[[1]], lfcThreshold=test_above_log2_fold_change)%>%
 					as_tibble(rownames = "transcript"),
 
 				# Multiple comparisons NOT USED AT THE MOMENT
@@ -1142,7 +1144,7 @@ get_differential_transcript_abundance_deseq2_SE <- function(.data,
 							~ 	deseq2_obj %>%
 
 								# select method
-								DESeq2::results(contrast = my_contrasts[[.x]])	%>%
+								DESeq2::results(contrast = my_contrasts[[.x]], lfcThreshold=test_above_log2_fold_change)	%>%
 
 								# Convert to tibble
 								as_tibble(rownames = "transcript") %>%

diff --git a/R/methods.R b/R/methods.R
@@ -2290,7 +2290,7 @@ setMethod("ensembl_to_symbol", "tidybulk", .ensembl_to_symbol)
 #'	) %>%
 #'
 #'	# DESeq2
-#'	DESeq2::DESeqDataSet( design = .formula) %>%
+#'	DESeq2::DESeqDataSet(design = .formula) %>%
 #'	DESeq2::DESeq() %>%
 #'	DESeq2::results()
 #'
@@ -2322,18 +2322,35 @@ setMethod("ensembl_to_symbol", "tidybulk", .ensembl_to_symbol)
 #' my_se_mini = tidybulk::se_mini
 #' my_se_mini$condition  = factor(my_se_mini$condition)
 #'
+#' # demontrating with `fitType` that you can access any arguments to DESeq()
 #' my_se_mini  |>
-#'  identify_abundant() |>
-#' 	test_differential_abundance( ~ condition, method="deseq2" )
+#'    identify_abundant() |>
+#'        test_differential_abundance( ~ condition, method="deseq2", fitType="local")
 #'
-#' 	# The function `test_differential_abundance` operates with contrasts too
+#' # testing above a log2 threshold, passes along value to lfcThreshold of results()
+#' res <- my_se_mini  |>
+#'    identify_abundant() |>
+#'         test_differential_abundance( ~ condition, method="deseq2",
+#'             fitType="local",
+#'             test_above_log2_fold_change=4 )
+#'
+#' # confirm that lfcThreshold was used
+#' \dontrun{
+#'     res |>
+#'         mcols() |>
+#'         DESeq2::DESeqResults() |>
+#'         DESeq2::plotMA()
+#' }
+#' 
+#' # The function `test_differential_abundance` operates with contrasts too
 #'
 #'  my_se_mini |>
 #'  identify_abundant() |>
 #'  test_differential_abundance(
 #' 	    ~ 0 + condition,
 #' 	    contrasts = list(c("condition", "TRUE", "FALSE")),
-#' 	    method="deseq2"
+#' 	    method="deseq2",
+#'          fitType="local"
 #'  )
 #'
 #' @docType methods
@@ -2492,6 +2509,7 @@ such as batch effects (if applicable) in the formula.
 				.abundance = !!.abundance,
 				.contrasts = contrasts,
 				method = method,
+				test_above_log2_fold_change = test_above_log2_fold_change,                                
 				scaling_method = scaling_method,
 				omit_contrast_in_colnames = omit_contrast_in_colnames,
 				prefix = prefix,

diff --git a/R/methods_SE.R b/R/methods_SE.R
@@ -1164,7 +1164,8 @@ such as batch effects (if applicable) in the formula.
 				test_above_log2_fold_change = test_above_log2_fold_change,
 				scaling_method = scaling_method,
 				omit_contrast_in_colnames = omit_contrast_in_colnames,
-				prefix = prefix
+                                prefix = prefix,
+                                ...
 			),
 
 			# DESeq2
@@ -1173,9 +1174,11 @@ such as batch effects (if applicable) in the formula.
 				.formula,
 				.contrasts = contrasts,
 				method = method,
+                                test_above_log2_fold_change = test_above_log2_fold_change,
 				scaling_method = scaling_method,
 				omit_contrast_in_colnames = omit_contrast_in_colnames,
-				prefix = prefix
+				prefix = prefix,
+                                ...
 			),
 
 			# Else error
@@ -1201,7 +1204,7 @@ such as batch effects (if applicable) in the formula.
 			tolower(method) == "limma_voom" ~ (.) %>% memorise_methods_used("voom"),
 			tolower(method) == "limma_voom_sample_weights" ~ (.) %>% memorise_methods_used("voom_sample_weights"),
 			tolower(method) == "deseq2" ~ (.) %>% memorise_methods_used("deseq2"),
-			~ stop("tidybulk says: method must be either \"correlation\" for dropping correlated elements or \"reduced_dimension\" to drop the closest pair according to two dimensions (e.g., PCA)")
+			~ stop("tidybulk says: method not supported")
 		) %>%
 
 	    when(

diff --git a/man/test_differential_abundance-methods.Rd b/man/test_differential_abundance-methods.Rd
diff --git a/tests/testthat/test-bulk_methods.R b/tests/testthat/test-bulk_methods.R
@@ -644,7 +644,7 @@ test_that("DESeq2 differential trancript abundance - no object",{
 		test_deseq2_df |>
 		tidybulk() |>
 		 identify_abundant(factor_of_interest = condition) |>
-		test_differential_abundance(~condition, method="DeSEQ2", action="get")
+		test_differential_abundance(~condition, method="DESeq2", action="get")
 
 
 	expect_equal(
@@ -790,6 +790,19 @@ test_that("DESeq2 differential trancript abundance - no object",{
 		0
 	)
 
+        # DESeq2 with lfcThreshold using test_above_log2_fold_change
+        res_lfc <- test_deseq2_df |>
+          keep_abundant(factor_of_interest = condition) |>
+          test_differential_abundance(~condition, method="DESeq2", test_above_log2_fold_change=1) |>
+          mcols()
+
+        # DESeq2::plotMA(DESeq2::DESeqResults(res_lfc))
+
+        # significant are outside of LFC 1 / -1
+        idx <- which(res_lfc$padj < .1)
+        expect_true(all(abs(res_lfc$log2FoldChange[idx]) > 1))
+
+
 })
 
 test_that("test prefix",{
@@ -801,7 +814,7 @@ test_that("test prefix",{
 
 	res_DeSEQ2 =
 		df |>
-		test_differential_abundance(~condition, method="DeSEQ2", action="only", prefix = "prefix_")
+		test_differential_abundance(~condition, method="DESeq2", action="only", prefix = "prefix_")
 
 	res_voom =
 		df |>