Mercurial > repos > eschen42 > w4mcorcov

comparison w4mcorcov_calc.R @ 1:0c2ad44b6c9c 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 01d4a951cf09e7b88fcec96b8043bc7568cc5c92
author eschen42
date Sun, 22 Oct 2017 18:47:57 -0400
parents 23f9fad4edfc
children e03582f26617
comparison
equal deleted inserted replaced
0:23f9fad4edfc 1:0c2ad44b6c9c
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_progress = print, x_env) {
11 off <- function(x) if (x_show_labels) x else 0 11 off <- function(x) if (x_show_labels == "0") x else 0
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.1))
62 , ylim=c(-1.0 - off(0.1), 1.0) 54 , ylim=c(-1.0 - off(0.1), 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)")
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)
74 text(x = high_x, y = 0.15, labels = fctr_lvl_2) 66 text(x = high_x, y = 0.05, labels = fctr_lvl_2)
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(head(sort(my_loado),n = n_labels))
81 , names(tail(sort(my_loadp),n = n_labels))
82 , names(tail(sort(my_loado),n = n_labels))
83 )
84 labels <- unname(sapply( X = tsv1$featureID, FUN = function(x) if( x %in% labels_to_show ) x else "" ))
76 text( 85 text(
77 y = minus_cor - 0.013 86 y = plus_cor - 0.013
78 , x = minus_cov + 0.020 87 , x = plus_cov + 0.020
79 , cex = 0.3 88 , cex = 0.3
80 , labels = tsv1$featureID 89 , labels = labels
81 , col = rgb(blue = blue, red = red, green = 0, alpha = 0.2 + 0.8 * alpha) 90 , col = rgb(blue = blue, red = red, green = 0, alpha = 0.2 + 0.8 * alpha)
82 , srt = -30 # slant 30 degrees downward 91 , srt = -30 # slant 30 degrees downward
83 , adj = 0 # left-justified 92 , adj = 0 # left-justified
84 ) 93 )
85 } 94 }
160 if (!(facC %in% names(smpl_metadata))) { 169 if (!(facC %in% names(smpl_metadata))) {
161 failure_action(sprintf("bad parameter! Factor name '%s' not found in sampleMetadata", facC)) 170 failure_action(sprintf("bad parameter! Factor name '%s' not found in sampleMetadata", facC))
162 return ( FALSE ) 171 return ( FALSE )
163 } 172 }
164 173
174 mz <- vrbl_metadata$mz
175 names(mz) <- vrbl_metadata$variableMetadata
176 mz_lookup <- function(feature) unname(mz[feature])
177
178 rt <- vrbl_metadata$rt
179 names(rt) <- vrbl_metadata$variableMetadata
180 rt_lookup <- function(feature) unname(rt[feature])
181
165 # calculate salience_df as data.frame(feature, max_level, max_median, max_rcv, mean_median, salience, salient_rcv) 182 # 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( 183 salience_df <- calc_env$salience_df <- w4msalience(
167 data_matrix = data_matrix 184 data_matrix = data_matrix
168 , sample_class = smpl_metadata[,facC] 185 , sample_class = smpl_metadata[,facC]
169 , failure_action = failure_action 186 , failure_action = failure_action
172 my_df <- data.frame( 189 my_df <- data.frame(
173 featureID = salience_df$feature 190 featureID = salience_df$feature
174 , salientLevel = salience_df$max_level 191 , salientLevel = salience_df$max_level
175 , salientRCV = salience_df$salient_rcv 192 , salientRCV = salience_df$salient_rcv
176 , salience = salience_df$salience 193 , salience = salience_df$salience
194 , mz = mz_lookup(salience_df$feature)
195 , rt = rt_lookup(salience_df$feature)
177 ) 196 )
178 my_df[order(-my_df$salience),] 197 my_df[order(-my_df$salience),]
179 }) 198 })
180 salience <- salience_df$salience 199
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 200 # transform wildcards to regexen
191 if (matchingC == "wildcard") { 201 if (matchingC == "wildcard") {
192 # strsplit(x = "hello,wild,world", split = ",") 202 # strsplit(x = "hello,wild,world", split = ",")
193 levCSV <- gsub("[.]", "[.]", levCSV) 203 levCSV <- gsub("[.]", "[.]", levCSV)
194 levCSV <- strsplit(x = levCSV, split = ",") 204 levCSV <- strsplit(x = levCSV, split = ",")
269 level_union <- c(level_union, col_match[2], col_match[3]) 279 level_union <- c(level_union, col_match[2], col_match[3])
270 } 280 }
271 } 281 }
272 } 282 }
273 level_union <- unique(sort(level_union)) 283 level_union <- unique(sort(level_union))
274 overall_significant <- 1 == vrbl_metadata[,intersample_sig_col] 284 overall_significant <- 1 == ( if (intersample_sig_col %in% colnames(vrbl_metadata)) vrbl_metadata[,intersample_sig_col] else TRUE )
275 if ( length(level_union) > 1 ) { 285 if ( length(level_union) > 2 ) {
276 chosen_samples <- smpl_metadata_facC %in% level_union 286 chosen_samples <- smpl_metadata_facC %in% level_union
277 chosen_facC <- as.character(smpl_metadata_facC[chosen_samples]) 287 chosen_facC <- as.character(smpl_metadata_facC[chosen_samples])
278 col_selector <- vrbl_metadata_names[ overall_significant ] 288 col_selector <- vrbl_metadata_names[ overall_significant ]
279 my_matrix <- scdm[ chosen_samples, col_selector, drop = FALSE ] 289 my_matrix <- scdm[ chosen_samples, col_selector, drop = FALSE ]
280 fctr_lvl_2 <- "other" 290 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)) 291 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 ) 292 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( 293 my_cor_cov <- do_detail_plot(
285 x_dataMatrix = my_matrix 294 x_dataMatrix = my_matrix
286 , x_predictor = predictor 295 , x_predictor = predictor
292 , x_env = calc_env 301 , x_env = calc_env
293 ) 302 )
294 if ( is.null(my_cor_cov) ) { 303 if ( is.null(my_cor_cov) ) {
295 progress_action("NOTHING TO PLOT.") 304 progress_action("NOTHING TO PLOT.")
296 } else { 305 } else {
297 tsv <- my_cor_cov$tsv1 306 my_tsv <- my_cor_cov$tsv1
298 # tsv$salientLevel <- salient_level_lookup(tsv$featureID) 307 my_tsv$mz <- mz_lookup(my_tsv$featureID)
299 # tsv$salientRCV <- salient_rcv_lookup(tsv$featureID) 308 my_tsv$rt <- rt_lookup(my_tsv$featureID)
300 # tsv$salience <- salience_lookup(tsv$featureID) 309 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 ] 310 tsv <<- my_tsv
302 corcov_tsv_action(tsv) 311 corcov_tsv_action(tsv)
303 did_plot <- TRUE 312 did_plot <<- TRUE
304 } 313 }
305 } 314 }
306 } 315 if ( length(level_union) != 2 ) {
307 316 fctr_lvl_2 <- "other"
308 ## next, contrast each selected level with each of the other levels individually ## 317 for ( fctr_lvl_1 in level_union[1:length(level_union)] ) {
309 # process columns matching the pattern 318 plot_action(fctr_lvl_1, fctr_lvl_2)
310 for (i in 1:length(col_matches)) { 319 }
311 # for each potential match of the pattern 320 } else {
312 col_match <- col_matches[[i]] 321 plot_action(fctr_lvl_1 = level_union[1], fctr_lvl_2 = level_union[2])
313 if (length(col_match) > 0) { 322 }
314 # it's an actual match; extract the pieces, e.g., k10_kruskal_k4.k3_sig 323 }
315 vrbl_metadata_col <- col_match[1] # ^^^^^^^^^^^^^^^^^^^^^ # Column name 324
316 fctr_lvl_1 <- col_match[2] # ^^ # Factor-level 1 325 if ( length(level_union) > 1 ) {
317 fctr_lvl_2 <- col_match[3] # ^^ # Factor-level 2 326 ## next, contrast each selected level with each of the other levels individually ##
318 # only process this column if both factors are members of lvlCSV 327 # process columns matching the pattern
319 is_match <- isLevelSelected(fctr_lvl_1) && isLevelSelected(fctr_lvl_2) 328 for (i in 1:length(col_matches)) {
320 progress_action(sprintf("calculating/plotting contrast of %s vs. %s", fctr_lvl_1, fctr_lvl_2)) 329 # for each potential match of the pattern
321 # TODO delete next line displaying currently-processed column 330 col_match <- col_matches[[i]]
322 # cat(sprintf("%s %s %s %s\n", vrbl_metadata_col, fctr_lvl_1, fctr_lvl_2, is_match)) 331 if (length(col_match) > 0) {
323 # choose only samples with one of the two factors for this column 332 # 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) 333 vrbl_metadata_col <- col_match[1] # ^^^^^^^^^^^^^^^^^^^^^ # Column name
325 predictor <- smpl_metadata_facC[chosen_samples] 334 fctr_lvl_1 <- col_match[2] # ^^ # Factor-level 1
326 # extract only the significantly-varying features and the chosen samples 335 fctr_lvl_2 <- col_match[3] # ^^ # Factor-level 2
327 fully_significant <- 1 == vrbl_metadata[,vrbl_metadata_col] * vrbl_metadata[,intersample_sig_col] 336 # only process this column if both factors are members of lvlCSV
328 col_selector <- vrbl_metadata_names[ if ( pairSigFeatOnly ) fully_significant else overall_significant ] 337 is_match <- isLevelSelected(fctr_lvl_1) && isLevelSelected(fctr_lvl_2)
329 my_matrix <- scdm[ chosen_samples, col_selector, drop = FALSE ] 338 progress_action(sprintf("calculating/plotting contrast of %s vs. %s", fctr_lvl_1, fctr_lvl_2))
330 my_cor_cov <- do_detail_plot( 339 # TODO delete next line displaying currently-processed column
331 x_dataMatrix = my_matrix 340 # cat(sprintf("%s %s %s %s\n", vrbl_metadata_col, fctr_lvl_1, fctr_lvl_2, is_match))
332 , x_predictor = predictor 341 # choose only samples with one of the two factors for this column
333 , x_is_match = is_match 342 chosen_samples <- smpl_metadata_facC %in% c(fctr_lvl_1, fctr_lvl_2)
334 , x_algorithm = algoC 343 predictor <- smpl_metadata_facC[chosen_samples]
335 , x_prefix = if (pairSigFeatOnly) "Significantly contrasting features" else "Significant features" 344 # extract only the significantly-varying features and the chosen samples
336 , x_show_labels = labelFeatures 345 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 346 col_selector <- vrbl_metadata_names[ if ( pairSigFeatOnly ) fully_significant else overall_significant ]
338 , x_env = calc_env 347 my_matrix <- scdm[ chosen_samples, col_selector, drop = FALSE ]
339 ) 348 my_cor_cov <- do_detail_plot(
340 if ( is.null(my_cor_cov) ) { 349 x_dataMatrix = my_matrix
341 progress_action("NOTHING TO PLOT.") 350 , x_predictor = predictor
342 } else { 351 , x_is_match = is_match
343 tsv <- my_cor_cov$tsv1 352 , x_algorithm = algoC
344 # tsv$salientLevel <- salient_level_lookup(tsv$featureID) 353 , x_prefix = if (pairSigFeatOnly) "Significantly contrasting features" else "Significant features"
345 # tsv$salientRCV <- salient_rcv_lookup(tsv$featureID) 354 , x_show_labels = labelFeatures
346 # tsv$salience <- salience_lookup(tsv$featureID) 355 , x_progress = progress_action
347 tsv["level1Level2Sig"] <- vrbl_metadata[ match(tsv$featureID, vrbl_metadata_names), vrbl_metadata_col ] 356 , x_env = calc_env
348 corcov_tsv_action(tsv) 357 )
349 did_plot <- TRUE 358 if ( is.null(my_cor_cov) ) {
359 progress_action("NOTHING TO PLOT.")
360 } else {
361 tsv <- my_cor_cov$tsv1
362 tsv$mz <- mz_lookup(tsv$featureID)
363 tsv$rt <- rt_lookup(tsv$featureID)
364 tsv["level1Level2Sig"] <- vrbl_metadata[ match(tsv$featureID, vrbl_metadata_names), vrbl_metadata_col ]
365 corcov_tsv_action(tsv)
366 did_plot <- TRUE
367 }
350 } 368 }
351 } 369 }
352 } 370 }
353 } else { # tesC == "none" 371 } else { # tesC == "none"
354 level_union <- unique(sort(smpl_metadata_facC)) 372 level_union <- unique(sort(smpl_metadata_facC))
355 if ( length(level_union) > 1 ) { 373 if ( length(level_union) > 1 ) {
356 ## pass 1 - contrast each selected level with all other levels combined into one "super-level" ## 374 if ( length(level_union) > 2 ) {
357 completed <- c() 375 ## pass 1 - contrast each selected level with all other levels combined into one "super-level" ##
358 lapply( 376 completed <- c()
359 X = level_union 377 lapply(
360 , FUN = function(x) { 378 X = level_union
361 fctr_lvl_1 <- x[1] 379 , FUN = function(x) {
362 fctr_lvl_2 <- { 380 fctr_lvl_1 <- x[1]
363 if ( fctr_lvl_1 %in% completed ) 381 fctr_lvl_2 <- {
364 return("DUMMY") 382 if ( fctr_lvl_1 %in% completed )
365 # strF(completed) 383 return("DUMMY")
366 completed <<- c(completed, fctr_lvl_1) 384 # strF(completed)
367 setdiff(level_union, fctr_lvl_1) 385 completed <<- c(completed, fctr_lvl_1)
386 setdiff(level_union, fctr_lvl_1)
387 }
388 chosen_samples <- smpl_metadata_facC %in% c(fctr_lvl_1, fctr_lvl_2)
389 fctr_lvl_2 <- "other"
390 # print( sprintf("sum(chosen_samples) %d, factor_level_2 %s", sum(chosen_samples), fctr_lvl_2) )
391 progress_action(sprintf("calculating/plotting contrast of %s vs. %s", fctr_lvl_1, fctr_lvl_2))
392 if (length(unique(chosen_samples)) < 1) {
393 progress_action("NOTHING TO PLOT...")
394 } else {
395 chosen_facC <- as.character(smpl_metadata_facC[chosen_samples])
396 predictor <- sapply( X = chosen_facC, FUN = function(fac) if ( fac == fctr_lvl_1 ) fctr_lvl_1 else fctr_lvl_2 )
397 my_matrix <- scdm[ chosen_samples, , drop = FALSE ]
398 # only process this column if both factors are members of lvlCSV
399 is_match <- isLevelSelected(fctr_lvl_1)
400 my_cor_cov <- do_detail_plot(
401 x_dataMatrix = my_matrix
402 , x_predictor = predictor
403 , x_is_match = is_match
404 , x_algorithm = algoC
405 , x_prefix = "Features"
406 , x_show_labels = labelFeatures
407 , x_progress = progress_action
408 , x_env = calc_env
409 )
410 if ( is.null(my_cor_cov) ) {
411 progress_action("NOTHING TO PLOT")
412 } else {
413 tsv <- my_cor_cov$tsv1
414 tsv$mz <- mz_lookup(tsv$featureID)
415 tsv$rt <- rt_lookup(tsv$featureID)
416 corcov_tsv_action(tsv)
417 did_plot <<- TRUE
418 }
419 }
420 #print("baz")
421 "dummy" # need to return a value; otherwise combn fails with an error
368 } 422 }
369 chosen_samples <- smpl_metadata_facC %in% c(fctr_lvl_1, fctr_lvl_2) 423 )
370 fctr_lvl_2 <- "other" 424 }
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 ## 425 ## pass 2 - contrast each selected level with each of the other levels individually ##
407 completed <- c() 426 completed <- c()
408 utils::combn( 427 utils::combn(
409 x = level_union 428 x = level_union
410 , m = 2 429 , m = 2
434 ) 453 )
435 if ( is.null(my_cor_cov) ) { 454 if ( is.null(my_cor_cov) ) {
436 progress_action("NOTHING TO PLOT") 455 progress_action("NOTHING TO PLOT")
437 } else { 456 } else {
438 tsv <- my_cor_cov$tsv1 457 tsv <- my_cor_cov$tsv1
439 # tsv$salientLevel <- salient_level_lookup(tsv$featureID) 458 tsv$mz <- mz_lookup(tsv$featureID)
440 # tsv$salientRCV <- salient_rcv_lookup(tsv$featureID) 459 tsv$rt <- rt_lookup(tsv$featureID)
441 # tsv$salience <- salience_lookup(tsv$featureID)
442 corcov_tsv_action(tsv) 460 corcov_tsv_action(tsv)
443 did_plot <<- TRUE 461 did_plot <<- TRUE
444 } 462 }
445 } 463 }
446 #print("baz") 464 #print("baz")
508 526
509 # Variant 4 of Variable Influence on Projection for OPLS from Galindo_Prieto_2014 527 # 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.) 528 # Length = number of features; labels = feature identifiers. (The same is true for $correlation and $covariance.)
511 result$vip4p <- as.numeric(ropls_x@vipVn) 529 result$vip4p <- as.numeric(ropls_x@vipVn)
512 result$vip4o <- as.numeric(ropls_x@orthoVipVn) 530 result$vip4o <- as.numeric(ropls_x@orthoVipVn)
531 result$loadp <- as.numeric(ropls_x@loadingMN)
532 result$loado <- as.numeric(ropls_x@orthoLoadingMN)
513 # get the level names 533 # get the level names
514 level_names <- sort(levels(as.factor(ropls_x@suppLs$y))) 534 level_names <- sort(levels(as.factor(ropls_x@suppLs$y)))
535 fctr_lvl_1 <- level_names[1]
536 fctr_lvl_2 <- level_names[2]
515 feature_count <- length(ropls_x@vipVn) 537 feature_count <- length(ropls_x@vipVn)
516 result$level1 <- rep.int(x = level_names[1], times = feature_count) 538 result$level1 <- rep.int(x = fctr_lvl_1, times = feature_count)
517 result$level2 <- rep.int(x = level_names[2], times = feature_count) 539 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))) 540 # print(sprintf("sd(covariance) = %f; sd(correlation) = %f", sd(result$covariance), sd(result$correlation)))
519 superresult <- list() 541 superresult <- list()
520 if (length(result$vip4o) == 0) result$vip4o <- NA 542 if (length(result$vip4o) == 0) result$vip4o <- NA
543 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( 544 superresult$tsv1 <- data.frame(
522 featureID = names(ropls_x@vipVn) 545 featureID = names(ropls_x@vipVn)
523 , factorLevel1 = result$level1 546 , factorLevel1 = result$level1
524 , factorLevel2 = result$level2 547 , factorLevel2 = result$level2
548 , greaterLevel = greaterLevel
525 , correlation = result$correlation 549 , correlation = result$correlation
526 , covariance = result$covariance 550 , covariance = result$covariance
527 , vip4p = result$vip4p 551 , vip4p = result$vip4p
528 , vip4o = result$vip4o 552 , vip4o = result$vip4o
553 , loadp = result$loadp
554 , loado = result$loado
529 , row.names = NULL 555 , row.names = NULL
530 ) 556 )
531 rownames(superresult$tsv1) <- superresult$tsv1$featureID 557 rownames(superresult$tsv1) <- superresult$tsv1$featureID
532 superresult$covariance <- result$covariance 558 superresult$covariance <- result$covariance
533 superresult$correlation <- result$correlation 559 superresult$correlation <- result$correlation
534 superresult$vip4p <- result$vip4p 560 superresult$vip4p <- result$vip4p
535 superresult$vip4o <- result$vip4o 561 superresult$vip4o <- result$vip4o
562 superresult$loadp <- result$loadp
563 superresult$loado <- result$loado
536 superresult$details <- result 564 superresult$details <- result
537 # #print(superresult$tsv1) 565 # #print(superresult$tsv1)
538 result$superresult <- superresult 566 result$superresult <- superresult
539 # Include thise in case future consumers of this routine want to use it in currently unanticipated ways 567 # Include thise in case future consumers of this routine want to use it in currently unanticipated ways
540 result$oplsda <- ropls_x 568 result$oplsda <- ropls_x