##
## stats.R
## Statistics functions
##
## Written by Marta R. Hidalgo, marta.hidalgo@outlook.es
##
## Code style by Hadley Wickham (http://r-pkgs.had.co.nz/style.html)
## https://www.bioconductor.org/developers/how-to/coding-style/
##
#' Normalize expression data from a SummarizedExperiment or matrix
#' to be used in \code{hipathia}
#'
#' Transforms the rank of the SummarizedExperiment or matrix of gene expression
#' to [0,1] in order
#' to be processed by \code{hipathia}. The transformation may be performed
#' in two different ways. If \code{percentil = FALSE}, the transformation
#' is a re-scaling of the rank of the matrix. If \code{percentil = TRUE},
#' the transformation is performed assigning to each cell its percentil in
#' the corresponding distribution. This option is recommended for
#' distributions with very long tails.
#'
#' This transformation may be applied either to the whole matrix
#' (by setting \code{by_gene = FALSE}), which we strongly recommend, or to
#' each of the rows (by setting \code{by_gene = TRUE}), allowing each gene
#' to have its own scale.
#'
#' A previous quantiles normalization may be applied by setting
#' \code{by_quantiles = TRUE}. This is recommended for noisy data.
#'
#' For distributions with extreme outlayer values, a percentil \code{p}
#' may be given to the parameter \code{truncation_percentil}. When provided,
#' values beyond percentil p are truncated to the value of percentil p, and
#' values beyond 1-p are truncated to percentil 1-p. This step is performed
#' before any other tranformation. By default no truncation is performed.
#'
#' @param data Either a SummarizedExperiment or a matrix of gene expression.
#' @param sel_assay Character or integer, indicating the assay to be normalized
#' in the SummarizedExperiment. Default is 1.
#' @param by_quantiles Boolean, whether to normalize the data by quantiles.
#' Default is FALSE.
#' @param by_gene Boolean, whether to transform the rank of each row of the
#' matrix to [0,1]. Default is FALSE.
#' @param percentil Boolean, whether to take as value the percentil of each
#' sample in the corresponding distribution.
#' @param truncation_percentil Real number p in [0,1]. When provided, values
#' beyond percentil p are truncated to the value of percentil p, and values
#' beyond 1-p are truncated to percentil 1-p. By default no truncation
#' is performed.
#'
#' @return Matrix of gene expression whose values are in [0,1].
#'
#' @examples data("brca_data")
#' trans_data <- translate_data(brca_data, "hsa")
#' exp_data <- normalize_data(trans_data)
#' exp_data <- normalize_data(trans_data, by_quantiles = TRUE,
#' truncation_percentil=0.95)
#'
#' @export
#' @import preprocessCore
#' @import SummarizedExperiment
#' @importFrom stats quantile
#' @importFrom stats ecdf
#' @importFrom methods is
#'
normalize_data <- function(data, sel_assay = 1, by_quantiles = FALSE,
by_gene = FALSE, percentil = FALSE,
truncation_percentil = NULL){
if(is(data, "SummarizedExperiment")){
se_flag <- TRUE
mat <- assay(data, sel_assay)
}else if(is(data, "matrix")){
se_flag <- FALSE
mat <- data
}else{
stop("Only SummarizedExperiment or matrix classes accepted as data")
}
norm_mat <- normalize_matrix(mat, by_quantiles = by_quantiles,
by_gene = by_gene, percentil = percentil,
truncation_percentil=truncation_percentil)
if(se_flag == TRUE){
norm_data <- SummarizedExperiment(list(norm = norm_mat),
colData = colData(data))
}else{
norm_data <- norm_mat
}
return(norm_data)
}
#' @import preprocessCore
#' @importFrom stats quantile
#' @importFrom stats ecdf
#' @importFrom DelayedArray rowMins
#' @importFrom DelayedArray rowMaxs
normalize_matrix <- function(mat, sel_assay = 1, by_quantiles = FALSE,
by_gene = FALSE,
percentil = FALSE, truncation_percentil = NULL){
# Normalize data matrix
norm_data <- data.matrix(mat)
if(!is.null(truncation_percentil)){
if( truncation_percentil >= 0 & truncation_percentil <= 1){
# Guarantees that truncation_percentil is in [0.5,1]
if(truncation_percentil < (1-truncation_percentil))
truncation_percentil <- 1-truncation_percentil
# Truncates by the percentil
norm_data <- t(apply(norm_data, 1, function(x){
quan_inf <- stats::quantile(x, 1 - truncation_percentil,
na.rm = TRUE)
x[x < quan_inf] <- quan_inf
quan_sup <- stats::quantile(x, truncation_percentil,
na.rm = TRUE)
x[x > quan_sup] <- quan_sup
x
}))
}
else{
stop("Parameter truncation_percentil must be in [0,1]")
}
}
if(by_quantiles == TRUE){
norm_data <- preprocessCore::normalize.quantiles(norm_data)
}
if(by_gene == TRUE){
if(percentil == TRUE){
norm_data <- t(apply(norm_data, 1, function(x){stats::ecdf(x)(x)}))
}else{
min <- rowMins(norm_data, na.rm = TRUE)
max <- rowMaxs(norm_data, na.rm = TRUE)
norm_data <- (norm_data - min)/(max - min)
}
} else {
if(percentil == TRUE){
emp <- stats::ecdf(norm_data)
norm_data <- t(apply(norm_data, 1, emp))
}else{
norm_data <- (norm_data - min(norm_data, na.rm = TRUE))/
(max(norm_data, na.rm = TRUE) - min(norm_data, na.rm = TRUE))
}
}
colnames(norm_data) <- colnames(mat)
rownames(norm_data) <- rownames(mat)
return(norm_data)
}
#'@importFrom stats p.adjust
do_anova_test <- function(data, group, adjust = TRUE){
tests <- apply(data, 1, anova_test_fun, group = group)
out <- do.call("rbind", lapply(tests, function(x){
c(summary(x$fit)[[1]][[1, "Pr(>F)"]])
}))
out[is.na(out)] <- 1
if(adjust == TRUE){
out <- cbind(out, stats::p.adjust(out[,1], method = "fdr"))
} else {
out <- cbind(out, out[,1])
}
out <- as.data.frame(out, stringsAsFactors = FALSE)
colnames(out) <- c("p.value", "adj.p.value")
return(out)
}
#' @importFrom stats aov
#' @importFrom stats TukeyHSD
anova_test_fun <- function(values, group, verbose = FALSE){
group <- factor(group)
fit <- stats::aov(values ~ group)
tuckey <- stats::TukeyHSD(fit)
if(verbose==TRUE){
summary(fit)
}
return(list(fit = fit,tuckey = tuckey))
}
#'@importFrom stats p.adjust
do_cor <- function(sel_vals, design, adjust = TRUE){
data <- sel_vals[,design$sample]
testData <- do.call("rbind",apply(data, 1, cor_test_fun, design$value))
if(adjust==TRUE){
fdrData <- stats::p.adjust(testData[,1], method = "fdr")
} else {
fdrData <- testData[,1]
}
data2 <- data.frame(testData[,seq_len(3)], fdrData,
stringsAsFactors = FALSE)
colnames(data2) <- c("p.value", "UP/DOWN", "correlation", "FDRp.value")
data2[data2$statistic>0,"UP/DOWN"] <- "UP"
data2[data2$statistic<0,"UP/DOWN"] <- "DOWN"
return(data2)
}
#'@importFrom stats cor.test
cor_test_fun <- function(x, values){
r <- try(stats::cor.test(x = as.numeric(x), y = as.numeric(values)))
result <- cor_data_frame(r)
return(result)
}
cor_data_frame <- function(wilcox){
if (is(wilcox, "try-error") | is.na(wilcox$p.value)){
pvalue <- 1
class <- "0"
esti <- 0
} else{
pvalue <- wilcox$p.value
esti <- wilcox$estimate[1]
if (esti < 0){
class <- "DOWN" ## regarding DISEASE
} else if (esti > 0){
class <- "UP" ## regarding DISEASE
} else if (esti == 0){
if (wilcox$conf.int[1] == 0){
class <- "UP"
} else if (wilcox$conf.int[2] == 0){
class <- "DOWN"
} else{
class <- 0
}
}
}
result <- data.frame(pvalue, class, esti, stringsAsFactors = FALSE)
return(result)
}
#' Apply Wilcoxon test
#'
#' Performs a Wilcoxon test for the values in \code{sel_vals} comparing
#' conditions \code{g1} and \code{g2}
#'
#' @param data Either a SummarizedExperiment object or a matrix, containing the
#' values. Columns represent samples.
#' @param group Either a character indicating the name of the column in colData
#' including the classes to compare, or a character vector with the class to
#' which each sample belongs.
#' Samples must be ordered as in \code{data}
#' @param g1 String, label of the first group to be compared
#' @param g2 String, label of the second group to be compared
#' @param paired Boolean, whether the samples to be compared are paired.
#' If TRUE, function \code{wilcoxsign_test} from package \code{coin} is
#' used. If FALSE, function \code{wilcox.test} from package \code{stats}
#' is used.
#' @param adjust Boolean, whether to adjust the p.value with
#' Benjamini-Hochberg FDR method
#' @param sel_assay Character or integer, indicating the assay to be normalized
#' in the SummarizedExperiment. Default is 1.
#' @param order Boolean, whether to order the results table by the
#' \code{FDRp.value} column. Default is FALSE.
#'
#' @return Dataframe with the result of the comparison
#'
#' @examples
#' data(path_vals)
#' data(brca_design)
#' sample_group <- brca_design[colnames(path_vals),"group"]
#' comp <- do_wilcoxon(path_vals, sample_group, g1 = "Tumor", g2 = "Normal")
#'
#' @export
#' @import SummarizedExperiment
#' @importFrom methods is
#'
do_wilcoxon <- function(data, group, g1, g2, paired = FALSE,
adjust = TRUE, sel_assay = 1, order = FALSE){
if(is(data, "SummarizedExperiment")){
se_flag <- TRUE
if(is(group, "character") & length(group) == 1)
if(group %in% colnames(colData(data))){
group <- colData(data)[[group]]
}else{
stop("Group variable must be a column in colData(data)")
}
vals <- assay(data, sel_assay)
}else if(is(data, "matrix")){
se_flag <- FALSE
vals <- data
}else{
stop("Only SummarizedExperiment or matrix classes accepted as data")
}
g1_indexes <- which(group == g1)
g2_indexes <- which(group == g2)
stat_vals <- suppressWarnings(
calculate_wilcox_test(vals, g2_indexes, g1_indexes, paired = paired,
adjust = adjust))
if(se_flag == TRUE && "feat.name" %in% colnames(rowData(data)))
stat_vals <- cbind(name = rowData(data)[["feat.name"]], stat_vals)
if(order == TRUE)
stat_vals <- stat_vals[order(stat_vals$FDRp.value, decreasing = FALSE),]
return(stat_vals)
}
#'@importFrom stats p.adjust
calculate_wilcox_test <- function(data, control, disease, paired, adjust=TRUE){
if(paired == TRUE){
dat <- apply(data, 1, wilcoxsign_test_fun, control, disease)
testData <- do.call("rbind", dat)
if(adjust == TRUE){
fdrData <- stats::p.adjust(testData[,1], method = "fdr")
} else {
fdrData <- testData[,1]
}
data2 <- data.frame(testData[,c(2,3,1)], fdrData,
stringsAsFactors = FALSE)
}else{
testData <- do.call("rbind",
apply(data, 1, wilcox_test_fun,
control, disease, paired))
if(adjust == TRUE){
fdrData <- stats::p.adjust(testData[,1], method = "fdr")
} else {
fdrData <- testData[,1]
}
# Standardize statistic
lc <- length(control)
ld <- length(disease)
m <- lc*ld/2
sigma <- sqrt(lc * ld * (lc + ld + 1)/12)
z <- (testData[,3] - m)/sigma
data2 <- data.frame(testData[,2], z, testData[,1], fdrData,
stringsAsFactors = FALSE)
rownames(data2) <- rownames(testData)
need0 <- which(data2$pvalue == 1 &
data2$class == "0" &
data2$fdrData == 1)
data2[need0, "z"] <- 0
}
colnames(data2) <- c("UP/DOWN", "statistic", "p.value", "FDRp.value")
data2[data2$statistic > 0,"UP/DOWN"] <- "UP"
data2[data2$statistic < 0,"UP/DOWN"] <- "DOWN"
return(data2)
}
wilcoxsign_test_fun <- function(x, control, disease){
r <- try(coin::wilcoxsign_test(
as.numeric(x[disease]) ~ as.numeric(x[control]), showWarnings = FALSE))
if (is(r, "try-error")){
pvalue <- 1
class <- "0"
stat <- 0
} else{
pvalue <- coin::pvalue(r)
stat <- coin::statistic(r, "standardized")[1]
if(stat > 0){
class <- "UP"
}else if (stat < 0){
class <- "DOWN"
}else{
class <- "0"
}
}
result <- data.frame(pvalue = pvalue,
class = class,
stat = stat,
stringsAsFactors = FALSE)
return(result)
}
#'@importFrom stats wilcox.test
wilcox_test_fun <- function(x, control, disease, paired){
r <- try(stats::wilcox.test(x = as.numeric(x[disease]),
y = as.numeric(x[control]),
conf.int = TRUE,
alternative = "two.sided",
paired = paired),
silent = TRUE)
result <- wilcox_data_frame(r)
return(result)
}
wilcox_data_frame <- function(wilcox){
if (is(wilcox, "try-error")){
pvalue <- 1
class <- "0"
stat <- 0
} else{
pvalue <- wilcox$p.value
esti <- wilcox$estimate
stat <- wilcox$statistic[[1]]
if (esti < 0){
class <- "DOWN" ## regarding DISEASE
} else if (esti > 0){
class <- "UP" ## regarding DISEASE
} else if (esti == 0){
if (wilcox$conf.int[1] == 0){
class <- "UP"
} else if (wilcox$conf.int[2] == 0){
class <- "DOWN"
} else{
class <- 0
}
}
}
result <- data.frame(pvalue, class, stat,stringsAsFactors=FALSE)
return(result)
}
#'
#' Performs a Principal Components Analysis
#'
#' @param data SummarizedExperiment or matrix of values to be analyzed. Samples
#' must be represented in the columns.
#' @param sel_assay Character or integer, indicating the assay to be normalized
#' in the SummarizedExperiment. Default is 1.
#' @param cor A logical value indicating whether the calculation should use
#' the correlation matrix or the covariance matrix. (The correlation matrix
#' can only be used if there are no constant variables.)
#'
#' @return \code{do_pca} returns a list with class \code{princomp}.
#'
#' @examples
#' data(path_vals)
#' pca_model <- do_pca(path_vals[seq_len(ncol(path_vals)),])
#'
#' @export
#' @importFrom stats princomp
#' @importFrom methods is
#'
do_pca <- function(data, sel_assay = 1, cor = FALSE){
if(is(data, "SummarizedExperiment")){
data <- assay(data, sel_assay)
}else if(is(data, "matrix")){
data <- data
}else{
stop("Only SummarizedExperiment or matrix classes accepted as data")
}
fit <- stats::princomp(t(data), cor = cor)
fit$var <- fit$sdev^2
fit$explain_var <- fit$var/sum(fit$var)
fit$acum_explain_var <- cumsum(fit$explain_var)
return(fit)
}
do_limma <- function(data, groups, expdes, g2 = NULL, sel_assay = 1,
order = FALSE){
if(is(data, "SummarizedExperiment")){
se_flag <- TRUE
if(is(groups, "character") & length(groups) == 1)
if(groups %in% colnames(colData(data))){
groups <- colData(data)[[groups]]
}else{
stop("Group variable must be a column in colData(data)")
}
vals <- assay(data, sel_assay)
}else if(is(data, "matrix")){
se_flag <- FALSE
vals <- data
}else{
stop("Only SummarizedExperiment or matrix classes accepted as data")
}
if(!is.null(g2))
expdes <- paste(expdes, "-", g2)
# Make contrasts
design <- stats::model.matrix(~0 + factor(groups))
colnames(design) <- levels(factor(groups))
rownames(design) <- colnames(vals)
cont_matrix <- limma::makeContrasts(contrasts = expdes, levels = design)
# Remove nodes with 0 variance
vars <- apply(vals, 1, var)
novar <- vals[-which(vars == 0),]
# Do analysis
fit <- limma::lmFit(novar, design)
fit1 <- limma::contrasts.fit(fit, cont_matrix)
fit2 <- limma::eBayes(fit1)
tt <- limma::topTable(fit2, coef = 1, number = "all", sort.by = "none")
# Add nodes with 0 variance to results table
ttnovar <- data.frame(path = names(vars[vars == 0]), logFC = 0, AveExpr = 1,
t = -8, P.Value = 1, adj.P.Val = 1, B = -10)
rownames(ttnovar) <- ttnovar$path
ttnovar <- ttnovar[,-1]
tt <- rbind(tt, ttnovar)
tt <- tt[rownames(vals),]
# Create results data frame
updown <- c("UP", "DOWN", "UP")
names(updown) <- c("1", "-1", "0")
ud <- updown[as.character(sign(tt$logFC))]
comp <- data.frame(ud, tt$t, tt$P.Value, tt$adj.P.Val,
stringsAsFactors = FALSE)
colnames(comp) <- c("UP/DOWN", "statistic", "p.value", "FDRp.value")
rownames(comp) <- rownames(tt)
if(se_flag == TRUE && "feat.name" %in% colnames(rowData(data)))
comp <- cbind(name = rowData(data)[["feat.name"]], comp)
if(order == TRUE)
comp <- comp[order(comp$p.value, decreasing = FALSE),]
return(comp)
}
#' Compute pathway summary
#'
#' Computes a summary of the results, summarizing the number and proportion
#' of up- and down-regulated subpathways in each pathway.
#'
#' @param comp Comparison data frame as returned by the \code{do_wilcoxon}
#' function.
#' @param metaginfo Pathways object
#' @param conf Level of significance of the comparison for the adjusted
#' p-value. Default is 0.05.
#'
#' @return Table with the summarized information for each of the pathways.
#' Rows are the analized pathways. Columns are:
#' * \code{num_total_paths} Number of total subpathways in which each pathway
#' is decomposed.
#' * \code{num_significant_paths} Number of significant subpathways in the
#' provided comparison.
#' * \code{percent_significant_paths} Percentage of significant subpathways
#' from the total number of subpathways in a pathway.
#' * \code{num_up_paths} Number of significant up-regulated subpathways in the
#' provided comparison.
#' * \code{percent_up_paths} Percentage of significant up-regulated subpathways
#' from the total number of subpathways in a pathway.
#' * \code{num_down_paths} Number of significant down-regulated subpathways in
#' the provided comparison.
#' * \code{percent_down_paths} Percentage of significant down-regulated
#' subpathways from the total number of subpathways in a pathway.
#'
#' @examples
#' data(comp)
#' pathways <- load_pathways(species = "hsa", pathways_list = c("hsa03320",
#' "hsa04012"))
#' get_pathways_summary(comp, pathways)
#'
#' @export
#'
get_pathways_summary <- function(comp, metaginfo, conf = 0.05){
comp$pathways <- sapply(strsplit(rownames(comp), split = "-"), "[[", 2)
local_paths <- unique(metaginfo$all.labelids[,c("path.id", "path.name")])
id_pathways <- local_paths[,"path.id"]
name_pathways <- local_paths[,"path.name"]
summ <- lapply(id_pathways, function(pathway){
minicomp <- comp[comp$pathways == pathway,]
num_total_paths <- nrow(minicomp)
num_sig_paths <- sum(minicomp$FDRp.value < conf)
percent_sig_paths <- round(num_sig_paths/nrow(minicomp)*100, digits = 2)
is_up_path <- minicomp$FDRp.value < conf & minicomp$`UP/DOWN` == "UP"
num_up_paths <- sum(is_up_path)
percent_up_paths <- round(num_up_paths/nrow(minicomp) * 100, digits = 2)
is_down_path <- minicomp$FDRp.value < conf &
minicomp$`UP/DOWN` == "DOWN"
num_down_paths <- sum(is_down_path)
percent_down_paths <- round(num_down_paths/nrow(minicomp) * 100,
digits = 2)
data.frame(num_total_paths = num_total_paths,
num_significant_paths = num_sig_paths,
percent_significant_paths = percent_sig_paths,
num_up_paths = num_up_paths,
percent_up_paths = percent_up_paths,
num_down_paths = num_down_paths,
percent_down_paths = percent_down_paths)
})
summ <- do.call("rbind", summ)
summ <- cbind(id_pathways, summ)
rownames(summ) <- name_pathways
summ <- summ[order(summ$percent_significant_paths, decreasing = TRUE),]
return(summ)
}
#' Computes pathway significance
#'
#' Performs a test for each pathway checking if the number of significant
#' paths is significant, compared to not having any of the paths as significant.
#'
#' @param comp Comparison data frame as returned by the \code{do_wilcoxon}
#' function.
#'
#' @return
#' Table with the names of the pathways and their p-value for the Fisher test
#' comparing the proportion of significant subpaths vs. 0.
#'
#' @examples
#' data(comp)
#' top_pathways(comp)
#'
#' @export
#'
top_pathways <- function(comp){
if("name" %in% colnames(comp)){
path_names <- as.character(comp$name)
comp$pathways <- sapply(strsplit(path_names, split = ":"), "[[", 1)
}else{
path_names <- rownames(comp)
comp$pathways <- sapply(strsplit(path_names, split = "-"), "[[", 2)
}
pathways <- unique(comp$pathways)
tests <- do.call(rbind, lapply(pathways, function(path) {
t1 <- table(comp[comp$pathways == path,"FDRp.value"] < 0.05)
t2 <- c(sum(t1), 0)
ft <- fisher.test(rbind(t1, t2))
data.frame(pathway = path, pval = ft$p.value, stringsAsFactors = FALSE)
}))
tests <- tests[order(tests$pval),]
return(tests)
}
