| Previous changeset 5:50f60f94c034 (2018-03-30) Next changeset 7:066b1f409e9f (2018-07-26) |
|
Commit message:
planemo upload for repository https://github.com/HegemanLab/w4mcorcov_galaxy_wrapper/tree/master commit cafda5095a79ce2376325b57337302f95137195d |
|
modified:
w4mcorcov.xml w4mcorcov_calc.R w4mcorcov_wrapper.R |
| b |
| diff -r 50f60f94c034 -r 7bd523ca1f9a w4mcorcov.xml --- a/w4mcorcov.xml Fri Mar 30 14:59:19 2018 -0400 +++ b/w4mcorcov.xml Wed Jul 18 12:35:55 2018 -0400 |
| [ |
| b'@@ -1,7 +1,7 @@\n-\xef\xbb\xbf<tool id="w4mcorcov" name="OPLS-DA_Contrasts" version="0.98.8">\n+\xef\xbb\xbf<tool id="w4mcorcov" name="OPLS-DA_Contrasts" version="0.98.10">\n \n <description>OPLS-DA Contrasts of Univariate Results</description>\n- \n+\n <macros>\n <xml name="paramPairSigFeatOnly">\n \t\t\t<param\n@@ -13,13 +13,35 @@\n \t\t\t\tlabel="Retain only pairwise-significant features"\n \t\t\t\thelp="When this option is set to \'Yes\', analysis will be performed including only features that differ significantly for the pair of levels being contrasted; when set to \'No\', any feature that varies significantly across all levels will be included (i.e., exclude any feature that is not significantly different across all levels). See examples below." />\n </xml>\n+ <xml name="cplots">\n+\t\t\t<param name="cplot_y" label="C-plot Y-axis" type="select" help="Choose the Y-axis for C-plots.">\n+\t\t\t\t<option value="correlation">Plot VIP versus correlation</option>\n+\t\t\t\t<option value="covariance">Plot VIP versus covariance</option>\n+\t\t\t</param>\n+\t\t\t<param\n+\t\t\t\tname="cplot_p"\n+\t\t\t\ttype="boolean"\n+\t\t\t\tchecked="TRUE"\n+\t\t\t\ttruevalue="TRUE"\n+\t\t\t\tfalsevalue="FALSE"\n+\t\t\t\tlabel="Produce predictor C-plot"\n+\t\t\t\thelp="When this option is set to \'Yes\', correlation will be plotted against vip4 for predictor loadings." />\n+\t\t\t<param\n+\t\t\t\tname="cplot_o"\n+\t\t\t\ttype="boolean"\n+\t\t\t\tchecked="TRUE"\n+\t\t\t\ttruevalue="TRUE"\n+\t\t\t\tfalsevalue="FALSE"\n+\t\t\t\tlabel="Produce orthogonal C-plot"\n+\t\t\t\thelp="When this option is set to \'Yes\', correlation will be plotted against vip4 for orthogonal loadings." />\n+ </xml>\n \t</macros>\n \n <requirements>\n <requirement type="package">r-batch</requirement>\n <requirement type="package">bioconductor-ropls</requirement>\n </requirements>\n- \n+\n <stdio>\n <exit_code range="1:" level="fatal" />\n </stdio>\n@@ -40,6 +62,15 @@\n levCSV \'$levCSV\'\n matchingC \'$matchingC\'\n labelFeatures \'$labelFeatures\'\n+ #if str( $xplots.expPlot ) == "none":\n+ cplot_p "FALSE"\n+ cplot_o "FALSE"\n+ cplot_y "correlation"\n+ #else if str( $xplots.expPlot ) == "cplot":\n+ cplot_p "$xplots.cplot_p"\n+ cplot_o "$xplots.cplot_o"\n+ cplot_y "$xplots.cplot_y"\n+ #end if\n contrast_detail \'$contrast_detail\'\n contrast_corcov \'$contrast_corcov\'\n contrast_salience \'$contrast_salience\'\n@@ -55,7 +86,7 @@\n \t\t\t\t<option value="none">none - Display all features from variableMetadata (rather than choosing a subset based on significance in univariate testing)</option>\n \t\t\t\t<option value="ttest">ttest - Student\'s t-test (parametric test, qualitative factor with exactly 2 levels)</option>\n \t\t\t\t<option value="anova">anova - Analysis of variance (parametric test, qualitative factor with more than 2 levels)</option>\n-\t\t\t\t<option value="wilcoxon">wilcoxon - Wilcoxon rank test (nonparametric test, qualitative factor with exactly 2 levels)</option> \n+\t\t\t\t<option value="wilcoxon">wilcoxon - Wilcoxon rank test (nonparametric test, qualitative factor with exactly 2 levels)</option>\n \t\t\t\t<option value="kruskal">kruskal - Kruskal-Wallis rank test (nonparametric test, qualitative factor with more than 2 levels)</option>\n \t\t\t</param>\n \t\t\t<when value="none" />\n@@ -104,6 +135,16 @@\n <option value="regex">use regular expressions for matching level-names</option>\n </param>\n <param name="labelFeatures" type="text" value="3" label="How many features having extreme loadings should be labelled on cov-vs.-cor plot" help="Specify the number of features at each of the loading-extremes that should be labelled (with the name of the feature) on the covariance-vs.-correlation plot; specify \'ALL\' to label all features or \'0\' to label no features; this choice has no effect on the OPLS-DA loadings plot."/>\n+ <conditional name="xplots">\n+\t\t\t<param name="expPlot" label="Extra plots to include" type="select" help="Choosing \'none\' hides further choices.">\n+\t\t\t\t<option value="none">none - Do not include additonal extra plots.</option>\n+\t\t\t\t<option val'..b'ion coefficient for the fitted data; black bars indicate performance in cross-validation tests (Th]]>é<![CDATA[venot, 2017)\n-- (top-right) OPLS-DA **scores-plot** for the two projections (Th]]>é<![CDATA[venot *et al.*, 2015)\n-- (bottom-right) OPLS-DA **loadings-plot** for the two projections (*ibid.*)\n+- (first row, left) **correlation-versus-covariance plot** of OPLS-DA results (a work-alike for the S-PLOT, computed using formula in Supplement to Wiklund, *op. cit.*); point-color becomes saturated as the "variable importance in projection to the predictive components" (VIP\\ :subscript:`4,p` from Galindo-Prieto *et al.* 2014) ranges from 0.83 and 1.21 (Mehmood *et al.* 2012), for use to identify features for consideration as biomarkers.\n+- (second row, left) **model-overview plot** for the two projections; grey bars are the correlation coefficient for the fitted data; black bars indicate performance in cross-validation tests (Th]]>é<![CDATA[venot, 2017)\n+- (first row, right) OPLS-DA **scores-plot** for the two projections (Th]]>é<![CDATA[venot *et al.*, 2015)\n+- (second row, right) **correlation-versus-covariance plot** of OPLS-DA results **orthogonal to the predictor** (see section "S-Plot of Orthogonal Component" in Wiklund, *op. cit.*, pp. 120-121; this characterizes features with the greatest variation independent of the predictor).\n+- (third row, left, when "**predictor C-plot**" is chosen under "Extra plots to include") plot of the covariance or correlation vs. the VIP for the projection *parallel to the predictor*, for use to identify features for consideration as biomarkers.\n+- (third row, right, when "**orthogonal C-plot**" is chosen under "Extra plots to include") plotof the covariance or correlation vs. the VIP for the projection *orthogonal to the predictor*, for use to identify features varying considerably without regard to the predictor.\n \n [OUT] Contrast Correlation-Covarinace data TABULAR\n | A tab-separated values file of metadata for each feature for each contrast in which it was included.\n@@ -475,7 +558,7 @@\n - **level1Level2Sig** - (Only present when a test other than "none" is chosen) \'1\' when feature varies significantly across all classes (i.e., not pair-wise); \'0\' otherwise\n \n [OUT] Feature "Salience" data TABULAR\n- | Metrics for the "salient level" for each feature, i.e., the level at which the feature is more prominent than any other level. This is *not* at all related to the SIMCA OPLS-DA S-PLOT; rather, it is intended as a potential (and unproven) way to identify features that may suggest potential biomarkers without dimensional reduction of data. This is a tab-separated values file having the following columns:\n+ | Metrics for the "salient level" for each feature, i.e., the level at which the feature is more prominent than any other level. This is *not* at all related to the SIMCA OPLS-DA S-PLOT; rather, it is intended as a potential way to discover features for consideration as potential biomarkers without dimensionally reducting the data. This is a tab-separated values file having the following columns:\n \n - **featureID** - feature identifier\n - **salientLevel** - salient level, i.e., for the feature, the class-level having the greatest median intensity\n@@ -651,6 +734,15 @@\n Release notes\n -------------\n \n+0.98.10\n+\n+- new feature: C-plots of VIP versus correlation or relative covariance.\n+- bug fix: Handle issue 2 - features now are only pareto-scaled per Wikland *op cit.*.\n+\n+0.98.9\n+\n+- bug fix: Handle issue 1 - handle features removed by ropls because variance is less than 2.2e-16.\n+\n 0.98.8\n \n - new feature: Replace loadings plot with correlation-versus-covariance plot for orthogonal features, i.e., the consistency of features influencing within-treatment variation (which is linearly related to the loading of the orthogonal projection) versus consistency. This eliminates the need for the parameter to suppress labels for features with extreme orthogonal loadings\n' |
| b |
| diff -r 50f60f94c034 -r 7bd523ca1f9a w4mcorcov_calc.R --- a/w4mcorcov_calc.R Fri Mar 30 14:59:19 2018 -0400 +++ b/w4mcorcov_calc.R Wed Jul 18 12:35:55 2018 -0400 |
| [ |
| b'@@ -7,9 +7,23 @@\n #### OPLS-DA\n algoC <- "nipals"\n \n-do_detail_plot <- function(x_dataMatrix, x_predictor, x_is_match, x_algorithm, x_prefix, x_show_labels, x_progress = print, x_env, x_crossval_i) {\n+do_detail_plot <- function(\n+ x_dataMatrix\n+, x_predictor\n+, x_is_match\n+, x_algorithm\n+, x_prefix\n+, x_show_labels\n+, x_progress = print\n+, x_env\n+, x_crossval_i\n+) {\n off <- function(x) if (x_show_labels == "0") 0 else x\n- if ( x_is_match && ncol(x_dataMatrix) > 0 && length(unique(x_predictor))> 1 && x_crossval_i < nrow(x_dataMatrix) ) {\n+ if ( x_is_match\n+ && ncol(x_dataMatrix) > 0\n+ && length(unique(x_predictor))> 1\n+ && x_crossval_i < nrow(x_dataMatrix)\n+ ) {\n my_oplsda <- opls(\n x = x_dataMatrix\n , y = x_predictor\n@@ -19,21 +33,47 @@\n , printL = FALSE\n , plotL = FALSE\n , crossvalI = x_crossval_i\n+ , scaleC = "none" # data have been pareto scaled outside this routine and must not be scaled again. This line fixes issue #2.\n )\n my_oplsda_suppLs_y_levels <- levels(as.factor(my_oplsda@suppLs$y))\n fctr_lvl_1 <- my_oplsda_suppLs_y_levels[1]\n fctr_lvl_2 <- my_oplsda_suppLs_y_levels[2]\n- do_s_plot <- function(parallel_x) {\n- my_cor_vs_cov <- cor_vs_cov(\n- matrix_x = x_dataMatrix\n- , ropls_x = my_oplsda\n- , parallel_x = parallel_x\n- )\n+ do_s_plot <- function(\n+ x_env\n+ , predictor_projection_x = TRUE\n+ , cplot_x = FALSE\n+ , cor_vs_cov_x = NULL\n+ )\n+ {\n+ #print(ls(x_env)) # "cplot_y" etc\n+ #print(str(x_env$cplot_y)) # chr "covariance"\n+ if (cplot_x) {\n+ #print(x_env$cplot_y) # "covariance"\n+ cplot_y_correlation <- (x_env$cplot_y == "correlation")\n+ #print(cplot_y_correlation) # FALSE\n+ }\n+ if (is.null(cor_vs_cov_x)) {\n+ my_cor_vs_cov <- cor_vs_cov(\n+ matrix_x = x_dataMatrix\n+ , ropls_x = my_oplsda\n+ , predictor_projection_x = predictor_projection_x\n+ )\n+ } else {\n+ my_cor_vs_cov <- cor_vs_cov_x\n+ }\n+ # print("str(my_cor_vs_cov)")\n+ # str(my_cor_vs_cov)\n+ if (sum(!is.na(my_cor_vs_cov$tsv1$covariance)) < 2) {\n+ x_progress("No cor_vs_cov data produced")\n+ plot(x=1, y=1, xaxt="n", yaxt="n", xlab="", ylab="", type="n")\n+ text(x=1, y=1, labels="too few covariance data")\n+ return(my_cor_vs_cov)\n+ }\n with(\n my_cor_vs_cov\n , {\n- min_x <- min(covariance)\n- max_x <- max(covariance)\n+ min_x <- min(covariance, na.rm = TRUE)\n+ max_x <- max(covariance, na.rm = TRUE)\n lim_x <- max(sapply(X=c(min_x, max_x), FUN=abs))\n covariance <- covariance / lim_x\n lim_x <- 1.2\n@@ -42,37 +82,78 @@\n # (Mehmood 2012 doi:10.1186/1748-7188-6-27)\n plus_cor <- correlation\n plus_cov <- covariance\n- cex <- 0.75\n+ cex <- 0.65\n which_projection <- if (projection == 1) "t1" else "o1"\n which_loading <- if (projection == 1) "parallel" else "orthogonal"\n- if (projection == 1) {\n- my_xlab <- "relative covariance(feature,t1)"\n- my_x <- plus_cov\n- my_ylab <- "correlation(feature,t1) [~ parallel loading]"\n- my_y <- plus_cor\n- my_xlim <- c( -lim_x - off(0.2), lim_x + off(0.2) )\n+ if (projection == 1) { # predictor-projection\n+ vipcp <- pmax(0, pmin(1,(vip4p-0.83)/(1.21-0.83)))\n+ if (!cplot_x) { # S-plot predictor-projection\n+ my_xlab <- "relative covariance(feature,t1)"\n+ my_x <- plus_cov\n+ my_ylab <- "correlation(feature,t1)"\n+ my_y <- plus_cor\n+ } else { # C-plot predictor-projection\n+ my_xlab <- "variable importance in predictor-projection"\n+ my_x <- vip4p\n+ if (cplot_y_correl'..b'd ropls_x to be of class ropls::opls but instead it was of class ", as.character(x_class) )\n+ if ( !( as.character(x_class) == "opls" ) ) { # || !( attr(class(x_class),"package") == "ropls" ) )\n+ stop(\n+ paste(\n+ "cor_vs_cov: Expected ropls_x to be of class ropls::opls but instead it was of class "\n+ , as.character(x_class)\n+ )\n+ )\n }\n result <- list()\n- result$projection <- projection <- if (parallel_x) 1 else 2\n+ result$projection <- projection <- if (predictor_projection_x) 1 else 2\n # suppLs$algoC - Character: algorithm used - "svd" for singular value decomposition; "nipals" for NIPALS\n if ( ropls_x@suppLs$algoC == "nipals") {\n # Equations (1) and (2) from *Supplement to* Wiklund 2008, doi:10.1021/ac0713510\n mag <- function(one_dimensional) sqrt(sum(one_dimensional * one_dimensional))\n mag_xi <- sapply(X = 1:ncol(matrix_x), FUN = function(x) mag(matrix_x[,x]))\n- if (parallel_x)\n+ if (predictor_projection_x)\n score_matrix <- ropls_x@scoreMN\n else\n score_matrix <- ropls_x@orthoScoreMN\n score_matrix_transposed <- t(score_matrix)\n score_matrix_magnitude <- mag(score_matrix)\n- result$covariance <- score_matrix_transposed %*% matrix_x / ( score_matrix_magnitude * score_matrix_magnitude )\n- result$correlation <- score_matrix_transposed %*% matrix_x / ( score_matrix_magnitude * mag_xi )\n+ result$covariance <-\n+ score_matrix_transposed %*% matrix_x / ( score_matrix_magnitude * score_matrix_magnitude )\n+ result$correlation <-\n+ score_matrix_transposed %*% matrix_x / ( score_matrix_magnitude * mag_xi )\n } else {\n # WARNING - untested code - I don\'t have test data to exercise this branch\n # Equations (1) and (2) from Wiklund 2008, doi:10.1021/ac0713510\n # scoreMN - Numerical matrix of x scores (T; dimensions: nrow(x) x predI) X = TP\' + E; Y = TC\' + F\n- if (parallel_x)\n+ if (predictor_projection_x)\n score_matrix <- ropls_x@scoreMN\n else\n score_matrix <- ropls_x@orthoScoreMN\n@@ -612,9 +863,20 @@\n result$level2 <- rep.int(x = fctr_lvl_2, times = feature_count)\n superresult <- list()\n if (length(result$vip4o) == 0) result$vip4o <- NA\n- greaterLevel <- sapply( X = result$correlation, FUN = function(my_corr) if ( my_corr < 0 ) fctr_lvl_1 else fctr_lvl_2 )\n- superresult$tsv1 <- data.frame(\n- featureID = names(ropls_x@vipVn)\n+ greaterLevel <- sapply(\n+ X = result$correlation\n+ , FUN = function(my_corr) if ( my_corr < 0 ) fctr_lvl_1 else fctr_lvl_2\n+ )\n+\n+ # begin fixes for https://github.com/HegemanLab/w4mcorcov_galaxy_wrapper/issues/1\n+ featureID <- names(ropls_x@vipVn)\n+ greaterLevel <- greaterLevel[featureID]\n+ result$correlation <- result$correlation[featureID]\n+ result$covariance <- result$covariance[featureID]\n+ # end fixes for https://github.com/HegemanLab/w4mcorcov_galaxy_wrapper/issues/1\n+\n+ tsv1 <- data.frame(\n+ featureID = featureID\n , factorLevel1 = result$level1\n , factorLevel2 = result$level2\n , greaterLevel = greaterLevel\n@@ -627,7 +889,10 @@\n , loado = result$loado\n , row.names = NULL\n )\n- rownames(superresult$tsv1) <- superresult$tsv1$featureID\n+ tsv1 <- tsv1[!is.na(tsv1$correlation),]\n+ tsv1 <- tsv1[!is.na(tsv1$covariance),]\n+ superresult$tsv1 <- tsv1\n+ rownames(superresult$tsv1) <- tsv1$featureID\n superresult$projection <- result$projection\n superresult$covariance <- result$covariance\n superresult$correlation <- result$correlation\n@@ -638,7 +903,7 @@\n superresult$details <- result\n result$superresult <- superresult\n # Include thise in case future consumers of this routine want to use it in currently unanticipated ways\n- result$oplsda <- ropls_x \n+ result$oplsda <- ropls_x\n result$predictor <- ropls_x@suppLs$y # in case future consumers of this routine want to use it in currently unanticipated ways\n return (superresult)\n }\n' |
| b |
| diff -r 50f60f94c034 -r 7bd523ca1f9a w4mcorcov_wrapper.R --- a/w4mcorcov_wrapper.R Fri Mar 30 14:59:19 2018 -0400 +++ b/w4mcorcov_wrapper.R Wed Jul 18 12:35:55 2018 -0400 |
| [ |
| @@ -33,7 +33,7 @@ ##--------- my_log <- function(x, ...) { cat(paste(iso8601.znow(), " ", x, ..., nl, sep=""))} -my_fatal <- function(x, ...) { +my_fatal <- function(x, ...) { my_log("ERROR: ", x, ...) quit(save = "no", status = 11, runLast = TRUE) } @@ -72,6 +72,9 @@ my_env$levCSV <- as.character(argVc["levCSV"]) my_env$matchingC <- as.character(argVc["matchingC"]) my_env$labelFeatures <- as.character(argVc["labelFeatures"]) # number of features to label at each extreme of the loadings or 'ALL' +my_env$cplot_o <- as.logical(argVc["cplot_o"]) # TRUE if orthogonal C-plot is requested +my_env$cplot_p <- as.logical(argVc["cplot_p"]) # TRUE if parallel C-plot is requested +my_env$cplot_y <- as.character(argVc["cplot_y"]) # Choice of covariance/correlation for Y-axis on C-plot label_features <- my_env$labelFeatures labelfeatures_check <- TRUE @@ -128,24 +131,70 @@ # compute and plot the correlation_vs_covariance details plot # The parameter settings here are generally taken from bioconductor ropls::plot.opls source. - marVn <- c(4.6, 4.1, 2.6, 1.6) - old_par <- par( - font = 2 # bold font face - , font.axis = 2 # bold font face for axis - , font.lab = 2 # bold font face for x and y labels - , lwd = 2 # line-width - interpretation is device spcific - , mar = marVn # margins - , pch = 18 # black diamond plot-character, see help for graphics::points - # , mfrow = c(2,2) # two rows by two columns - , pty = "s" # force plots to be square - ) + if ( my_env$cplot_p || my_env$cplot_o ) { + old_par <- par( + font = 2 # bold font face + , font.axis = 2 # bold font face for axis + , font.lab = 2 # bold font face for x and y labels + , lwd = 2 # line-width - interpretation is device spcific + , pch = 18 # black diamond plot-character, see help for graphics::points + , pty = "m" # do not force plots to be square + , no.readonly = TRUE # only save writable parameters + ) + pdf_height <- 12 + pdf_width <- 8 + my_layout <- function() { + # lay out 2 columns by 3 rows with extra width at the margin of individual plots + layout( + matrix( + # blank row plot 1 & 2 blank row plot 3 & 4 blank row plot 5 & 6 blank row + c(0,0,0,0,0, 0,1,0,2,0, 0,0,0,0,0, 0,3,0,4,0, 0,0,0,0,0, 0,5,0,6,0, 0,0,0,0,0) + , nrow = 7 + , ncol = 5 + , byrow = TRUE + ) + # slim columns 1, 3, and 5 + , widths = c(0.1, 0.9, 0.1, 0.9, 0.1) + # slim rows 1, 3, 5, and 7 + , heights = c(0.1, 0.9, 0.1, 0.9, 0.1, 0.9, 0.1) + ) + } + } else { + old_par <- par( + font = 2 # bold font face + , font.axis = 2 # bold font face for axis + , font.lab = 2 # bold font face for x and y labels + , lwd = 2 # line-width - interpretation is device spcific + , pch = 18 # black diamond plot-character, see help for graphics::points + , pty = "m" # do not force plots to be square + , no.readonly = TRUE # only save writable parameters + ) + pdf_height <- 8 + pdf_width <- 8 + my_layout <- function() { + # lay out 2 columns by 2 rows with extra width at the margin of individual plots + layout( + matrix( + # blank row plot 1 & 2 blank row plot 3 & 4 blank row + c(0,0,0,0,0, 0,1,0,2,0, 0,0,0,0,0, 0,3,0,4,0, 0,0,0,0,0) + , nrow = 5 + , ncol = 5 + , byrow = TRUE + ) + # slim columns 1, 3, and 5 + , widths = c(0.1, 0.9, 0.1, 0.9, 0.1) + # slim rows 1, 3, and 5 + , heights = c(0.1, 0.9, 0.1, 0.9, 0.1) + ) + } + } plot2pdf( file.name = my_env$contrast_detail - , width = 8 - , height = 8 + , width = pdf_width + , height = pdf_height , plot.function = function() { - # plot layout four plots per page - layout(matrix(1:4, byrow = TRUE, nrow = 2)) + # plot layout four or six plots per page + my_layout() my_result <<- corcov_calc( calc_env = my_env , failure_action = my_fatal @@ -156,7 +205,7 @@ } ) par(old_par) - + my_log( "-------------------------- Finished data processing --------------------------") } |