Skip to content

Commit

Permalink
Add PomaBatch function
Browse files Browse the repository at this point in the history
  • Loading branch information
pcastellanoescuder committed Jan 19, 2024
1 parent 89be7f6 commit cf351da
Show file tree
Hide file tree
Showing 12 changed files with 195 additions and 158 deletions.
2 changes: 1 addition & 1 deletion DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
Package: POMA
Title: Tools for Omics Data Analysis
Version: 1.13.10
Version: 1.13.11
Authors@R:
c(person(given = "Pol",
family = "Castellano-Escuder",
Expand Down
1 change: 1 addition & 0 deletions NAMESPACE
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
# Generated by roxygen2: do not edit by hand

export("%>%")
export(PomaBatch)
export(PomaBoxplots)
export(PomaClust)
export(PomaCorr)
Expand Down
2 changes: 1 addition & 1 deletion NEWS.md
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
# POMA 1.13.10
# POMA 1.13.11

* New POMA theme and colorblind-friendly palette
* Available sample normalization (sum and quantile)
Expand Down
56 changes: 56 additions & 0 deletions R/PomaBatch.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,56 @@

#' Batch Correction
#'
#' @description `PomaBatch` performs batch correction on a `SummarizedExperiment` object given a batch factor variable.
#'
#' @param data A `SummarizedExperiment` object.
#' @param batch Character. The name of the column in `colData` that contains the batch information.
#' @param mod Character vector. Indicates the names of `colData` columns to be included as covariates. Default is NULL (no covariates).
#'
#' @export
#'
#' @return A `SummarizedExperiment` object with batch-corrected data.
#' @references Leek JT, Johnson WE, Parker HS, Fertig EJ, Jaffe AE, Zhang Y, Storey JD, Torres LC (2023). sva: Surrogate Variable Analysis. doi:10.18129/B9.bioc.sva <https://doi.org/10.18129/B9.bioc.sva>
#' @author Pol Castellano-Escuder
#'
#' @importFrom magrittr %>%
#'
#' @examples
#' data("st000284")
#'
#' st000284 %>%
#' PomaImpute(method = "knn") %>%
#' PomaBatch(batch = "gender")
PomaBatch <- function(data,
batch,
mod = NULL) {

if (!is(data, "SummarizedExperiment")) {
stop("data is not a SummarizedExperiment object. \nSee POMA::PomaCreateObject or SummarizedExperiment::SummarizedExperiment")
}
if (!batch %in% names(SummarizedExperiment::colData(data))) {
stop("Specified batch column not found in colData")
}

batch_info <- as.factor(SummarizedExperiment::colData(data)[, batch])
to_batch <- SummarizedExperiment::assay(data)

if (!is.null(mod)) {
covariates <- SummarizedExperiment::colData(data) %>%
as.data.frame() %>%
dplyr::select_at(dplyr::vars(dplyr::matches(mod)))

mod_matrix <- stats::model.matrix(~ 0 + ., data = covariates)
} else {
mod_matrix <- NULL
}

corrected_data <- sva::ComBat(dat = to_batch, batch = batch_info, mod = mod_matrix)

data <- SummarizedExperiment::SummarizedExperiment(assays = corrected_data,
colData = SummarizedExperiment::colData(data))

if (validObject(data))
return(data)
}

109 changes: 29 additions & 80 deletions R/PomaVolcano.R
Original file line number Diff line number Diff line change
@@ -1,106 +1,55 @@

#' Volcano Plot
#'
#' @description `PomaVolcano` generates a volcano plot. Data should not contain negative values.
#' @description `PomaVolcano` creates a volcano plot from a given dataset. This function is designed to visualize the statistical significance (p-value) against the magnitude of change (log2 fold change) for features.
#'
#' @param data A `SummarizedExperiment` object.
#' @param method Character. Indicates the statistical method to use. Options are "ttest", "mann", "limma", and "DESeq".
#' @param pval Character. Indicates the pvalue type to generate the volcano plot. Options are: "raw" and "adjusted".
#' @param pval_cutoff Numeric. Indicated the pvalue cutoff (horizontal line).
#' @param adjust Character, Indicates the multiple correction method. Options are: "fdr", "holm", "hochberg", "hommel", "bonferroni", "BH" and "BY".
#' @param log2fc_cutoff Numeric. Indicates the log2 fold change cutoff (vertical lines).
#' @param labels Logical. Indicates if significant labels should be plotted. Defaul is FALSE.
#' @param paired Logical. Indicates if the data is paired or not. Default is FALSE.
#' @param var_equal Logical. Indicates if the data variances are assumed to be equal or not. Default is FALSE.
#'
#' @export
#' @param data A data frame with at least three columns: feature names, statistical significance values, and magnitude of change values. These should be the first three columns of the data, in this exact order.
#' @param pval_cutoff Numeric. Specifies the p-value threshold for significance in the volcano plot. The default is set to 0.05. This parameter determines the horizontal line in the plot indicating the threshold for statistical significance.
#' @param log2fc_cutoff Numeric. Specifies the log2 fold change cutoff for the volcano plot. If not provided, the cutoff is set to the 75th percentile of the absolute log2 fold changes in the data. This parameter determines the vertical lines in the plot indicating the magnitude of change threshold.
#' @param labels Logical. Indicates whether to plot labels for significant features.
#' @param x_label Character. Custom label for the x-axis.
#' @param y_label Character. Custom label for the y-axis.
#'
#' @return A `ggplot` object.
#' @export
#'
#' @return A `ggplot` object representing the volcano plot.
#' @author Pol Castellano-Escuder
#'
#' @importFrom magrittr %>%
#'
#' @examples
#' data("st000336")
#' @importFrom magrittr %>%
#'
#' @examples
#' st000336 %>%
#' PomaImpute() %>%
#' PomaVolcano()
#'
#' st000336 %>%
#' PomaImpute() %>%
#' PomaVolcano(labels = TRUE)
#'
#' st000336 %>%
#' PomaImpute() %>%
#' PomaVolcano(labels = TRUE, pval = "adjusted")
#' PomaImpute() %>%
#' PomaUnivariate() %>%
#' magrittr::extract2("result") %>%
#' dplyr::select(feature, fold_change, pvalue) %>%
#' PomaVolcano()
PomaVolcano <- function(data,
method = "ttest",
pval = "raw",
pval_cutoff = 0.05,
adjust = "fdr",
log2fc_cutoff = NULL,
labels = FALSE,
paired = FALSE,
var_equal = FALSE) {
x_label = "log2 (Fold Change)",
y_label = "-log10 (P-value)") {

if(!is(data, "SummarizedExperiment")){
stop("data is not a SummarizedExperiment object. \nSee POMA::PomaCreateObject or SummarizedExperiment::SummarizedExperiment")
}
if (length(table(SummarizedExperiment::colData(data)[,1])) != 2) {
stop("Grouping factor must have exactly 2 levels (first column of the metadata file)")
}
if (!(pval %in% c("raw", "adjusted"))) {
stop('Incorrect value for pval argument. Options are: "raw" and "adjusted"')
}
if (!(adjust %in% c("fdr", "holm", "hochberg", "hommel", "bonferroni", "BH", "BY"))) {
stop("Incorrect value for adjust argument.")
}

if (method == "ttest") {
volcano_res <- data %>%
PomaUnivariate(method = "ttest", adjust = adjust, paired = paired, var_equal = var_equal) %>%
magrittr::extract2("result") %>%
dplyr::mutate(logFC = log2(fold_change))
}
else if (method == "mann") {
volcano_res <- data %>%
PomaUnivariate(method = "mann", adjust = adjust, paired = paired, var_equal = var_equal) %>%
magrittr::extract2("result") %>%
dplyr::mutate(logFC = log2(fold_change))
}
else if (method == "limma") {
contrast <- paste0(rev(names(table(SummarizedExperiment::colData(data)[,1]))), collapse = "-")

volcano_res <- data %>%
PomaLimma(adjust = adjust, contrast = contrast)
}
else if (method == "DESeq") {
volcano_res <- data %>%
PomaDESeq(adjust = adjust)
}

if (pval == "raw") {
volcano_res <- volcano_res %>%
dplyr::select(feature, logFC, pvalue = pvalue)
} else {
volcano_res <- volcano_res %>%
dplyr::select(feature, logFC, pvalue = adj_pvalue)
}

to_volcano <- data %>%
as.data.frame() %>%
dplyr::select(feature = 1, logFC = 2, pvalue = 3)

if (is.null(log2fc_cutoff)) {
log2fc_cutoff <- quantile(abs(volcano_res$logFC), 0.75, na.rm = TRUE)
log2fc_cutoff <- quantile(abs(to_volcano$logFC), 0.75, na.rm = TRUE)
}

plot_complete <- volcano_res %>%
plot_complete <- to_volcano %>%
ggplot2::ggplot(ggplot2::aes(logFC, -log10(pvalue), label = feature)) +
ggplot2::geom_point(fill = "gray", size = 3, pch = 21, alpha = 0.8) +
ggplot2::geom_point(data = volcano_res[volcano_res$pvalue < pval_cutoff & abs(volcano_res$logFC) > log2fc_cutoff ,], fill = "red", size = 3, pch = 21,) +
{if(labels) ggrepel::geom_label_repel(data = volcano_res[volcano_res$pvalue < pval_cutoff & abs(volcano_res$logFC) > log2fc_cutoff ,], color = "black", size = 4)} +
ggplot2::geom_point(data = to_volcano[to_volcano$pvalue < pval_cutoff & abs(to_volcano$logFC) > log2fc_cutoff ,], fill = "red", size = 3, pch = 21,) +
{if(labels) ggrepel::geom_label_repel(data = to_volcano[to_volcano$pvalue < pval_cutoff & abs(to_volcano$logFC) > log2fc_cutoff ,], color = "black", size = 4)} +
ggplot2::geom_vline(xintercept = c(-log2fc_cutoff, log2fc_cutoff), linetype = "dashed", color = "orange") +
ggplot2::geom_hline(yintercept = -log10(pval_cutoff), linetype = "dashed", color = "orange") +
ggplot2::labs(x = "log2 (Fold Change)",
y = ifelse(pval == "raw", "-log10 (P-value)", "-log10 (Adjusted P-value)")) +
ggplot2::labs(x = x_label,
y = y_label) +
theme_poma()

return(plot_complete)
Expand Down
2 changes: 1 addition & 1 deletion README.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ output: github_document
| _BioC_ branch | Status | Version | Dependencies | Rank |
|- |- |- |- |- |
| [Release](http://bioconductor.org/packages/release/bioc/html/POMA.html) | [![Bioc release status](https://bioconductor.org/shields/build/release/bioc/POMA.svg)](https://bioconductor.org/checkResults/release/bioc-LATEST/POMA/) | [![BioC released version](https://img.shields.io/badge/release%20version-1.6.0-blue.svg)](https://www.bioconductor.org/packages/POMA) | [![Dependencies](http://bioconductor.org/shields/dependencies/release/POMA.svg)](http://bioconductor.org/packages/release/bioc/html/POMA.html#since) | [![Rank](http://www.bioconductor.org/shields/downloads/release/POMA.svg)](https://bioconductor.org/packages/stats/bioc/POMA) |
| [Devel](http://bioconductor.org/packages/devel/bioc/html/POMA.html) | [![Bioc devel status](https://bioconductor.org/shields/build/devel/bioc/POMA.svg)](https://bioconductor.org/checkResults/devel/bioc-LATEST/POMA/) | [![BioC devel version](https://img.shields.io/badge/devel%20version-1.13.10-blue.svg)](https://bioconductor.org/packages/devel/bioc/html/POMA.html) | [![Dependencies](http://bioconductor.org/shields/dependencies/devel/POMA.svg)](http://bioconductor.org/packages/devel/bioc/html/POMA.html#since) | [![Rank](http://www.bioconductor.org/shields/downloads/devel/POMA.svg)](https://bioconductor.org/packages/stats/bioc/POMA) |
| [Devel](http://bioconductor.org/packages/devel/bioc/html/POMA.html) | [![Bioc devel status](https://bioconductor.org/shields/build/devel/bioc/POMA.svg)](https://bioconductor.org/checkResults/devel/bioc-LATEST/POMA/) | [![BioC devel version](https://img.shields.io/badge/devel%20version-1.13.11-blue.svg)](https://bioconductor.org/packages/devel/bioc/html/POMA.html) | [![Dependencies](http://bioconductor.org/shields/dependencies/devel/POMA.svg)](http://bioconductor.org/packages/devel/bioc/html/POMA.html#since) | [![Rank](http://www.bioconductor.org/shields/downloads/devel/POMA.svg)](https://bioconductor.org/packages/stats/bioc/POMA) |

<!-- badges: end -->

Expand Down
2 changes: 1 addition & 1 deletion README.md
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/li
| *BioC* branch | Status | Version | Dependencies | Rank |
|-------------------------------------------------------------------------|---------------------------------------------------------------------------------------------------------------------------------------------------------|-----------------------------------------------------------------------------------------------------------------------------------------------------|------------------------------------------------------------------------------------------------------------------------------------------------------|------------------------------------------------------------------------------------------------------------------------------|
| [Release](http://bioconductor.org/packages/release/bioc/html/POMA.html) | [![Bioc release status](https://bioconductor.org/shields/build/release/bioc/POMA.svg)](https://bioconductor.org/checkResults/release/bioc-LATEST/POMA/) | [![BioC released version](https://img.shields.io/badge/release%20version-1.6.0-blue.svg)](https://www.bioconductor.org/packages/POMA) | [![Dependencies](http://bioconductor.org/shields/dependencies/release/POMA.svg)](http://bioconductor.org/packages/release/bioc/html/POMA.html#since) | [![Rank](http://www.bioconductor.org/shields/downloads/release/POMA.svg)](https://bioconductor.org/packages/stats/bioc/POMA) |
| [Devel](http://bioconductor.org/packages/devel/bioc/html/POMA.html) | [![Bioc devel status](https://bioconductor.org/shields/build/devel/bioc/POMA.svg)](https://bioconductor.org/checkResults/devel/bioc-LATEST/POMA/) | [![BioC devel version](https://img.shields.io/badge/devel%20version-1.13.10-blue.svg)](https://bioconductor.org/packages/devel/bioc/html/POMA.html) | [![Dependencies](http://bioconductor.org/shields/dependencies/devel/POMA.svg)](http://bioconductor.org/packages/devel/bioc/html/POMA.html#since) | [![Rank](http://www.bioconductor.org/shields/downloads/devel/POMA.svg)](https://bioconductor.org/packages/stats/bioc/POMA) |
| [Devel](http://bioconductor.org/packages/devel/bioc/html/POMA.html) | [![Bioc devel status](https://bioconductor.org/shields/build/devel/bioc/POMA.svg)](https://bioconductor.org/checkResults/devel/bioc-LATEST/POMA/) | [![BioC devel version](https://img.shields.io/badge/devel%20version-1.13.11-blue.svg)](https://bioconductor.org/packages/devel/bioc/html/POMA.html) | [![Dependencies](http://bioconductor.org/shields/dependencies/devel/POMA.svg)](http://bioconductor.org/packages/devel/bioc/html/POMA.html#since) | [![Rank](http://www.bioconductor.org/shields/downloads/devel/POMA.svg)](https://bioconductor.org/packages/stats/bioc/POMA) |

<!-- badges: end -->

Expand Down
34 changes: 34 additions & 0 deletions man/PomaBatch.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

46 changes: 15 additions & 31 deletions man/PomaVolcano.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

26 changes: 26 additions & 0 deletions tests/testthat/test-PomaBatch.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,26 @@

test_that("PomaBatch handles valid SummarizedExperiment objects", {
data <- create_mock_summarized_experiment()
SummarizedExperiment::colData(data)$batch <- factor(rep(c("Batch1", "Batch2"), each = ncol(data)/2))
corrected_data <- PomaBatch(data, batch = "batch")
expect_is(corrected_data, "SummarizedExperiment")
})

test_that("PomaBatch stops with non-SummarizedExperiment objects", {
data <- data.frame(matrix(runif(100), ncol = 10))
expect_error(PomaBatch(data, batch = "batch"))
})

test_that("PomaBatch stops if batch column is not found", {
data <- create_mock_summarized_experiment()
expect_error(PomaBatch(data, batch = "nonexistent_batch"))
})

test_that("PomaBatch handles additional covariates", {
data <- create_mock_summarized_experiment()
SummarizedExperiment::colData(data)$batch <- factor(rep(c("Batch1", "Batch2"), each = ncol(data)/2))
SummarizedExperiment::colData(data)$covariate <- rnorm(nrow(SummarizedExperiment::colData(data)))
corrected_data <- PomaBatch(data, batch = "batch", mod = c("covariate"))
expect_is(corrected_data, "SummarizedExperiment")
})

Loading

0 comments on commit cf351da

Please sign in to comment.