R/proprotionality.R

In abodein/timeOmics: Time-Course Multi-Omics data integration

#' Proportionality Distance
#' 
#' \code{proportionality} is a wrapper that compute proportionality distance for 
#' a clustering result (\code{pca}, \code{spca}, \code{pls}, \code{spls}, \code{block.pls}, \code{block.spls}).
#' and it performs a u-test to compare the median within a cluster to the median of the entire background set.
#' 
#' @param X an object of the class: \code{pca}, \code{spca}, \code{pls}, \code{spls}, \code{block.pls} or \code{block.spls}
#' 
#' @return 
#' Return a list containing the following components:
#'   \item{propr.distance}{Square matrix with proportionality distance between pairs of features}
#'   \item{propr.distance.w.cluster}{distance between pairs with cluster label}
#'   \item{pvalue}{Wilcoxon U-test p-value comparing the medians within clusters and with the entire background set}
#' 
#' 
#' @references 
#' Lovell, D., Pawlowsky-Glahn, V., Egozcue, J. J., Marguerat, S., Bähler, J. (2015). Proportionality: a valid alternative to correlation for relative data. PLoS Comput. Biol. 11, e1004075. doi: 10.1371/journal.pcbi.1004075
#' 
#' Quinn, T. P., Richardson, M. F., Lovell, D., Crowley, T. M. (2017). propr: an r-package for identifying proportionally abundant features using compositional data analysis. Sci. Rep. 7, 16252. doi: 10.1038/s41598-017-16520-0
#' 
#' @examples 
#' demo <- suppressWarnings(get_demo_cluster())
#' 
#' # pca
#' X <- demo$pca
#' propr.res <- proportionality(X)
#' plot(propr.res)
#' 
#' # pls
#' X <- demo$spls
#' propr.res <- proportionality(X)
#' plot(propr.res)
#' 
#' # block.pls
#' X <- demo$block.spls
#' propr.res <- proportionality(X)
#' plot(propr.res)
#' 
#' @importFrom dplyr mutate filter rename left_join
#' @importFrom magrittr %>%
#' @importFrom tibble rownames_to_column
#' @importFrom tidyr pivot_longer
#' 
#' @export
proportionality <- function(X){
    #stopifnot(is(X, c("pca", "spca", "mixo_pls", "mixo_spls", "block.pls", "block.spls")))
    # modif bioc 3.16
    stopifnot(any(class(X) %in% c("pca", "spca", "mixo_pls", "mixo_spls", "block.pls", "block.spls")))
    # 1. get cluster
    cluster.info <- getCluster(X) %>%
        dplyr::select(molecule, cluster)
    
    # 2. extract data and add cluster
    
    # pca / spca
    if(any(class(X) %in% c("pca", "spca"))){
        # unscaling + positive value
        data <- unscale(X$X) %>% `+`(abs(min(.))) 
    } else if(any(class(X) %in% c("mixo_pls", "mixo_spls"))){
        # unscale X, unscale Y, cat
        data.X <- unscale(X$X) %>%
            `+`(abs(min(.)))
        data.Y <- unscale(X$Y) %>%
            `+`(abs(min(.)))
        data <- cbind(data.X, data.Y)
    } else if(any(class(X) %in% c("block.pls", "block.spls"))){
        # if(is.null(X$Y)){  ## no need: Y is passed to X
        data <- lapply(X$X, function(x) x %>%
                           unscale %>%
                           `+`(abs(min(.)))) %>%
            do.call(what="cbind")
        # } else {
        #     data.X <- lapply(X$X, function(x) x %>% unscale %>% `+`(abs(min(.)))) %>%
        #         do.call(what="cbind")
        #     data.Y <- unscale(X$Y) %>% `+`(abs(min(.)))
        #     data <- cbind(data.X, data.Y)
        # }
    }
    
    # 3. compute phi_s
    
    # update 16 march 2023: propr deprecated;
    # implementation of phi_s only based on Quinn paper
    # data.propr <- suppressMessages(as.data.frame(propr::propr(data, metric = "phs")@matrix))
    
    # gather, add cluster info, compute basic stats
    data.propr <- get_phi_s(data)
    data.propr.gather <- data.propr %>%
        tibble::rownames_to_column("feature1") %>%
        tidyr::pivot_longer(-feature1, names_to = "feature2", values_to = 'value') %>%
        dplyr::filter(feature1 %in% cluster.info$molecule) %>%
        dplyr::left_join(cluster.info, by = c("feature1"="molecule")) %>%
        dplyr::rename("cluster1"="cluster") %>%
        dplyr::filter(feature2 %in% cluster.info$molecule) %>%
        dplyr::left_join(cluster.info, by = c("feature2"="molecule")) %>%
        dplyr::rename("cluster2"="cluster") %>%
        dplyr::mutate(insideout = ifelse(cluster1 == cluster2, "inside", "outside"))
        
    # compute stat (median, u-test pval, adj.pval)
    res.stat <- stat_median(data.propr.gather) %>%
        na.omit()
    
    res <- list()
    res[["propr.distance"]] <- data.propr
    res[["propr.distance.w.cluster"]] <- data.propr.gather
    res[["pvalue"]] <- res.stat
    class(res) <- "proportionality"
    return(res)
}

