# center with 'colMeans()' - ref: http://gastonsanchez.com/visually-enforced/how-to/2014/01/15/Center-data-in-R/ center_colmeans <- function(x) { xcenter = colMeans(x) x - rep(xcenter, rep.int(nrow(x), ncol(x))) } #### OPLS-DA algoC <- "nipals" do_detail_plot <- function(x_dataMatrix, x_predictor, x_is_match, x_algorithm, x_prefix, x_show_labels, x_progress = print, x_env) { off <- function(x) if (x_show_labels) x else 0 salience_lookup <- x_env$salience_lookup salient_rcv_lookup <- x_env$salient_rcv_lookup # x_progress("head(salience_df): ", head(salience_df)) # x_progress("head(salience): ", head(salience)) if (x_is_match && ncol(x_dataMatrix) > 0 && length(unique(x_predictor))> 1) { my_oplsda <- opls( x = x_dataMatrix , y = x_predictor , algoC = x_algorithm , predI = 1 , orthoI = if (ncol(x_dataMatrix) > 1) 1 else 0 , printL = FALSE , plotL = FALSE ) my_oplsda_suppLs_y_levels <- levels(as.factor(my_oplsda@suppLs$y)) fctr_lvl_1 <- my_oplsda_suppLs_y_levels[1] fctr_lvl_2 <- my_oplsda_suppLs_y_levels[2] my_cor_vs_cov <- cor_vs_cov( matrix_x = x_dataMatrix , ropls_x = my_oplsda ) with( my_cor_vs_cov , { min_x <- min(covariance) max_x <- max(covariance) lim_x <- max(sapply(X=c(min_x, max_x), FUN=abs)) covariance <- covariance / lim_x lim_x <- 1.2 main_label <- sprintf("%s for levels %s versus %s", x_prefix, fctr_lvl_1, fctr_lvl_2) # print("main_label") # print(main_label) main_cex <- min(1.0, 46.0/nchar(main_label)) # "It is generally accepted that a variable should be selected if vj>1, [27–29], # but a proper threshold between 0.83 and 1.21 can yield more relevant variables according to [28]." # (Mehmood 2012 doi:10.1186/1748-7188-6-27) vipco <- pmax(0, pmin(1,(vip4p-0.83)/(1.21-0.83))) alpha <- 0.1 + 0.4 * vipco red <- as.numeric(correlation < 0) * vipco blue <- as.numeric(correlation > 0) * vipco minus_cor <- -correlation minus_cov <- -covariance # cex <- salience_lookup(feature = names(minus_cor)) # cex <- 0.25 + (1.25 * cex / max(cex)) cex <- 0.75 plot( y = minus_cor , x = minus_cov , type="p" , xlim=c(-lim_x, lim_x + off(0.1)) , ylim=c(-1.0 - off(0.1), 1.0) , xlab = sprintf("relative covariance(feature,t1)") , ylab = sprintf("correlation(feature,t1)") , main = main_label , cex.main = main_cex , cex = cex , pch = 16 , col = rgb(blue = blue, red = red, green = 0, alpha = alpha) ) low_x <- -0.7 * lim_x high_x <- 0.7 * lim_x text(x = low_x, y = -0.15, labels = fctr_lvl_1) text(x = high_x, y = 0.15, labels = fctr_lvl_2) if (x_show_labels) { text( y = minus_cor - 0.013 , x = minus_cov + 0.020 , cex = 0.3 , labels = tsv1$featureID , col = rgb(blue = blue, red = red, green = 0, alpha = 0.2 + 0.8 * alpha) , srt = -30 # slant 30 degrees downward , adj = 0 # left-justified ) } } ) typeVc <- c("correlation", # 1 "outlier", # 2 "overview", # 3 "permutation", # 4 "predict-train", # 5 "predict-test", # 6 "summary", # 7 = c(2,3,4,9) "x-loading", # 8 "x-score", # 9 "x-variance", # 10 "xy-score", # 11 "xy-weight" # 12 ) # [c(3,8,9)] # [c(4,3,8,9)] if ( length(my_oplsda@orthoVipVn) > 0 ) { my_typevc <- typeVc[c(9,3,8)] } else { my_typevc <- c("(dummy)","overview","(dummy)") } for (my_type in my_typevc) { if (my_type %in% typeVc) { # print(sprintf("plotting type %s", my_type)) plot( x = my_oplsda , typeVc = my_type , parCexN = 0.4 , parDevNewL = FALSE , parLayL = TRUE , parEllipsesL = TRUE ) } else { # print("plotting dummy graph") plot(x=1, y=1, xaxt="n", yaxt="n", xlab="", ylab="", type="n") text(x=1, y=1, labels="no orthogonal projection is possible") } } return (my_cor_vs_cov) } else { # x_progress(sprintf("x_is_match = %s, ncol(x_dataMatrix) = %d, length(unique(x_predictor)) = %d",x_is_match, ncol(x_dataMatrix), length(unique(x_predictor)))) return (NULL) } } # S-PLOT and OPLS reference: Wiklund_2008 doi:10.1021/ac0713510 corcov_calc <- function(calc_env, failure_action = stop, progress_action = function(x){}, corcov_tsv_action = function(t){}, salience_tsv_action = function(t){}) { if ( ! is.environment(calc_env) ) { failure_action("corcov_calc: fatal error - 'calc_env' is not an environment") return ( FALSE ) } if ( ! is.function(corcov_tsv_action) ) { failure_action("corcov_calc: fatal error - 'corcov_tsv_action' is not a function") return ( FALSE ) } if ( ! is.function(salience_tsv_action) ) { failure_action("salience_calc: fatal error - 'salience_tsv_action' is not a function") return ( FALSE ) } # extract parameters from the environment vrbl_metadata <- calc_env$vrbl_metadata vrbl_metadata_names <- vrbl_metadata[,1] smpl_metadata <- calc_env$smpl_metadata data_matrix <- calc_env$data_matrix pairSigFeatOnly <- calc_env$pairSigFeatOnly facC <- calc_env$facC tesC <- calc_env$tesC # extract the levels from the environment originalLevCSV <- levCSV <- calc_env$levCSV # matchingC is one of { "none", "wildcard", "regex" } matchingC <- calc_env$matchingC labelFeatures <- calc_env$labelFeatures # arg/env checking if (!(facC %in% names(smpl_metadata))) { failure_action(sprintf("bad parameter! Factor name '%s' not found in sampleMetadata", facC)) return ( FALSE ) } # calculate salience_df as data.frame(feature, max_level, max_median, max_rcv, mean_median, salience, salient_rcv) salience_df <- calc_env$salience_df <- w4msalience( data_matrix = data_matrix , sample_class = smpl_metadata[,facC] , failure_action = failure_action ) salience_tsv_action({ my_df <- data.frame( featureID = salience_df$feature , salientLevel = salience_df$max_level , salientRCV = salience_df$salient_rcv , salience = salience_df$salience ) my_df[order(-my_df$salience),] }) salience <- salience_df$salience names(salience) <- salience_df$feature salience_lookup <- calc_env$salience_lookup <- function(feature) unname(salience[feature]) salient_rcv <- salience_df$salient_rcv names(salient_rcv) <- salience_df$feature salient_rcv_lookup <- calc_env$salient_rcv_lookup <- function(feature) unname(salient_rcv[feature]) salient_level <- salience_df$max_level names(salient_level) <- salience_df$feature salient_level_lookup <- calc_env$salient_level_lookup <- function(feature) unname(salient_level[feature]) # transform wildcards to regexen if (matchingC == "wildcard") { # strsplit(x = "hello,wild,world", split = ",") levCSV <- gsub("[.]", "[.]", levCSV) levCSV <- strsplit(x = levCSV, split = ",") levCSV <- sapply(levCSV, utils::glob2rx, trim.tail = FALSE) levCSV <- paste(levCSV, collapse = ",") } # function to determine whether level is a member of levCSV isLevelSelected <- function(lvl) { matchFun <- if (matchingC != "none") grepl else `==` return( Reduce( f = "||" , x = sapply( X = strsplit( x = levCSV , split = "," , fixed = TRUE )[[1]] , FUN = matchFun , lvl ) ) ) } # transpose matrix because ropls matrix is the transpose of XCMS matrix # Wiklund_2008 centers and pareto-scales data before OPLS-DA S-plot # center cdm <- center_colmeans(t(data_matrix)) # pareto-scale my_scale <- sqrt(apply(cdm, 2, sd, na.rm=TRUE)) scdm <- sweep(cdm, 2, my_scale, "/") # pattern to match column names like k10_kruskal_k4.k3_sig col_pattern <- sprintf('^%s_%s_(.*)[.](.*)_sig$', facC, tesC) # column name like k10_kruskal_sig intersample_sig_col <- sprintf('%s_%s_sig', facC, tesC) # get the facC column from sampleMetadata, dropping to one dimension smpl_metadata_facC <- smpl_metadata[,facC] # allocate a slot in the environment for the contrast_list, each element of which will be a data.frame with this structure: # - feature ID # - value1 # - value2 # - Wiklund_2008 correlation # - Wiklund_2008 covariance # - Wiklund_2008 VIP calc_env$contrast_list <- list() did_plot <- FALSE if (tesC != "none") { # for each column name, extract the parts of the name matched by 'col_pattern', if any the_colnames <- colnames(vrbl_metadata) if ( ! Reduce( f = "||", x = grepl(tesC, the_colnames) ) ) { failure_action(sprintf("bad parameter! variableMetadata must contain results of W4M Univariate test '%s'.", tesC)) return ( FALSE ) } col_matches <- lapply( X = the_colnames, FUN = function(x) { regmatches( x, regexec(col_pattern, x) )[[1]] } ) ## first contrast each selected level with all other levels combined into one "super-level" ## # process columns matching the pattern level_union <- c() for (i in 1:length(col_matches)) { col_match <- col_matches[[i]] if (length(col_match) > 0) { # it's an actual match; extract the pieces, e.g., k10_kruskal_k4.k3_sig vrbl_metadata_col <- col_match[1] # ^^^^^^^^^^^^^^^^^^^^^ # Column name fctr_lvl_1 <- col_match[2] # ^^ # Factor-level 1 fctr_lvl_2 <- col_match[3] # ^^ # Factor-level 2 # only process this column if both factors are members of lvlCSV is_match <- isLevelSelected(fctr_lvl_1) && isLevelSelected(fctr_lvl_2) if (is_match) { level_union <- c(level_union, col_match[2], col_match[3]) } } } level_union <- unique(sort(level_union)) overall_significant <- 1 == vrbl_metadata[,intersample_sig_col] if ( length(level_union) > 1 ) { chosen_samples <- smpl_metadata_facC %in% level_union chosen_facC <- as.character(smpl_metadata_facC[chosen_samples]) col_selector <- vrbl_metadata_names[ overall_significant ] my_matrix <- scdm[ chosen_samples, col_selector, drop = FALSE ] fctr_lvl_2 <- "other" for ( fctr_lvl_1 in level_union[1:length(level_union)] ) { progress_action(sprintf("calculating/plotting contrast of %s vs. %s", fctr_lvl_1, fctr_lvl_2)) predictor <- sapply( X = chosen_facC, FUN = function(fac) if ( fac == fctr_lvl_1 ) fctr_lvl_1 else fctr_lvl_2 ) my_cor_cov <- do_detail_plot( x_dataMatrix = my_matrix , x_predictor = predictor , x_is_match = is_match , x_algorithm = algoC , x_prefix = if (pairSigFeatOnly) "Significantly contrasting features" else "Significant features" , x_show_labels = labelFeatures , x_progress = progress_action , x_env = calc_env ) if ( is.null(my_cor_cov) ) { progress_action("NOTHING TO PLOT.") } else { tsv <- my_cor_cov$tsv1 # tsv$salientLevel <- salient_level_lookup(tsv$featureID) # tsv$salientRCV <- salient_rcv_lookup(tsv$featureID) # tsv$salience <- salience_lookup(tsv$featureID) tsv["level1Level2Sig"] <- vrbl_metadata[ match(tsv$featureID, vrbl_metadata_names), vrbl_metadata_col ] corcov_tsv_action(tsv) did_plot <- TRUE } } } ## next, contrast each selected level with each of the other levels individually ## # process columns matching the pattern for (i in 1:length(col_matches)) { # for each potential match of the pattern col_match <- col_matches[[i]] if (length(col_match) > 0) { # it's an actual match; extract the pieces, e.g., k10_kruskal_k4.k3_sig vrbl_metadata_col <- col_match[1] # ^^^^^^^^^^^^^^^^^^^^^ # Column name fctr_lvl_1 <- col_match[2] # ^^ # Factor-level 1 fctr_lvl_2 <- col_match[3] # ^^ # Factor-level 2 # only process this column if both factors are members of lvlCSV is_match <- isLevelSelected(fctr_lvl_1) && isLevelSelected(fctr_lvl_2) progress_action(sprintf("calculating/plotting contrast of %s vs. %s", fctr_lvl_1, fctr_lvl_2)) # TODO delete next line displaying currently-processed column # cat(sprintf("%s %s %s %s\n", vrbl_metadata_col, fctr_lvl_1, fctr_lvl_2, is_match)) # choose only samples with one of the two factors for this column chosen_samples <- smpl_metadata_facC %in% c(fctr_lvl_1, fctr_lvl_2) predictor <- smpl_metadata_facC[chosen_samples] # extract only the significantly-varying features and the chosen samples fully_significant <- 1 == vrbl_metadata[,vrbl_metadata_col] * vrbl_metadata[,intersample_sig_col] col_selector <- vrbl_metadata_names[ if ( pairSigFeatOnly ) fully_significant else overall_significant ] my_matrix <- scdm[ chosen_samples, col_selector, drop = FALSE ] my_cor_cov <- do_detail_plot( x_dataMatrix = my_matrix , x_predictor = predictor , x_is_match = is_match , x_algorithm = algoC , x_prefix = if (pairSigFeatOnly) "Significantly contrasting features" else "Significant features" , x_show_labels = labelFeatures , x_progress = progress_action , x_env = calc_env ) if ( is.null(my_cor_cov) ) { progress_action("NOTHING TO PLOT.") } else { tsv <- my_cor_cov$tsv1 # tsv$salientLevel <- salient_level_lookup(tsv$featureID) # tsv$salientRCV <- salient_rcv_lookup(tsv$featureID) # tsv$salience <- salience_lookup(tsv$featureID) tsv["level1Level2Sig"] <- vrbl_metadata[ match(tsv$featureID, vrbl_metadata_names), vrbl_metadata_col ] corcov_tsv_action(tsv) did_plot <- TRUE } } } } else { # tesC == "none" level_union <- unique(sort(smpl_metadata_facC)) if ( length(level_union) > 1 ) { ## pass 1 - contrast each selected level with all other levels combined into one "super-level" ## completed <- c() lapply( X = level_union , FUN = function(x) { fctr_lvl_1 <- x[1] fctr_lvl_2 <- { if ( fctr_lvl_1 %in% completed ) return("DUMMY") # strF(completed) completed <<- c(completed, fctr_lvl_1) setdiff(level_union, fctr_lvl_1) } chosen_samples <- smpl_metadata_facC %in% c(fctr_lvl_1, fctr_lvl_2) fctr_lvl_2 <- "other" # print( sprintf("sum(chosen_samples) %d, factor_level_2 %s", sum(chosen_samples), fctr_lvl_2) ) progress_action(sprintf("calculating/plotting contrast of %s vs. %s", fctr_lvl_1, fctr_lvl_2)) if (length(unique(chosen_samples)) < 1) { progress_action("NOTHING TO PLOT...") } else { chosen_facC <- as.character(smpl_metadata_facC[chosen_samples]) predictor <- sapply( X = chosen_facC, FUN = function(fac) if ( fac == fctr_lvl_1 ) fctr_lvl_1 else fctr_lvl_2 ) my_matrix <- scdm[ chosen_samples, , drop = FALSE ] # only process this column if both factors are members of lvlCSV is_match <- isLevelSelected(fctr_lvl_1) my_cor_cov <- do_detail_plot( x_dataMatrix = my_matrix , x_predictor = predictor , x_is_match = is_match , x_algorithm = algoC , x_prefix = "Features" , x_show_labels = labelFeatures , x_progress = progress_action , x_env = calc_env ) if ( is.null(my_cor_cov) ) { progress_action("NOTHING TO PLOT") } else { tsv <- my_cor_cov$tsv1 # tsv$salientLevel <- salient_level_lookup(tsv$featureID) # tsv$salientRCV <- salient_rcv_lookup(tsv$featureID) # tsv$salience <- salience_lookup(tsv$featureID) corcov_tsv_action(tsv) did_plot <<- TRUE } } #print("baz") "dummy" # need to return a value; otherwise combn fails with an error } ) ## pass 2 - contrast each selected level with each of the other levels individually ## completed <- c() utils::combn( x = level_union , m = 2 , FUN = function(x) { fctr_lvl_1 <- x[1] fctr_lvl_2 <- x[2] chosen_samples <- smpl_metadata_facC %in% c(fctr_lvl_1, fctr_lvl_2) # print( sprintf("sum(chosen_samples) %d, factor_level_2 %s", sum(chosen_samples), fctr_lvl_2) ) progress_action(sprintf("calculating/plotting contrast of %s vs. %s", fctr_lvl_1, fctr_lvl_2)) if (length(unique(chosen_samples)) < 1) { progress_action("NOTHING TO PLOT...") } else { chosen_facC <- as.character(smpl_metadata_facC[chosen_samples]) predictor <- chosen_facC my_matrix <- scdm[ chosen_samples, , drop = FALSE ] # only process this column if both factors are members of lvlCSV is_match <- isLevelSelected(fctr_lvl_1) && isLevelSelected(fctr_lvl_2) my_cor_cov <- do_detail_plot( x_dataMatrix = my_matrix , x_predictor = predictor , x_is_match = is_match , x_algorithm = algoC , x_prefix = "Features" , x_show_labels = labelFeatures , x_progress = progress_action , x_env = calc_env ) if ( is.null(my_cor_cov) ) { progress_action("NOTHING TO PLOT") } else { tsv <- my_cor_cov$tsv1 # tsv$salientLevel <- salient_level_lookup(tsv$featureID) # tsv$salientRCV <- salient_rcv_lookup(tsv$featureID) # tsv$salience <- salience_lookup(tsv$featureID) corcov_tsv_action(tsv) did_plot <<- TRUE } } #print("baz") "dummy" # need to return a value; otherwise combn fails with an error } ) } else { progress_action("NOTHING TO PLOT....") } } if (!did_plot) { failure_action(sprintf("bad parameter! sampleMetadata must have at least two levels of factor '%s' matching '%s'", facC, originalLevCSV)) return ( FALSE ) } return ( TRUE ) } # Calculate data for correlation-versus-covariance plot # Adapted from: # Wiklund_2008 doi:10.1021/ac0713510 # Galindo_Prieto_2014 doi:10.1002/cem.2627 # https://github.com/HegemanLab/extra_tools/blob/master/generic_PCA.R cor_vs_cov <- function(matrix_x, ropls_x) { x_class <- class(ropls_x) if ( !( as.character(x_class) == "opls" ) ) { # || !( attr(class(x_class),"package") == "ropls" ) ) stop( "cor_vs_cov: Expected ropls_x to be of class ropls::opls but instead it was of class ", as.character(x_class) ) } result <- list() # suppLs$algoC - Character: algorithm used - "svd" for singular value decomposition; "nipals" for NIPALS if ( ropls_x@suppLs$algoC == "nipals") { # Equations (1) and (2) from *Supplement to* Wiklund 2008, doi:10.1021/ac0713510 mag <- function(one_dimensional) sqrt(sum(one_dimensional * one_dimensional)) mag_xi <- sapply(X = 1:ncol(matrix_x), FUN = function(x) mag(matrix_x[,x])) score_matrix <- ropls_x@scoreMN score_matrix_transposed <- t(score_matrix) score_matrix_magnitude <- mag(score_matrix) result$covariance <- score_matrix_transposed %*% matrix_x / ( score_matrix_magnitude * score_matrix_magnitude ) result$correlation <- score_matrix_transposed %*% matrix_x / ( score_matrix_magnitude * mag_xi ) } else { # WARNING - untested code - I don't have test data to exercise this branch # Equations (1) and (2) from Wiklund 2008, doi:10.1021/ac0713510 # scoreMN - Numerical matrix of x scores (T; dimensions: nrow(x) x predI) X = TP' + E; Y = TC' + F score_matrix <- ropls_x@scoreMN score_matrix_transposed <- t(score_matrix) cov_divisor <- nrow(matrix_x) - 1 result$covariance <- sapply( X = 1:ncol(matrix_x) , FUN = function(x) score_matrix_transposed %*% matrix_x[,x] / cov_divisor ) score_sd <- sapply( X = 1:ncol(score_matrix) , FUN = function(x) sd(score_matrix[,x]) ) # xSdVn - Numerical vector: variable standard deviations of the 'x' matrix xSdVn <- ropls_x@xSdVn result$correlation <- sapply( X = 1:ncol(matrix_x) , FUN = function(x) { ( score_matrix_transposed / score_sd ) %*% ( matrix_x[,x] / (xSdVn[x] * cov_divisor) ) } ) } result$correlation <- result$correlation[1,,drop = TRUE] result$covariance <- result$covariance[1,,drop = TRUE] # Variant 4 of Variable Influence on Projection for OPLS from Galindo_Prieto_2014 # Length = number of features; labels = feature identifiers. (The same is true for $correlation and $covariance.) result$vip4p <- as.numeric(ropls_x@vipVn) result$vip4o <- as.numeric(ropls_x@orthoVipVn) # get the level names level_names <- sort(levels(as.factor(ropls_x@suppLs$y))) feature_count <- length(ropls_x@vipVn) result$level1 <- rep.int(x = level_names[1], times = feature_count) result$level2 <- rep.int(x = level_names[2], times = feature_count) # print(sprintf("sd(covariance) = %f; sd(correlation) = %f", sd(result$covariance), sd(result$correlation))) superresult <- list() if (length(result$vip4o) == 0) result$vip4o <- NA superresult$tsv1 <- data.frame( featureID = names(ropls_x@vipVn) , factorLevel1 = result$level1 , factorLevel2 = result$level2 , correlation = result$correlation , covariance = result$covariance , vip4p = result$vip4p , vip4o = result$vip4o , row.names = NULL ) rownames(superresult$tsv1) <- superresult$tsv1$featureID superresult$covariance <- result$covariance superresult$correlation <- result$correlation superresult$vip4p <- result$vip4p superresult$vip4o <- result$vip4o superresult$details <- result # #print(superresult$tsv1) result$superresult <- superresult # Include thise in case future consumers of this routine want to use it in currently unanticipated ways result$oplsda <- ropls_x result$predictor <- ropls_x@suppLs$y # in case future consumers of this routine want to use it in currently unanticipated ways return (superresult) }