Mercurial > repos > eschen42 > w4mcorcov

comparison w4mcorcov_calc.R @ 2:e03582f26617 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 7682e8e7ae2bfb926d94b414b9a1649389f33582
author eschen42
date Sun, 12 Nov 2017 19:45:36 -0500
parents 0c2ad44b6c9c
children 5aaab36bc523
comparison
equal deleted inserted replaced
1:0c2ad44b6c9c 2:e03582f26617
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 == "0") x else 0 11 off <- function(x) if (x_show_labels == "0") 0 else x
12 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) {
13 my_oplsda <- opls( 13 my_oplsda <- opls(
14 x = x_dataMatrix 14 x = x_dataMatrix
15 , y = x_predictor 15 , y = x_predictor
16 , algoC = x_algorithm 16 , algoC = x_algorithm
32 min_x <- min(covariance) 32 min_x <- min(covariance)
33 max_x <- max(covariance) 33 max_x <- max(covariance)
34 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))
35 covariance <- covariance / lim_x 35 covariance <- covariance / lim_x
36 lim_x <- 1.2 36 lim_x <- 1.2
37 main_label <- sprintf("%s for levels %s versus %s", x_prefix, fctr_lvl_1, fctr_lvl_2) 37 main_label <- sprintf("%s for level %s versus %s", x_prefix, fctr_lvl_1, fctr_lvl_2)
38 main_cex <- min(1.0, 46.0/nchar(main_label)) 38 main_cex <- min(1.0, 46.0/nchar(main_label))
39 # "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],
40 # 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]."
41 # (Mehmood 2012 doi:10.1186/1748-7188-6-27) 41 # (Mehmood 2012 doi:10.1186/1748-7188-6-27)
42 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)))
48 cex <- 0.75 48 cex <- 0.75
49 plot( 49 plot(
50 y = plus_cor 50 y = plus_cor
51 , x = plus_cov 51 , x = plus_cov
52 , type="p" 52 , type="p"
53 , xlim=c(-lim_x, lim_x + off(0.1)) 53 , xlim=c( -lim_x - off(0.2), lim_x + off(0.2) )
54 , ylim=c(-1.0 - off(0.1), 1.0) 54 , ylim=c( -1.0 - off(0.2), 1.0 + off(0.2) )
55 , xlab = sprintf("relative covariance(feature,t1)") 55 , xlab = sprintf("relative covariance(feature,t1)")
56 , ylab = sprintf("correlation(feature,t1)") 56 , ylab = sprintf("correlation(feature,t1)")
57 , main = main_label 57 , main = main_label
58 , cex.main = main_cex 58 , cex.main = main_cex
59 , cex = cex 59 , cex = cex
60 , pch = 16 60 , pch = 16
61 , col = rgb(blue = blue, red = red, green = 0, alpha = alpha) 61 , col = rgb(blue = blue, red = red, green = 0, alpha = alpha)
62 ) 62 )
63 low_x <- -0.7 * lim_x 63 low_x <- -0.7 * lim_x
64 high_x <- 0.7 * lim_x 64 high_x <- 0.7 * lim_x
65 text(x = low_x, y = -0.05, labels = fctr_lvl_1) 65 text(x = low_x, y = -0.05, labels = fctr_lvl_1, col = "blue")
66 text(x = high_x, y = 0.05, labels = fctr_lvl_2) 66 text(x = high_x, y = 0.05, labels = fctr_lvl_2, col = "red")
67 if ( x_show_labels != "0" ) { 67 if ( x_show_labels != "0" ) {
68 my_loadp <- loadp 68 my_loadp <- loadp
69 my_loado <- loado 69 my_loado <- loado
70 names(my_loadp) <- tsv1$featureID 70 names(my_loadp) <- tsv1$featureID
71 names(my_loado) <- tsv1$featureID 71 names(my_loado) <- tsv1$featureID
75 n_labels <- as.numeric(x_show_labels) 75 n_labels <- as.numeric(x_show_labels)
76 } 76 }
77 n_labels <- min( n_labels, (1 + length(loadp)) / 2 ) 77 n_labels <- min( n_labels, (1 + length(loadp)) / 2 )
78 labels_to_show <- c( 78 labels_to_show <- c(
79 names(head(sort(my_loadp),n = n_labels)) 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)) 80 , names(tail(sort(my_loadp),n = n_labels))
82 , names(tail(sort(my_loado),n = n_labels))
83 ) 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 }
84 labels <- unname(sapply( X = tsv1$featureID, FUN = function(x) if( x %in% labels_to_show ) x else "" )) 89 labels <- unname(sapply( X = tsv1$featureID, FUN = function(x) if( x %in% labels_to_show ) x else "" ))
85 text( 90 text(
86 y = plus_cor - 0.013 91 y = plus_cor - 0.013
87 , x = plus_cov + 0.020 92 , x = plus_cov + 0.020
88 , cex = 0.3 93 , cex = 0.4
89 , labels = labels 94 , labels = labels
90 , col = rgb(blue = blue, red = red, green = 0, alpha = 0.2 + 0.8 * alpha) 95 , col = rgb(blue = 0, red = 0, green = 0, alpha = 0.5) # rgb(blue = blue, red = red, green = 0, alpha = 0.2 + 0.8 * alpha)
91 , srt = -30 # slant 30 degrees downward 96 , srt = -30 # slant 30 degrees downward
92 , adj = 0 # left-justified 97 , adj = 0 # left-justified
93 ) 98 )
94 } 99 }
95 } 100 }
162 # extract the levels from the environment 167 # extract the levels from the environment
163 originalLevCSV <- levCSV <- calc_env$levCSV 168 originalLevCSV <- levCSV <- calc_env$levCSV
164 # matchingC is one of { "none", "wildcard", "regex" } 169 # matchingC is one of { "none", "wildcard", "regex" }
165 matchingC <- calc_env$matchingC 170 matchingC <- calc_env$matchingC
166 labelFeatures <- calc_env$labelFeatures 171 labelFeatures <- calc_env$labelFeatures
172 labelOrthoFeatures <- calc_env$labelOrthoFeatures
167 173
168 # arg/env checking 174 # arg/env checking
169 if (!(facC %in% names(smpl_metadata))) { 175 if (!(facC %in% names(smpl_metadata))) {
170 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))
171 return ( FALSE ) 177 return ( FALSE )
295 , x_predictor = predictor 301 , x_predictor = predictor
296 , x_is_match = is_match 302 , x_is_match = is_match
297 , x_algorithm = algoC 303 , x_algorithm = algoC
298 , x_prefix = if (pairSigFeatOnly) "Significantly contrasting features" else "Significant features" 304 , x_prefix = if (pairSigFeatOnly) "Significantly contrasting features" else "Significant features"
299 , x_show_labels = labelFeatures 305 , x_show_labels = labelFeatures
306 , x_show_loado_labels = labelOrthoFeatures
300 , x_progress = progress_action 307 , x_progress = progress_action
301 , x_env = calc_env 308 , x_env = calc_env
302 ) 309 )
303 if ( is.null(my_cor_cov) ) { 310 if ( is.null(my_cor_cov) ) {
304 progress_action("NOTHING TO PLOT.") 311 progress_action("NOTHING TO PLOT.")
350 , x_predictor = predictor 357 , x_predictor = predictor
351 , x_is_match = is_match 358 , x_is_match = is_match
352 , x_algorithm = algoC 359 , x_algorithm = algoC
353 , x_prefix = if (pairSigFeatOnly) "Significantly contrasting features" else "Significant features" 360 , x_prefix = if (pairSigFeatOnly) "Significantly contrasting features" else "Significant features"
354 , x_show_labels = labelFeatures 361 , x_show_labels = labelFeatures
362 , x_show_loado_labels = labelOrthoFeatures
355 , x_progress = progress_action 363 , x_progress = progress_action
356 , x_env = calc_env 364 , x_env = calc_env
357 ) 365 )
358 if ( is.null(my_cor_cov) ) { 366 if ( is.null(my_cor_cov) ) {
359 progress_action("NOTHING TO PLOT.") 367 progress_action("NOTHING TO PLOT.")
402 , x_predictor = predictor 410 , x_predictor = predictor
403 , x_is_match = is_match 411 , x_is_match = is_match
404 , x_algorithm = algoC 412 , x_algorithm = algoC
405 , x_prefix = "Features" 413 , x_prefix = "Features"
406 , x_show_labels = labelFeatures 414 , x_show_labels = labelFeatures
415 , x_show_loado_labels = labelOrthoFeatures
407 , x_progress = progress_action 416 , x_progress = progress_action
408 , x_env = calc_env 417 , x_env = calc_env
409 ) 418 )
410 if ( is.null(my_cor_cov) ) { 419 if ( is.null(my_cor_cov) ) {
411 progress_action("NOTHING TO PLOT") 420 progress_action("NOTHING TO PLOT")
446 , x_predictor = predictor 455 , x_predictor = predictor
447 , x_is_match = is_match 456 , x_is_match = is_match
448 , x_algorithm = algoC 457 , x_algorithm = algoC
449 , x_prefix = "Features" 458 , x_prefix = "Features"
450 , x_show_labels = labelFeatures 459 , x_show_labels = labelFeatures
460 , x_show_loado_labels = labelOrthoFeatures
451 , x_progress = progress_action 461 , x_progress = progress_action
452 , x_env = calc_env 462 , x_env = calc_env
453 ) 463 )
454 if ( is.null(my_cor_cov) ) { 464 if ( is.null(my_cor_cov) ) {
455 progress_action("NOTHING TO PLOT") 465 progress_action("NOTHING TO PLOT")