R/select_rho_partial.R

In ressomlab/INDEED: Interactive Visualization of Integrated Differential Expression and Differential Network Analysis for Biomarker Candidate Selection

#' @title Data preprocessing for partial correlation analysis
#' @description A method that integrates differential expression (DE) analysis
#'     and differential network (DN) analysis to select biomarker candidates for
#'     cancer studies. select_rho_partial is the pre-processing step for INDEED
#'     partial differential analysis.
#' @param data This is a p*n dataframe that contains the expression levels for all biomolecules and samples.
#' @param class_label This is a 1*n dataframe that contains the class label with 0 for group 1 and 1 for group 2.
#' @param id This is a p*1 dataframe that contains the ID for each biomolecule.
#' @param error_curve This is a boolean value indicating whether to plot the error curve (TRUE) or not (FALSE). 
#'     The default is TRUE.
#' @examples select_rho_partial(data = Met_GU, class_label = Met_Group_GU, id = Met_name_GU, 
#'     error_curve = TRUE)
#' @return A list of processed data for the next step, and an error curve to select optimal rho value
#'     for graphical lasso.
#' @import devtools
#' @importFrom glasso glasso
#' @importFrom stats qnorm cor quantile var sd glm 
#' @importFrom graphics abline title plot lines 
#' @export

select_rho_partial <- function(data = NULL, class_label = NULL, id = NULL, error_curve = TRUE){
    data_bind <- rbind(data, class_label)
    # Group 1: p*n1
    raw_group_1 <- data_bind[,data_bind[nrow(data_bind),] == 0][1:(nrow(data_bind) - 1),]
    # Group 2: p*n2
    raw_group_2 <- data_bind[,data_bind[nrow(data_bind),] == 1][1:(nrow(data_bind) - 1),]  
    p <- nrow(raw_group_2)
    n_group_1 <- ncol(raw_group_1)
    n_group_2 <- ncol(raw_group_2)
    

    # Z-transform the data for group-specific normalization
    data_group_1 <- scale(t(raw_group_1)) # Group 1: n1*p
    data_group_2 <- scale(t(raw_group_2)) # Group 2: n2*p
    cov_group_1 <- var(data_group_1)
    cov_group_2 <- var(data_group_2)
    
    ## Apply grapihcal LASSO given data_group_1 and data_group_2
    #  Group 1 first
    n_fold <- 5 # number of folds
    rho <- exp(seq(log(1), log(0.01), length.out = 20))

    # group 1
    # draw error curve
    error_group1 <- choose_rho(data_group_1, n_fold, rho)
    if(error_curve == TRUE){
        par(mfrow = c(1, 2))
        plot(rho, error_group1$log.cv, xlab = expression(rho), ylab = "Error")
        lines(rho, error_group1$log.cv)
        title(main = paste("Group 1 error curve", "using cross validation", sep = "\n"))
        # chosse optimal rho
        abline(v = rho[error_group1$log.cv == min(error_group1$log.cv)], col = "red", lty = 3)
        # one standard error rule
        abline(h = min(error_group1$log.cv) + 
                   error_group1$log.rho[error_group1$log.cv == min(error_group1$log.cv)], 
               col = "blue")
    }
    rho[error_group1$log.cv == min(error_group1$log.cv)] # rho based on minimum rule

    # rhos that are under blue line
    rho_under_blue <- rho[error_group1$log.cv < 
                              min(error_group1$log.cv) + 
                              error_group1$log.rho[error_group1$log.cv == 
                                                       min(error_group1$log.cv)]]

    # rhos that are on the right of the red line
    rho_right_red <- rho[rho > rho[error_group1$log.cv == min(error_group1$log.cv)]]

    #rhos that are under blue line and to the right of the red line
    rho_one_std_group1 <- max(intersect(rho_under_blue, rho_right_red ))
    rho_min_rule_group1 <- rho[error_group1$log.cv == min(error_group1$log.cv)]

    # group 2
    # draw error curve
    error_group2 <- choose_rho(data_group_2, n_fold, rho)
    if(error_curve == TRUE){
        plot(rho, error_group2$log.cv, xlab = expression(rho), ylab = "Error")
        lines(rho, error_group2$log.cv)
        title(main = paste("Group 2 error curve", "using cross validation", sep="\n"))
        # choOse optimal rho
        abline(v = rho[error_group2$log.cv == min(error_group2$log.cv)], col = "red", lty = 3)
        # one standard error rule
        abline(h = min(error_group2$log.cv) + 
                   error_group2$log.rho[error_group2$log.cv == 
                                            min(error_group2$log.cv)], col = "blue")
    }
    rho[error_group2$log.cv == min(error_group2$log.cv)] # rho based on minimum rule

    # rhos that are under blue line
    rho_under_blue <- rho[error_group2$log.cv < min(error_group2$log.cv) + 
                              error_group2$log.rho[error_group2$log.cv == 
                                                       min(error_group2$log.cv)]]

    # rhos that are on the right of the red line
    rho_right_red <- rho[rho > rho[error_group2$log.cv == min(error_group2$log.cv)]]

    #rhos that are under blue line and to the right of the red line
    rho_one_std_group2 <- max(intersect(rho_under_blue, rho_right_red )) # export as global
    rho_min_rule_group2 <- rho[error_group2$log.cv == min(error_group2$log.cv)] # export as global

    rho_df <- data.frame(c(rho_one_std_group1, rho_min_rule_group1), c(rho_one_std_group2, 
                                                                       rho_min_rule_group2), 
                         row.names = c("one standard error","minimum"))
    colnames(rho_df)<-c('group1', "group2")
    
    data_list <- list(p = p, rho_table = rho_df, cov_group_1 = cov_group_1, 
                      cov_group_2 = cov_group_2, n_group_1 = n_group_1, n_group_2 = n_group_2, 
                      data_group_1 = data_group_1, data_group_2 = data_group_2, 
                      class_label = class_label, id = id, data = data)
    return(data_list)
}