Skip to content

Commit c2b019d

Browse files
stemangiolastemangiola
authored
Merge pull request #268 from mikelove/master
Add LFC threshold tests for DESeq2
2 parents 279b53c + 4babb46 commit c2b019d

File tree

6 files changed

+107
-54
lines changed

6 files changed

+107
-54
lines changed

‎R/functions.R‎

Lines changed: 21 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -873,19 +873,20 @@ get_differential_transcript_abundance_bulk_voom <- function(.data,
873873
#' @param omit_contrast_in_colnames If just one contrast is specified you can choose to omit the contrast label in the colnames.
874874
#' @param ... Additional arguments for DESeq2
875875
#'
876-
#' @return A tibble with edgeR results
876+
#' @return A tibble with DESeq2 results
877877
#'
878878
get_differential_transcript_abundance_deseq2 <- function(.data,
879-
.formula,
880-
.sample = NULL,
881-
.transcript = NULL,
882-
.abundance = NULL,
883-
.contrasts = NULL,
884-
method = "edgeR_quasi_likelihood",
885-
scaling_method = "TMM",
886-
omit_contrast_in_colnames = FALSE,
887-
prefix = "",
888-
...) {
879+
.formula,
880+
.sample = NULL,
881+
.transcript = NULL,
882+
.abundance = NULL,
883+
.contrasts = NULL,
884+
method = "deseq2",
885+
test_above_log2_fold_change = NULL,
886+
scaling_method = "TMM",
887+
omit_contrast_in_colnames = FALSE,
888+
prefix = "",
889+
...) {
889890

890891
# Check if contrasts are of the same form
891892
if(
@@ -905,11 +906,6 @@ get_differential_transcript_abundance_deseq2 <- function(.data,
905906
omit_contrast_in_colnames = FALSE
906907
}
907908

908-
if (find.package("acepack", quiet = TRUE) %>% length %>% equals(0)) {
909-
message("tidybulk says: Installing acepack needed for analyses")
910-
install.packages("acepack", repos = "https://cloud.r-project.org")
911-
}
912-
913909
# Check if package is installed, otherwise install
914910
if (find.package("DESeq2", quiet = TRUE) %>% length %>% equals(0)) {
915911
message("tidybulk says: Installing DESeq2 needed for differential transcript abundance analyses")
@@ -918,6 +914,10 @@ get_differential_transcript_abundance_deseq2 <- function(.data,
918914
BiocManager::install("DESeq2", ask = FALSE)
919915
}
920916

917+
if (is.null(test_above_log2_fold_change)) {
918+
test_above_log2_fold_change <- 0
919+
}
920+
921921
# # Print the design column names in case I want contrasts
922922
# message(
923923
# sprintf(
@@ -952,7 +952,7 @@ get_differential_transcript_abundance_deseq2 <- function(.data,
952952
tidybulk_to_SummarizedExperiment(!!.sample, !!.transcript, counts) %>%
953953

954954
# DESeq2
955-
DESeq2::DESeqDataSet( design = .formula) %>%
955+
DESeq2::DESeqDataSet(design = .formula) %>%
956956
DESeq2::DESeq(...)
957957

958958
# Read ft object
@@ -966,7 +966,7 @@ get_differential_transcript_abundance_deseq2 <- function(.data,
966966
(deseq2_object@colData[,parse_formula(.formula)[1]] %>%
967967
class %in% c("numeric", "integer", "double")) ~
968968
(.) %>%
969-
DESeq2::results() %>%
969+
DESeq2::results(lfcThreshold=test_above_log2_fold_change) %>%
970970
as_tibble(rownames = quo_name(.transcript)),
971971

972972
# Simple comparison discrete
@@ -976,13 +976,13 @@ get_differential_transcript_abundance_deseq2 <- function(.data,
976976
parse_formula(.formula)[1],
977977
deseq2_object@colData[,parse_formula(.formula)[1]] %>% levels %>% .[2],
978978
deseq2_object@colData[,parse_formula(.formula)[1]] %>% levels %>% .[1]
979-
)) %>%
979+
), lfcThreshold=test_above_log2_fold_change) %>%
980980
as_tibble(rownames = quo_name(.transcript)),
981981

982982
# Simple comparison discrete
983983
my_contrasts %>% is.null %>% not() & omit_contrast_in_colnames ~
984984
(.) %>%
985-
DESeq2::results(contrast = my_contrasts[[1]])%>%
985+
DESeq2::results(contrast = my_contrasts[[1]], lfcThreshold=test_above_log2_fold_change)%>%
986986
as_tibble(rownames = quo_name(.transcript)),
987987

988988
# Multiple comparisons NOT USED AT THE MOMENT
@@ -994,7 +994,7 @@ get_differential_transcript_abundance_deseq2 <- function(.data,
994994
~ deseq2_obj %>%
995995

996996
# select method
997-
DESeq2::results(contrast = my_contrasts[[.x]]) %>%
997+
DESeq2::results(contrast = my_contrasts[[.x]], lfcThreshold=test_above_log2_fold_change) %>%
998998

999999
# Convert to tibble
10001000
as_tibble(rownames = quo_name(.transcript)) %>%

‎R/functions_SE.R‎

Lines changed: 20 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -1052,16 +1052,19 @@ get_differential_transcript_abundance_bulk_voom_SE <- function(.data,
10521052
#' @param omit_contrast_in_colnames If just one contrast is specified you can choose to omit the contrast label in the colnames.
10531053
#' @param ... Additional arguments for DESeq2
10541054
#'
1055-
#' @return A tibble with edgeR results
1055+
#' @return A tibble with DESeq2 results
10561056
#'
10571057
get_differential_transcript_abundance_deseq2_SE <- function(.data,
1058-
.formula,
1059-
.contrasts = NULL,
1060-
method = "edgeR_quasi_likelihood",
1061-
scaling_method = "TMM",
1062-
omit_contrast_in_colnames = FALSE,
1063-
prefix = "",
1064-
...) {
1058+
.formula,
1059+
.contrasts = NULL,
1060+
method = "deseq2",
1061+
1062+
test_above_log2_fold_change = NULL,
1063+
1064+
scaling_method = "TMM",
1065+
omit_contrast_in_colnames = FALSE,
1066+
prefix = "",
1067+
...) {
10651068

10661069

10671070
# Check if contrasts are of the same form
@@ -1077,11 +1080,6 @@ get_differential_transcript_abundance_deseq2_SE <- function(.data,
10771080
omit_contrast_in_colnames = FALSE
10781081
}
10791082

1080-
if (find.package("acepack", quiet = TRUE) %>% length %>% equals(0)) {
1081-
message("Installing acepack needed for analyses")
1082-
install.packages("acepack", repos = "https://cloud.r-project.org")
1083-
}
1084-
10851083
# Check if package is installed, otherwise install
10861084
if (find.package("DESeq2", quiet = TRUE) %>% length %>% equals(0)) {
10871085
message("Installing DESeq2 needed for differential transcript abundance analyses")
@@ -1090,14 +1088,18 @@ get_differential_transcript_abundance_deseq2_SE <- function(.data,
10901088
BiocManager::install("DESeq2", ask = FALSE)
10911089
}
10921090

1091+
if (is.null(test_above_log2_fold_change)) {
1092+
test_above_log2_fold_change <- 0
1093+
}
1094+
10931095
my_contrasts = .contrasts
10941096

10951097
deseq2_object =
10961098
.data %>%
10971099

10981100
# DESeq2
10991101
DESeq2::DESeqDataSet( design = .formula) %>%
1100-
DESeq2::DESeq()
1102+
DESeq2::DESeq(...)
11011103

11021104
# Return
11031105
list(
@@ -1114,7 +1116,7 @@ get_differential_transcript_abundance_deseq2_SE <- function(.data,
11141116
(deseq2_object@colData[,parse_formula(.formula)[1]] %>%
11151117
class %in% c("numeric", "integer", "double")) ~
11161118
(.) %>%
1117-
DESeq2::results() %>%
1119+
DESeq2::results(lfcThreshold=test_above_log2_fold_change) %>%
11181120
as_tibble(rownames = "transcript"),
11191121

11201122
# Simple comparison discrete
@@ -1124,13 +1126,13 @@ get_differential_transcript_abundance_deseq2_SE <- function(.data,
11241126
parse_formula(.formula)[1],
11251127
deseq2_object@colData[,parse_formula(.formula)[1]] %>% as.factor() %>% levels %>% .[2],
11261128
deseq2_object@colData[,parse_formula(.formula)[1]] %>% as.factor() %>% levels %>% .[1]
1127-
)) %>%
1129+
), lfcThreshold=test_above_log2_fold_change) %>%
11281130
as_tibble(rownames = "transcript"),
11291131

11301132
# Simple comparison discrete
11311133
my_contrasts %>% is.null %>% not() & omit_contrast_in_colnames ~
11321134
(.) %>%
1133-
DESeq2::results(contrast = my_contrasts[[1]])%>%
1135+
DESeq2::results(contrast = my_contrasts[[1]], lfcThreshold=test_above_log2_fold_change)%>%
11341136
as_tibble(rownames = "transcript"),
11351137

11361138
# Multiple comparisons NOT USED AT THE MOMENT
@@ -1142,7 +1144,7 @@ get_differential_transcript_abundance_deseq2_SE <- function(.data,
11421144
~ deseq2_obj %>%
11431145

11441146
# select method
1145-
DESeq2::results(contrast = my_contrasts[[.x]]) %>%
1147+
DESeq2::results(contrast = my_contrasts[[.x]], lfcThreshold=test_above_log2_fold_change) %>%
11461148

11471149
# Convert to tibble
11481150
as_tibble(rownames = "transcript") %>%

‎R/methods.R‎

Lines changed: 23 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -2290,7 +2290,7 @@ setMethod("ensembl_to_symbol", "tidybulk", .ensembl_to_symbol)
22902290
#' ) %>%
22912291
#'
22922292
#' # DESeq2
2293-
#' DESeq2::DESeqDataSet( design = .formula) %>%
2293+
#' DESeq2::DESeqDataSet(design = .formula) %>%
22942294
#' DESeq2::DESeq() %>%
22952295
#' DESeq2::results()
22962296
#'
@@ -2322,18 +2322,35 @@ setMethod("ensembl_to_symbol", "tidybulk", .ensembl_to_symbol)
23222322
#' my_se_mini = tidybulk::se_mini
23232323
#' my_se_mini$condition = factor(my_se_mini$condition)
23242324
#'
2325+
#' # demontrating with `fitType` that you can access any arguments to DESeq()
23252326
#' my_se_mini |>
2326-
#' identify_abundant() |>
2327-
#' test_differential_abundance( ~ condition, method="deseq2" )
2327+
#' identify_abundant() |>
2328+
#' test_differential_abundance( ~ condition, method="deseq2", fitType="local")
23282329
#'
2329-
#' # The function `test_differential_abundance` operates with contrasts too
2330+
#' # testing above a log2 threshold, passes along value to lfcThreshold of results()
2331+
#' res <- my_se_mini |>
2332+
#' identify_abundant() |>
2333+
#' test_differential_abundance( ~ condition, method="deseq2",
2334+
#' fitType="local",
2335+
#' test_above_log2_fold_change=4 )
2336+
#'
2337+
#' # confirm that lfcThreshold was used
2338+
#' \dontrun{
2339+
#' res |>
2340+
#' mcols() |>
2341+
#' DESeq2::DESeqResults() |>
2342+
#' DESeq2::plotMA()
2343+
#' }
2344+
#'
2345+
#' # The function `test_differential_abundance` operates with contrasts too
23302346
#'
23312347
#' my_se_mini |>
23322348
#' identify_abundant() |>
23332349
#' test_differential_abundance(
23342350
#' ~ 0 + condition,
23352351
#' contrasts = list(c("condition", "TRUE", "FALSE")),
2336-
#' method="deseq2"
2352+
#' method="deseq2",
2353+
#' fitType="local"
23372354
#' )
23382355
#'
23392356
#' @docType methods
@@ -2492,6 +2509,7 @@ such as batch effects (if applicable) in the formula.
24922509
.abundance = !!.abundance,
24932510
.contrasts = contrasts,
24942511
method = method,
2512+
test_above_log2_fold_change = test_above_log2_fold_change,
24952513
scaling_method = scaling_method,
24962514
omit_contrast_in_colnames = omit_contrast_in_colnames,
24972515
prefix = prefix,

‎R/methods_SE.R‎

Lines changed: 6 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1164,7 +1164,8 @@ such as batch effects (if applicable) in the formula.
11641164
test_above_log2_fold_change = test_above_log2_fold_change,
11651165
scaling_method = scaling_method,
11661166
omit_contrast_in_colnames = omit_contrast_in_colnames,
1167-
prefix = prefix
1167+
prefix = prefix,
1168+
...
11681169
),
11691170

11701171
# DESeq2
@@ -1173,9 +1174,11 @@ such as batch effects (if applicable) in the formula.
11731174
.formula,
11741175
.contrasts = contrasts,
11751176
method = method,
1177+
test_above_log2_fold_change = test_above_log2_fold_change,
11761178
scaling_method = scaling_method,
11771179
omit_contrast_in_colnames = omit_contrast_in_colnames,
1178-
prefix = prefix
1180+
prefix = prefix,
1181+
...
11791182
),
11801183

11811184
# Else error
@@ -1201,7 +1204,7 @@ such as batch effects (if applicable) in the formula.
12011204
tolower(method) == "limma_voom" ~ (.) %>% memorise_methods_used("voom"),
12021205
tolower(method) == "limma_voom_sample_weights" ~ (.) %>% memorise_methods_used("voom_sample_weights"),
12031206
tolower(method) == "deseq2" ~ (.) %>% memorise_methods_used("deseq2"),
1204-
~ 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)")
1207+
~ stop("tidybulk says: method not supported")
12051208
) %>%
12061209

12071210
when(

‎man/test_differential_abundance-methods.Rd‎

Lines changed: 22 additions & 5 deletions
Some generated files are not rendered by default. Learn more about customizing how changed files appear on GitHub.

‎tests/testthat/test-bulk_methods.R‎

Lines changed: 15 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -644,7 +644,7 @@ test_that("DESeq2 differential trancript abundance - no object",{
644644
test_deseq2_df |>
645645
tidybulk() |>
646646
identify_abundant(factor_of_interest = condition) |>
647-
test_differential_abundance(~condition, method="DeSEQ2", action="get")
647+
test_differential_abundance(~condition, method="DESeq2", action="get")
648648

649649

650650
expect_equal(
@@ -790,6 +790,19 @@ test_that("DESeq2 differential trancript abundance - no object",{
790790
0
791791
)
792792

793+
# DESeq2 with lfcThreshold using test_above_log2_fold_change
794+
res_lfc <- test_deseq2_df |>
795+
keep_abundant(factor_of_interest = condition) |>
796+
test_differential_abundance(~condition, method="DESeq2", test_above_log2_fold_change=1) |>
797+
mcols()
798+
799+
# DESeq2::plotMA(DESeq2::DESeqResults(res_lfc))
800+
801+
# significant are outside of LFC 1 / -1
802+
idx <- which(res_lfc$padj < .1)
803+
expect_true(all(abs(res_lfc$log2FoldChange[idx]) > 1))
804+
805+
793806
})
794807

795808
test_that("test prefix",{
@@ -801,7 +814,7 @@ test_that("test prefix",{
801814

802815
res_DeSEQ2 =
803816
df |>
804-
test_differential_abundance(~condition, method="DeSEQ2", action="only", prefix = "prefix_")
817+
test_differential_abundance(~condition, method="DESeq2", action="only", prefix = "prefix_")
805818

806819
res_voom =
807820
df |>

0 commit comments

Comments
 (0)