#' @importFrom dplyr mutate filter pull
#' @importFrom purrr set_names
#' @importFrom magrittr %>%
stat_median <- function(res.phs.X){
    i = 1
    res.pval <- matrix(ncol = 4, nrow = 4) %>%
        as.data.frame() %>%
        purrr::set_names("cluster", "median_inside", "median_outside", "Pvalue")
    for(clu in unique(res.phs.X$cluster1)){
        inside <- res.phs.X %>%
            filter(cluster1 == clu) %>% 
            dplyr::filter(cluster2==clu) %>% 
            dplyr::pull(value)
        outside <- res.phs.X %>% 
            filter(cluster1 == clu) %>% 
            dplyr::filter(cluster2!=clu) %>% 
            dplyr::pull(value)
        
        #ttest.pval <- t.test(inside, outside)$p.value
        utest.pval <- stats::wilcox.test(inside, outside)$p.value
        
        res.pval[i,] <- c(clu, round(median(inside), digits = 2),
                          round(median(outside), digits = 2), utest.pval)
        i = i+1
    }
    as.data.frame(res.pval) %>% 
        dplyr::mutate("Adj.Pvalue" = stats::p.adjust(Pvalue, method = "fdr")) %>%
        na.omit
    return(res.pval)
}

#' @export
#' @import ggplot2
#' @importFrom mixOmics color.mixo
plot.proportionality <- function(x, ...){
    ggplot2::ggplot(data = x$propr.distance.w.cluster, 
                    aes(x=as.factor(cluster1), y=value, col=insideout)) + 
        geom_boxplot() + theme_bw() + 
        xlab("Cluster ID") + 
        ylab("Proportionality distance") + 
        labs(color = "Proportionality distance") +
        scale_color_manual(values = mixOmics::color.mixo(1:2))
}

# based on the works of Thom Quinn (https://github.com/tpq/propr)
get_phi_s <- function(counts){
    ct <- counts
    
    if (any(as.matrix(counts) == 0)) {
        #message("Alert: Replacing 0s with next smallest value.")
        zeros <- ct == 0
        ct[zeros] <- min(ct[!zeros])
    }
    #ivar = "clr"
    #use <- propr::ivar2index(ct, ivar)
    use <- 1:ncol(counts)
    
    
    logX <- log(ct)
    logSet <- logX[, use, drop = FALSE]
    ref <- rowMeans(logSet)
    lr <- sweep(logX, 1, ref, "-")

    mat <- as.data.frame(to_lr2phs(lr))
    colnames(mat) <- colnames(lr)
    rownames(mat) <- colnames(lr)
    return(mat)
}


to_lr2phs <- function(lr){
    # Calculate phs = var(a-b)/var(a+b)
    nfeats <- ncol(lr);
    mat <- matrix(nrow = nfeats, ncol = nfeats)
    for(i in 1:nfeats){
        for(j in 1:nfeats){
            if(i == j){
                mat[i,j] <- 0
            } else {
                a <- lr[,i]; b <- lr[,j]
                mat[i,j] <- var(a-b)/var(a+b)
            }
        }
    }
    # return the same matrix as propr
    return(mat)
}