Skip to content

Commit

Permalink
Merge pull request #274 from stemangiola/allow-custom-assay-for-SE-DE
Browse files Browse the repository at this point in the history 
Allow custom assay for se de
  • Loading branch information
@stemangiola
stemangiola authored Jul 6, 2023
2 parents 43219c8 + 35b4fb5 commit 136461b
Show file tree
Hide file tree
Showing 4 changed files with 227 additions and 32 deletions.
1 change: 1 addition & 0 deletions NAMESPACE
View file
Original file line number Diff line number Diff line change
Expand Up @@ -125,6 +125,7 @@ importFrom(magrittr,"%$%")
importFrom(magrittr,"%>%")
importFrom(magrittr,divide_by)
importFrom(magrittr,equals)
importFrom(magrittr,extract2)
importFrom(magrittr,multiply_by)
importFrom(magrittr,set_colnames)
importFrom(magrittr,set_rownames)
Expand Down
138 changes: 120 additions & 18 deletions R/functions_SE.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -701,6 +701,7 @@ remove_redundancy_elements_though_reduced_dimensions_SE <-
#' @importFrom stats model.matrix
#' @importFrom utils install.packages
#' @importFrom purrr when
#' @importFrom magrittr extract2
#'
#'
#' @param .data A tibble
Expand All @@ -715,6 +716,7 @@ remove_redundancy_elements_though_reduced_dimensions_SE <-
#'
get_differential_transcript_abundance_bulk_SE <- function(.data,
.formula,
.abundance = NULL,
sample_annotation,
.contrasts = NULL,
method = "edgeR_quasi_likelihood",
Expand All @@ -724,6 +726,8 @@ get_differential_transcript_abundance_bulk_SE <- function(.data,
prefix = "",
...) {

.abundance = enquo(.abundance)

# Check if omit_contrast_in_colnames is correctly setup
if(omit_contrast_in_colnames & length(.contrasts) > 1){
warning("tidybulk says: you can omit contrasts in column names only when maximum one contrast is present")
Expand Down Expand Up @@ -759,15 +763,19 @@ get_differential_transcript_abundance_bulk_SE <- function(.data,
BiocManager::install("edgeR", ask = FALSE)
}

# If no assay is specified take first
my_assay = ifelse(
quo_is_symbol(.abundance),
quo_name(.abundance),
.data |>
assayNames() |>
extract2(1)
)

edgeR_object =
.data %>%

# Extract assay
assays() %>%
as.list() %>%
.[[1]] %>%

edgeR::DGEList(counts = .) %>%
.data |>
assay(my_assay) |>
edgeR::DGEList() %>%

# Scale data if method is not "none"
when(
Expand Down Expand Up @@ -881,6 +889,7 @@ get_differential_transcript_abundance_bulk_SE <- function(.data,
#'
get_differential_transcript_abundance_bulk_voom_SE <- function(.data,
.formula,
.abundance = NULL,
sample_annotation,
.contrasts = NULL,
method = NULL,
Expand All @@ -889,6 +898,7 @@ get_differential_transcript_abundance_bulk_voom_SE <- function(.data,
omit_contrast_in_colnames = FALSE,
prefix = "") {

.abundance = enquo(.abundance)

# Check if omit_contrast_in_colnames is correctly setup
if(omit_contrast_in_colnames & length(.contrasts) > 1){
Expand Down Expand Up @@ -926,15 +936,20 @@ get_differential_transcript_abundance_bulk_voom_SE <- function(.data,
BiocManager::install("limma", ask = FALSE)
}

# If no assay is specified take first
my_assay = ifelse(
quo_is_symbol(.abundance),
quo_name(.abundance),
.data |>
assayNames() |>
extract2(1)
)

voom_object =
.data %>%

# Extract assay
assays() %>%
as.list() %>%
.[[1]] %>%

edgeR::DGEList() %>%
assay(my_assay) |>
edgeR::DGEList() %>%

# Scale data if method is not "none"
when(
Expand Down Expand Up @@ -1027,6 +1042,79 @@ get_differential_transcript_abundance_bulk_voom_SE <- function(.data,
)


}


#' Get differential transcription information to a tibble using glmmSeq
#'
#' @keywords internal
#' @noRd
#'
#'
#'
#' @import tibble
#' @importFrom magrittr set_colnames
#' @importFrom stats model.matrix
#' @importFrom utils install.packages
#' @importFrom purrr when
#'
#'
#' @param .data A tibble
#' @param .formula a formula with no response variable, referring only to numeric variables
#' @param .contrasts A character vector. See edgeR makeContrasts specification for the parameter `contrasts`. If contrasts are not present the first covariate is the one the model is tested against (e.g., ~ factor_of_interest)
#' @param method A string character. Either "edgeR_quasi_likelihood" (i.e., QLF), "edgeR_likelihood_ratio" (i.e., LRT)
#' @param scaling_method A character string. The scaling method passed to the backend function (i.e., edgeR::calcNormFactors; "TMM","TMMwsp","RLE","upperquartile")
#' @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 glmmSeq
#'
#' @return A tibble with glmmSeq results
#'
get_differential_transcript_abundance_glmmSeq_SE <- function(.data,
.formula,
.abundance = NULL,
.contrasts = NULL,
sample_annotation ,
method = "deseq2",

test_above_log2_fold_change = NULL,

scaling_method = "TMM",
omit_contrast_in_colnames = FALSE,
prefix = "",
...) {

.abundance = enquo(.abundance)

# Check if contrasts are of the same form
if(
.contrasts %>% is.null %>% not() &
.contrasts %>% class %>% equals("list") %>% not()
)
stop("tidybulk says: for DESeq2 the list of constrasts should be given in the form list(c(\"condition_column\",\"condition1\",\"condition2\")) i.e. list(c(\"genotype\",\"knockout\",\"wildtype\"))")

# Check if omit_contrast_in_colnames is correctly setup
if(omit_contrast_in_colnames & length(.contrasts) > 1){
warning("tidybulk says: you can omit contrasts in column names only when maximum one contrast is present")
omit_contrast_in_colnames = FALSE
}

# 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")
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager", repos = "https://cloud.r-project.org")
BiocManager::install("DESeq2", ask = FALSE)
}

if (is.null(test_above_log2_fold_change)) {
test_above_log2_fold_change <- 0
}






}


Expand Down Expand Up @@ -1056,6 +1144,7 @@ get_differential_transcript_abundance_bulk_voom_SE <- function(.data,
#'
get_differential_transcript_abundance_deseq2_SE <- function(.data,
.formula,
.abundance = NULL,
.contrasts = NULL,
method = "deseq2",

Expand All @@ -1066,7 +1155,8 @@ get_differential_transcript_abundance_deseq2_SE <- function(.data,
prefix = "",
...) {


.abundance = enquo(.abundance)

# Check if contrasts are of the same form
if(
.contrasts %>% is.null %>% not() &
Expand Down Expand Up @@ -1094,11 +1184,23 @@ get_differential_transcript_abundance_deseq2_SE <- function(.data,

my_contrasts = .contrasts

# If no assay is specified take first
my_assay = ifelse(
quo_is_symbol(.abundance),
quo_name(.abundance),
.data |>
assayNames() |>
extract2(1)
)

deseq2_object =
.data %>%


# DESeq2
DESeq2::DESeqDataSet( design = .formula) %>%
DESeq2::DESeqDataSetFromMatrix(
countData = .data |> assay(my_assay),
colData = colData(.data),
design = .formula
) %>%
DESeq2::DESeq(...)

# Return
Expand Down
29 changes: 25 additions & 4 deletions R/methods_SE.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1130,6 +1130,7 @@ setMethod(
#' @importFrom rlang inform
.test_differential_abundance_se = function(.data,
.formula,
.abundance = NULL,
contrasts = NULL,
method = "edgeR_quasi_likelihood",
test_above_log2_fold_change = NULL,
Expand All @@ -1139,6 +1140,8 @@ setMethod(
...)
{

.abundance = enquo(.abundance)

# Fix NOTEs
. = NULL

Expand Down Expand Up @@ -1179,8 +1182,9 @@ such as batch effects (if applicable) in the formula.
get_differential_transcript_abundance_bulk_SE(
.,
.formula,
.abundance = !!.abundance,
.contrasts = contrasts,
colData(.data),
sample_annotation = colData(.data),
method = method,
test_above_log2_fold_change = test_above_log2_fold_change,
scaling_method = scaling_method,
Expand All @@ -1193,8 +1197,9 @@ such as batch effects (if applicable) in the formula.
grepl("voom", method) ~ get_differential_transcript_abundance_bulk_voom_SE(
.,
.formula,
.abundance = !!.abundance,
.contrasts = contrasts,
colData(.data),
sample_annotation = colData(.data),
method = method,
test_above_log2_fold_change = test_above_log2_fold_change,
scaling_method = scaling_method,
Expand All @@ -1207,15 +1212,31 @@ such as batch effects (if applicable) in the formula.
tolower(method)=="deseq2" ~ get_differential_transcript_abundance_deseq2_SE(
.,
.formula,
.abundance = !!.abundance,
.contrasts = contrasts,
method = method,
test_above_log2_fold_change = test_above_log2_fold_change,
test_above_log2_fold_change = test_above_log2_fold_change,
scaling_method = scaling_method,
omit_contrast_in_colnames = omit_contrast_in_colnames,
prefix = prefix,
...
...
),

# glmmseq
tolower(method) %in% c("glmmseq_lme4", "glmmseq_glmmTMB") ~ get_differential_transcript_abundance_glmmSeq_SE(
.,
.formula,
.abundance = !!.abundance,
.contrasts = contrasts,
sample_annotation = colData(.data),
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,
...
),

# Else error
TRUE ~ stop("tidybulk says: the only methods supported at the moment are \"edgeR_quasi_likelihood\" (i.e., QLF), \"edgeR_likelihood_ratio\" (i.e., LRT), \"limma_voom\", \"limma_voom_sample_weights\", \"DESeq2\"")
)
Expand Down
Loading

0 comments on commit 136461b

Please sign in to comment.