Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Allow custom assay for se de #274

Merged
merged 2 commits into from
Jul 6, 2023
Merged

Allow custom assay for se de #274

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
d600ef4
added .abundance argument
stemangiola Jul 6, 2023
35b4fb5
fix test
stemangiola Jul 6, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions NAMESPACE
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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