Mercurial > repos > eschen42 > w4mcorcov

comparison w4mcorcov_calc.R @ 2:a06344808ffc draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/HegemanLab/w4mcorcov_galaxy_wrapper/tree/master commit ce5178ce51b80f242d24db555044e6afc530ac99
author eschen42
date Sat, 11 Nov 2017 00:08:20 -0500
parents e25fd8a13665
children 61935618f92c
comparison
equal deleted inserted replaced
1:e25fd8a13665 2:a06344808ffc
5 } 5 }
6 6
7 #### OPLS-DA 7 #### OPLS-DA
8 algoC <- "nipals" 8 algoC <- "nipals"
9 9
10 do_detail_plot <- function(x_dataMatrix, x_predictor, x_is_match, x_algorithm, x_prefix, x_show_labels, x_progress = print, x_env) { 10 do_detail_plot <- function(x_dataMatrix, x_predictor, x_is_match, x_algorithm, x_prefix, x_show_labels, x_show_loado_labels, x_progress = print, x_env) {
11 off <- function(x) if (x_show_labels) x else 0 11 off <- function(x) if (x_show_labels == "0") 0 else x
12 salience_lookup <- x_env$salience_lookup
13 salient_rcv_lookup <- x_env$salient_rcv_lookup
14 # x_progress("head(salience_df): ", head(salience_df))
15 # x_progress("head(salience): ", head(salience))
16 if (x_is_match && ncol(x_dataMatrix) > 0 && length(unique(x_predictor))> 1) { 12 if (x_is_match && ncol(x_dataMatrix) > 0 && length(unique(x_predictor))> 1) {
17 my_oplsda <- opls( 13 my_oplsda <- opls(
18 x = x_dataMatrix 14 x = x_dataMatrix
19 , y = x_predictor 15 , y = x_predictor
20 , algoC = x_algorithm 16 , algoC = x_algorithm
37 max_x <- max(covariance) 33 max_x <- max(covariance)
38 lim_x <- max(sapply(X=c(min_x, max_x), FUN=abs)) 34 lim_x <- max(sapply(X=c(min_x, max_x), FUN=abs))
39 covariance <- covariance / lim_x 35 covariance <- covariance / lim_x
40 lim_x <- 1.2 36 lim_x <- 1.2
41 main_label <- sprintf("%s for levels %s versus %s", x_prefix, fctr_lvl_1, fctr_lvl_2) 37 main_label <- sprintf("%s for levels %s versus %s", x_prefix, fctr_lvl_1, fctr_lvl_2)
42 # print("main_label")
43 # print(main_label)
44 main_cex <- min(1.0, 46.0/nchar(main_label)) 38 main_cex <- min(1.0, 46.0/nchar(main_label))
45 # "It is generally accepted that a variable should be selected if vj>1, [27–29], 39 # "It is generally accepted that a variable should be selected if vj>1, [27–29],
46 # but a proper threshold between 0.83 and 1.21 can yield more relevant variables according to [28]." 40 # but a proper threshold between 0.83 and 1.21 can yield more relevant variables according to [28]."
47 # (Mehmood 2012 doi:10.1186/1748-7188-6-27) 41 # (Mehmood 2012 doi:10.1186/1748-7188-6-27)
48 vipco <- pmax(0, pmin(1,(vip4p-0.83)/(1.21-0.83))) 42 vipco <- pmax(0, pmin(1,(vip4p-0.83)/(1.21-0.83)))
49 alpha <- 0.1 + 0.4 * vipco 43 alpha <- 0.1 + 0.4 * vipco
50 red <- as.numeric(correlation < 0) * vipco 44 red <- as.numeric(correlation > 0) * vipco
51 blue <- as.numeric(correlation > 0) * vipco 45 blue <- as.numeric(correlation < 0) * vipco
52 minus_cor <- -correlation 46 plus_cor <- correlation
53 minus_cov <- -covariance 47 plus_cov <- covariance
54 # cex <- salience_lookup(feature = names(minus_cor))
55 # cex <- 0.25 + (1.25 * cex / max(cex))
56 cex <- 0.75 48 cex <- 0.75
57 plot( 49 plot(
58 y = minus_cor 50 y = plus_cor
59 , x = minus_cov 51 , x = plus_cov
60 , type="p" 52 , type="p"
61 , xlim=c(-lim_x, lim_x + off(0.1)) 53 , xlim=c(-lim_x, lim_x + off(0.2))
62 , ylim=c(-1.0 - off(0.1), 1.0) 54 , ylim=c(-1.0 - off(0.2), 1.0)
63 , xlab = sprintf("relative covariance(feature,t1)") 55 , xlab = sprintf("relative covariance(feature,t1)")
64 , ylab = sprintf("correlation(feature,t1)") 56 , ylab = sprintf("correlation(feature,t1)")
65 , main = main_label 57 , main = main_label
66 , cex.main = main_cex 58 , cex.main = main_cex
67 , cex = cex 59 , cex = cex
68 , pch = 16 60 , pch = 16
69 , col = rgb(blue = blue, red = red, green = 0, alpha = alpha) 61 , col = rgb(blue = blue, red = red, green = 0, alpha = alpha)
70 ) 62 )
71 low_x <- -0.7 * lim_x 63 low_x <- -0.7 * lim_x
72 high_x <- 0.7 * lim_x 64 high_x <- 0.7 * lim_x
73 text(x = low_x, y = -0.15, labels = fctr_lvl_1) 65 text(x = low_x, y = -0.05, labels = fctr_lvl_1, col = "blue")
74 text(x = high_x, y = 0.15, labels = fctr_lvl_2) 66 text(x = high_x, y = 0.05, labels = fctr_lvl_2, col = "red")
75 if (x_show_labels) { 67 if ( x_show_labels != "0" ) {
68 my_loadp <- loadp
69 my_loado <- loado
70 names(my_loadp) <- tsv1$featureID
71 names(my_loado) <- tsv1$featureID
72 if ( x_show_labels == "ALL" ) {
73 n_labels <- length(loadp)
74 } else {
75 n_labels <- as.numeric(x_show_labels)
76 }
77 n_labels <- min( n_labels, (1 + length(loadp)) / 2 )
78 labels_to_show <- c(
79 names(head(sort(my_loadp),n = n_labels))
80 , names(tail(sort(my_loadp),n = n_labels))
81 )
82 if ( x_show_loado_labels ) {
83 labels_to_show <- c(
84 labels_to_show
85 , names(head(sort(my_loado),n = n_labels))
86 , names(tail(sort(my_loado),n = n_labels))
87 )
88 }
89 labels <- unname(sapply( X = tsv1$featureID, FUN = function(x) if( x %in% labels_to_show ) x else "" ))
76 text( 90 text(
77 y = minus_cor - 0.013 91 y = plus_cor - 0.013
78 , x = minus_cov + 0.020 92 , x = plus_cov + 0.020
79 , cex = 0.3 93 , cex = 0.3
80 , labels = tsv1$featureID 94 , labels = labels
81 , col = rgb(blue = blue, red = red, green = 0, alpha = 0.2 + 0.8 * alpha) 95 , col = rgb(blue = blue, red = red, green = 0, alpha = 0.2 + 0.8 * alpha)
82 , srt = -30 # slant 30 degrees downward 96 , srt = -30 # slant 30 degrees downward
83 , adj = 0 # left-justified 97 , adj = 0 # left-justified
84 ) 98 )
85 } 99 }
153 # extract the levels from the environment 167 # extract the levels from the environment
154 originalLevCSV <- levCSV <- calc_env$levCSV 168 originalLevCSV <- levCSV <- calc_env$levCSV
155 # matchingC is one of { "none", "wildcard", "regex" } 169 # matchingC is one of { "none", "wildcard", "regex" }
156 matchingC <- calc_env$matchingC 170 matchingC <- calc_env$matchingC
157 labelFeatures <- calc_env$labelFeatures 171 labelFeatures <- calc_env$labelFeatures
172 labelOrthoFeatures <- calc_env$labelOrthoFeatures
158 173
159 # arg/env checking 174 # arg/env checking
160 if (!(facC %in% names(smpl_metadata))) { 175 if (!(facC %in% names(smpl_metadata))) {
161 failure_action(sprintf("bad parameter! Factor name '%s' not found in sampleMetadata", facC)) 176 failure_action(sprintf("bad parameter! Factor name '%s' not found in sampleMetadata", facC))
162 return ( FALSE ) 177 return ( FALSE )
163 } 178 }
164 179
180 mz <- vrbl_metadata$mz
181 names(mz) <- vrbl_metadata$variableMetadata
182 mz_lookup <- function(feature) unname(mz[feature])
183
184 rt <- vrbl_metadata$rt
185 names(rt) <- vrbl_metadata$variableMetadata
186 rt_lookup <- function(feature) unname(rt[feature])
187
165 # calculate salience_df as data.frame(feature, max_level, max_median, max_rcv, mean_median, salience, salient_rcv) 188 # calculate salience_df as data.frame(feature, max_level, max_median, max_rcv, mean_median, salience, salient_rcv)
166 salience_df <- calc_env$salience_df <- w4msalience( 189 salience_df <- calc_env$salience_df <- w4msalience(
167 data_matrix = data_matrix 190 data_matrix = data_matrix
168 , sample_class = smpl_metadata[,facC] 191 , sample_class = smpl_metadata[,facC]
169 , failure_action = failure_action 192 , failure_action = failure_action
172 my_df <- data.frame( 195 my_df <- data.frame(
173 featureID = salience_df$feature 196 featureID = salience_df$feature
174 , salientLevel = salience_df$max_level 197 , salientLevel = salience_df$max_level
175 , salientRCV = salience_df$salient_rcv 198 , salientRCV = salience_df$salient_rcv
176 , salience = salience_df$salience 199 , salience = salience_df$salience
200 , mz = mz_lookup(salience_df$feature)
201 , rt = rt_lookup(salience_df$feature)
177 ) 202 )
178 my_df[order(-my_df$salience),] 203 my_df[order(-my_df$salience),]
179 }) 204 })
180 salience <- salience_df$salience 205
181 names(salience) <- salience_df$feature
182 salience_lookup <- calc_env$salience_lookup <- function(feature) unname(salience[feature])
183 salient_rcv <- salience_df$salient_rcv
184 names(salient_rcv) <- salience_df$feature
185 salient_rcv_lookup <- calc_env$salient_rcv_lookup <- function(feature) unname(salient_rcv[feature])
186 salient_level <- salience_df$max_level
187 names(salient_level) <- salience_df$feature
188 salient_level_lookup <- calc_env$salient_level_lookup <- function(feature) unname(salient_level[feature])
189
190 # transform wildcards to regexen 206 # transform wildcards to regexen
191 if (matchingC == "wildcard") { 207 if (matchingC == "wildcard") {
192 # strsplit(x = "hello,wild,world", split = ",") 208 # strsplit(x = "hello,wild,world", split = ",")
193 levCSV <- gsub("[.]", "[.]", levCSV) 209 levCSV <- gsub("[.]", "[.]", levCSV)
194 levCSV <- strsplit(x = levCSV, split = ",") 210 levCSV <- strsplit(x = levCSV, split = ",")
269 level_union <- c(level_union, col_match[2], col_match[3]) 285 level_union <- c(level_union, col_match[2], col_match[3])
270 } 286 }
271 } 287 }
272 } 288 }
273 level_union <- unique(sort(level_union)) 289 level_union <- unique(sort(level_union))
274 overall_significant <- 1 == vrbl_metadata[,intersample_sig_col] 290 overall_significant <- 1 == ( if (intersample_sig_col %in% colnames(vrbl_metadata)) vrbl_metadata[,intersample_sig_col] else TRUE )
275 if ( length(level_union) > 1 ) { 291 if ( length(level_union) > 2 ) {
276 chosen_samples <- smpl_metadata_facC %in% level_union 292 chosen_samples <- smpl_metadata_facC %in% level_union
277 chosen_facC <- as.character(smpl_metadata_facC[chosen_samples]) 293 chosen_facC <- as.character(smpl_metadata_facC[chosen_samples])
278 col_selector <- vrbl_metadata_names[ overall_significant ] 294 col_selector <- vrbl_metadata_names[ overall_significant ]
279 my_matrix <- scdm[ chosen_samples, col_selector, drop = FALSE ] 295 my_matrix <- scdm[ chosen_samples, col_selector, drop = FALSE ]
280 fctr_lvl_2 <- "other" 296 plot_action <- function(fctr_lvl_1, fctr_lvl_2) {
281 for ( fctr_lvl_1 in level_union[1:length(level_union)] ) {
282 progress_action(sprintf("calculating/plotting contrast of %s vs. %s", fctr_lvl_1, fctr_lvl_2)) 297 progress_action(sprintf("calculating/plotting contrast of %s vs. %s", fctr_lvl_1, fctr_lvl_2))
283 predictor <- sapply( X = chosen_facC, FUN = function(fac) if ( fac == fctr_lvl_1 ) fctr_lvl_1 else fctr_lvl_2 ) 298 predictor <- sapply( X = chosen_facC, FUN = function(fac) if ( fac == fctr_lvl_1 ) fctr_lvl_1 else fctr_lvl_2 )
284 my_cor_cov <- do_detail_plot( 299 my_cor_cov <- do_detail_plot(
285 x_dataMatrix = my_matrix 300 x_dataMatrix = my_matrix
286 , x_predictor = predictor 301 , x_predictor = predictor
287 , x_is_match = is_match 302 , x_is_match = is_match
288 , x_algorithm = algoC 303 , x_algorithm = algoC
289 , x_prefix = if (pairSigFeatOnly) "Significantly contrasting features" else "Significant features" 304 , x_prefix = if (pairSigFeatOnly) "Significantly contrasting features" else "Significant features"
290 , x_show_labels = labelFeatures 305 , x_show_labels = labelFeatures
306 , x_show_loado_labels = labelOrthoFeatures
291 , x_progress = progress_action 307 , x_progress = progress_action
292 , x_env = calc_env 308 , x_env = calc_env
293 ) 309 )
294 if ( is.null(my_cor_cov) ) { 310 if ( is.null(my_cor_cov) ) {
295 progress_action("NOTHING TO PLOT.") 311 progress_action("NOTHING TO PLOT.")
296 } else { 312 } else {
297 tsv <- my_cor_cov$tsv1 313 my_tsv <- my_cor_cov$tsv1
298 # tsv$salientLevel <- salient_level_lookup(tsv$featureID) 314 my_tsv$mz <- mz_lookup(my_tsv$featureID)
299 # tsv$salientRCV <- salient_rcv_lookup(tsv$featureID) 315 my_tsv$rt <- rt_lookup(my_tsv$featureID)
300 # tsv$salience <- salience_lookup(tsv$featureID) 316 my_tsv["level1Level2Sig"] <- vrbl_metadata[ match(my_tsv$featureID, vrbl_metadata_names), vrbl_metadata_col ]
301 tsv["level1Level2Sig"] <- vrbl_metadata[ match(tsv$featureID, vrbl_metadata_names), vrbl_metadata_col ] 317 tsv <<- my_tsv
302 corcov_tsv_action(tsv) 318 corcov_tsv_action(tsv)
303 did_plot <- TRUE 319 did_plot <<- TRUE
304 } 320 }
305 } 321 }
306 } 322 if ( length(level_union) != 2 ) {
307 323 fctr_lvl_2 <- "other"
308 ## next, contrast each selected level with each of the other levels individually ## 324 for ( fctr_lvl_1 in level_union[1:length(level_union)] ) {
309 # process columns matching the pattern 325 plot_action(fctr_lvl_1, fctr_lvl_2)
310 for (i in 1:length(col_matches)) { 326 }
311 # for each potential match of the pattern 327 } else {
312 col_match <- col_matches[[i]] 328 plot_action(fctr_lvl_1 = level_union[1], fctr_lvl_2 = level_union[2])
313 if (length(col_match) > 0) { 329 }
314 # it's an actual match; extract the pieces, e.g., k10_kruskal_k4.k3_sig 330 }
315 vrbl_metadata_col <- col_match[1] # ^^^^^^^^^^^^^^^^^^^^^ # Column name 331
316 fctr_lvl_1 <- col_match[2] # ^^ # Factor-level 1 332 if ( length(level_union) > 1 ) {
317 fctr_lvl_2 <- col_match[3] # ^^ # Factor-level 2 333 ## next, contrast each selected level with each of the other levels individually ##
318 # only process this column if both factors are members of lvlCSV 334 # process columns matching the pattern
319 is_match <- isLevelSelected(fctr_lvl_1) && isLevelSelected(fctr_lvl_2) 335 for (i in 1:length(col_matches)) {
320 progress_action(sprintf("calculating/plotting contrast of %s vs. %s", fctr_lvl_1, fctr_lvl_2)) 336 # for each potential match of the pattern
321 # TODO delete next line displaying currently-processed column 337 col_match <- col_matches[[i]]
322 # cat(sprintf("%s %s %s %s\n", vrbl_metadata_col, fctr_lvl_1, fctr_lvl_2, is_match)) 338 if (length(col_match) > 0) {
323 # choose only samples with one of the two factors for this column 339 # it's an actual match; extract the pieces, e.g., k10_kruskal_k4.k3_sig
324 chosen_samples <- smpl_metadata_facC %in% c(fctr_lvl_1, fctr_lvl_2) 340 vrbl_metadata_col <- col_match[1] # ^^^^^^^^^^^^^^^^^^^^^ # Column name
325 predictor <- smpl_metadata_facC[chosen_samples] 341 fctr_lvl_1 <- col_match[2] # ^^ # Factor-level 1
326 # extract only the significantly-varying features and the chosen samples 342 fctr_lvl_2 <- col_match[3] # ^^ # Factor-level 2
327 fully_significant <- 1 == vrbl_metadata[,vrbl_metadata_col] * vrbl_metadata[,intersample_sig_col] 343 # only process this column if both factors are members of lvlCSV
328 col_selector <- vrbl_metadata_names[ if ( pairSigFeatOnly ) fully_significant else overall_significant ] 344 is_match <- isLevelSelected(fctr_lvl_1) && isLevelSelected(fctr_lvl_2)
329 my_matrix <- scdm[ chosen_samples, col_selector, drop = FALSE ] 345 progress_action(sprintf("calculating/plotting contrast of %s vs. %s", fctr_lvl_1, fctr_lvl_2))
330 my_cor_cov <- do_detail_plot( 346 # TODO delete next line displaying currently-processed column
331 x_dataMatrix = my_matrix 347 # cat(sprintf("%s %s %s %s\n", vrbl_metadata_col, fctr_lvl_1, fctr_lvl_2, is_match))
332 , x_predictor = predictor 348 # choose only samples with one of the two factors for this column
333 , x_is_match = is_match 349 chosen_samples <- smpl_metadata_facC %in% c(fctr_lvl_1, fctr_lvl_2)
334 , x_algorithm = algoC 350 predictor <- smpl_metadata_facC[chosen_samples]
335 , x_prefix = if (pairSigFeatOnly) "Significantly contrasting features" else "Significant features" 351 # extract only the significantly-varying features and the chosen samples
336 , x_show_labels = labelFeatures 352 fully_significant <- 1 == vrbl_metadata[,vrbl_metadata_col] * ( if (intersample_sig_col %in% colnames(vrbl_metadata)) vrbl_metadata[,intersample_sig_col] else TRUE )
337 , x_progress = progress_action 353 col_selector <- vrbl_metadata_names[ if ( pairSigFeatOnly ) fully_significant else overall_significant ]
338 , x_env = calc_env 354 my_matrix <- scdm[ chosen_samples, col_selector, drop = FALSE ]
339 ) 355 my_cor_cov <- do_detail_plot(
340 if ( is.null(my_cor_cov) ) { 356 x_dataMatrix = my_matrix
341 progress_action("NOTHING TO PLOT.") 357 , x_predictor = predictor
342 } else { 358 , x_is_match = is_match
343 tsv <- my_cor_cov$tsv1 359 , x_algorithm = algoC
344 # tsv$salientLevel <- salient_level_lookup(tsv$featureID) 360 , x_prefix = if (pairSigFeatOnly) "Significantly contrasting features" else "Significant features"
345 # tsv$salientRCV <- salient_rcv_lookup(tsv$featureID) 361 , x_show_labels = labelFeatures
346 # tsv$salience <- salience_lookup(tsv$featureID) 362 , x_show_loado_labels = labelOrthoFeatures
347 tsv["level1Level2Sig"] <- vrbl_metadata[ match(tsv$featureID, vrbl_metadata_names), vrbl_metadata_col ] 363 , x_progress = progress_action
348 corcov_tsv_action(tsv) 364 , x_env = calc_env
349 did_plot <- TRUE 365 )
366 if ( is.null(my_cor_cov) ) {
367 progress_action("NOTHING TO PLOT.")
368 } else {
369 tsv <- my_cor_cov$tsv1
370 tsv$mz <- mz_lookup(tsv$featureID)
371 tsv$rt <- rt_lookup(tsv$featureID)
372 tsv["level1Level2Sig"] <- vrbl_metadata[ match(tsv$featureID, vrbl_metadata_names), vrbl_metadata_col ]
373 corcov_tsv_action(tsv)
374 did_plot <- TRUE
375 }
350 } 376 }
351 } 377 }
352 } 378 }
353 } else { # tesC == "none" 379 } else { # tesC == "none"
354 level_union <- unique(sort(smpl_metadata_facC)) 380 level_union <- unique(sort(smpl_metadata_facC))
355 if ( length(level_union) > 1 ) { 381 if ( length(level_union) > 1 ) {
356 ## pass 1 - contrast each selected level with all other levels combined into one "super-level" ## 382 if ( length(level_union) > 2 ) {
357 completed <- c() 383 ## pass 1 - contrast each selected level with all other levels combined into one "super-level" ##
358 lapply( 384 completed <- c()
359 X = level_union 385 lapply(
360 , FUN = function(x) { 386 X = level_union
361 fctr_lvl_1 <- x[1] 387 , FUN = function(x) {
362 fctr_lvl_2 <- { 388 fctr_lvl_1 <- x[1]
363 if ( fctr_lvl_1 %in% completed ) 389 fctr_lvl_2 <- {
364 return("DUMMY") 390 if ( fctr_lvl_1 %in% completed )
365 # strF(completed) 391 return("DUMMY")
366 completed <<- c(completed, fctr_lvl_1) 392 # strF(completed)
367 setdiff(level_union, fctr_lvl_1) 393 completed <<- c(completed, fctr_lvl_1)
394 setdiff(level_union, fctr_lvl_1)
395 }
396 chosen_samples <- smpl_metadata_facC %in% c(fctr_lvl_1, fctr_lvl_2)
397 fctr_lvl_2 <- "other"
398 # print( sprintf("sum(chosen_samples) %d, factor_level_2 %s", sum(chosen_samples), fctr_lvl_2) )
399 progress_action(sprintf("calculating/plotting contrast of %s vs. %s", fctr_lvl_1, fctr_lvl_2))
400 if (length(unique(chosen_samples)) < 1) {
401 progress_action("NOTHING TO PLOT...")
402 } else {
403 chosen_facC <- as.character(smpl_metadata_facC[chosen_samples])
404 predictor <- sapply( X = chosen_facC, FUN = function(fac) if ( fac == fctr_lvl_1 ) fctr_lvl_1 else fctr_lvl_2 )
405 my_matrix <- scdm[ chosen_samples, , drop = FALSE ]
406 # only process this column if both factors are members of lvlCSV
407 is_match <- isLevelSelected(fctr_lvl_1)
408 my_cor_cov <- do_detail_plot(
409 x_dataMatrix = my_matrix
410 , x_predictor = predictor
411 , x_is_match = is_match
412 , x_algorithm = algoC
413 , x_prefix = "Features"
414 , x_show_labels = labelFeatures
415 , x_show_loado_labels = labelOrthoFeatures
416 , x_progress = progress_action
417 , x_env = calc_env
418 )
419 if ( is.null(my_cor_cov) ) {
420 progress_action("NOTHING TO PLOT")
421 } else {
422 tsv <- my_cor_cov$tsv1
423 tsv$mz <- mz_lookup(tsv$featureID)
424 tsv$rt <- rt_lookup(tsv$featureID)
425 corcov_tsv_action(tsv)
426 did_plot <<- TRUE
427 }
428 }
429 #print("baz")
430 "dummy" # need to return a value; otherwise combn fails with an error
368 } 431 }
369 chosen_samples <- smpl_metadata_facC %in% c(fctr_lvl_1, fctr_lvl_2) 432 )
370 fctr_lvl_2 <- "other" 433 }
371 # print( sprintf("sum(chosen_samples) %d, factor_level_2 %s", sum(chosen_samples), fctr_lvl_2) )
372 progress_action(sprintf("calculating/plotting contrast of %s vs. %s", fctr_lvl_1, fctr_lvl_2))
373 if (length(unique(chosen_samples)) < 1) {
374 progress_action("NOTHING TO PLOT...")
375 } else {
376 chosen_facC <- as.character(smpl_metadata_facC[chosen_samples])
377 predictor <- sapply( X = chosen_facC, FUN = function(fac) if ( fac == fctr_lvl_1 ) fctr_lvl_1 else fctr_lvl_2 )
378 my_matrix <- scdm[ chosen_samples, , drop = FALSE ]
379 # only process this column if both factors are members of lvlCSV
380 is_match <- isLevelSelected(fctr_lvl_1)
381 my_cor_cov <- do_detail_plot(
382 x_dataMatrix = my_matrix
383 , x_predictor = predictor
384 , x_is_match = is_match
385 , x_algorithm = algoC
386 , x_prefix = "Features"
387 , x_show_labels = labelFeatures
388 , x_progress = progress_action
389 , x_env = calc_env
390 )
391 if ( is.null(my_cor_cov) ) {
392 progress_action("NOTHING TO PLOT")
393 } else {
394 tsv <- my_cor_cov$tsv1
395 # tsv$salientLevel <- salient_level_lookup(tsv$featureID)
396 # tsv$salientRCV <- salient_rcv_lookup(tsv$featureID)
397 # tsv$salience <- salience_lookup(tsv$featureID)
398 corcov_tsv_action(tsv)
399 did_plot <<- TRUE
400 }
401 }
402 #print("baz")
403 "dummy" # need to return a value; otherwise combn fails with an error
404 }
405 )
406 ## pass 2 - contrast each selected level with each of the other levels individually ## 434 ## pass 2 - contrast each selected level with each of the other levels individually ##
407 completed <- c() 435 completed <- c()
408 utils::combn( 436 utils::combn(
409 x = level_union 437 x = level_union
410 , m = 2 438 , m = 2
427 , x_predictor = predictor 455 , x_predictor = predictor
428 , x_is_match = is_match 456 , x_is_match = is_match
429 , x_algorithm = algoC 457 , x_algorithm = algoC
430 , x_prefix = "Features" 458 , x_prefix = "Features"
431 , x_show_labels = labelFeatures 459 , x_show_labels = labelFeatures
460 , x_show_loado_labels = labelOrthoFeatures
432 , x_progress = progress_action 461 , x_progress = progress_action
433 , x_env = calc_env 462 , x_env = calc_env
434 ) 463 )
435 if ( is.null(my_cor_cov) ) { 464 if ( is.null(my_cor_cov) ) {
436 progress_action("NOTHING TO PLOT") 465 progress_action("NOTHING TO PLOT")
437 } else { 466 } else {
438 tsv <- my_cor_cov$tsv1 467 tsv <- my_cor_cov$tsv1
439 # tsv$salientLevel <- salient_level_lookup(tsv$featureID) 468 tsv$mz <- mz_lookup(tsv$featureID)
440 # tsv$salientRCV <- salient_rcv_lookup(tsv$featureID) 469 tsv$rt <- rt_lookup(tsv$featureID)
441 # tsv$salience <- salience_lookup(tsv$featureID)
442 corcov_tsv_action(tsv) 470 corcov_tsv_action(tsv)
443 did_plot <<- TRUE 471 did_plot <<- TRUE
444 } 472 }
445 } 473 }
446 #print("baz") 474 #print("baz")
508 536
509 # Variant 4 of Variable Influence on Projection for OPLS from Galindo_Prieto_2014 537 # Variant 4 of Variable Influence on Projection for OPLS from Galindo_Prieto_2014
510 # Length = number of features; labels = feature identifiers. (The same is true for $correlation and $covariance.) 538 # Length = number of features; labels = feature identifiers. (The same is true for $correlation and $covariance.)
511 result$vip4p <- as.numeric(ropls_x@vipVn) 539 result$vip4p <- as.numeric(ropls_x@vipVn)
512 result$vip4o <- as.numeric(ropls_x@orthoVipVn) 540 result$vip4o <- as.numeric(ropls_x@orthoVipVn)
541 result$loadp <- as.numeric(ropls_x@loadingMN)
542 result$loado <- as.numeric(ropls_x@orthoLoadingMN)
513 # get the level names 543 # get the level names
514 level_names <- sort(levels(as.factor(ropls_x@suppLs$y))) 544 level_names <- sort(levels(as.factor(ropls_x@suppLs$y)))
545 fctr_lvl_1 <- level_names[1]
546 fctr_lvl_2 <- level_names[2]
515 feature_count <- length(ropls_x@vipVn) 547 feature_count <- length(ropls_x@vipVn)
516 result$level1 <- rep.int(x = level_names[1], times = feature_count) 548 result$level1 <- rep.int(x = fctr_lvl_1, times = feature_count)
517 result$level2 <- rep.int(x = level_names[2], times = feature_count) 549 result$level2 <- rep.int(x = fctr_lvl_2, times = feature_count)
518 # print(sprintf("sd(covariance) = %f; sd(correlation) = %f", sd(result$covariance), sd(result$correlation))) 550 # print(sprintf("sd(covariance) = %f; sd(correlation) = %f", sd(result$covariance), sd(result$correlation)))
519 superresult <- list() 551 superresult <- list()
520 if (length(result$vip4o) == 0) result$vip4o <- NA 552 if (length(result$vip4o) == 0) result$vip4o <- NA
553 greaterLevel <- sapply( X = result$correlation, FUN = function(my_corr) if ( my_corr < 0 ) fctr_lvl_1 else fctr_lvl_2 )
521 superresult$tsv1 <- data.frame( 554 superresult$tsv1 <- data.frame(
522 featureID = names(ropls_x@vipVn) 555 featureID = names(ropls_x@vipVn)
523 , factorLevel1 = result$level1 556 , factorLevel1 = result$level1
524 , factorLevel2 = result$level2 557 , factorLevel2 = result$level2
558 , greaterLevel = greaterLevel
525 , correlation = result$correlation 559 , correlation = result$correlation
526 , covariance = result$covariance 560 , covariance = result$covariance
527 , vip4p = result$vip4p 561 , vip4p = result$vip4p
528 , vip4o = result$vip4o 562 , vip4o = result$vip4o
563 , loadp = result$loadp
564 , loado = result$loado
529 , row.names = NULL 565 , row.names = NULL
530 ) 566 )
531 rownames(superresult$tsv1) <- superresult$tsv1$featureID 567 rownames(superresult$tsv1) <- superresult$tsv1$featureID
532 superresult$covariance <- result$covariance 568 superresult$covariance <- result$covariance
533 superresult$correlation <- result$correlation 569 superresult$correlation <- result$correlation
534 superresult$vip4p <- result$vip4p 570 superresult$vip4p <- result$vip4p
535 superresult$vip4o <- result$vip4o 571 superresult$vip4o <- result$vip4o
572 superresult$loadp <- result$loadp
573 superresult$loado <- result$loado
536 superresult$details <- result 574 superresult$details <- result
537 # #print(superresult$tsv1) 575 # #print(superresult$tsv1)
538 result$superresult <- superresult 576 result$superresult <- superresult
539 # Include thise in case future consumers of this routine want to use it in currently unanticipated ways 577 # Include thise in case future consumers of this routine want to use it in currently unanticipated ways
540 result$oplsda <- ropls_x 578 result$oplsda <- ropls_x