Mercurial > repos > eschen42 > w4mcorcov
changeset 6:7bd523ca1f9a draft
planemo upload for repository https://github.com/HegemanLab/w4mcorcov_galaxy_wrapper/tree/master commit cafda5095a79ce2376325b57337302f95137195d
| author | eschen42 |
|---|---|
| date | Wed, 18 Jul 2018 12:35:55 -0400 |
| parents | 50f60f94c034 |
| children | 066b1f409e9f |
| files | w4mcorcov.xml w4mcorcov_calc.R w4mcorcov_wrapper.R |
| diffstat | 3 files changed, 574 insertions(+), 168 deletions(-) [+] |
line wrap: on
line diff
--- a/w4mcorcov.xml Fri Mar 30 14:59:19 2018 -0400 +++ b/w4mcorcov.xml Wed Jul 18 12:35:55 2018 -0400 @@ -1,7 +1,7 @@ -<tool id="w4mcorcov" name="OPLS-DA_Contrasts" version="0.98.8"> +<tool id="w4mcorcov" name="OPLS-DA_Contrasts" version="0.98.10"> <description>OPLS-DA Contrasts of Univariate Results</description> - + <macros> <xml name="paramPairSigFeatOnly"> <param @@ -13,13 +13,35 @@ label="Retain only pairwise-significant features" help="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." /> </xml> + <xml name="cplots"> + <param name="cplot_y" label="C-plot Y-axis" type="select" help="Choose the Y-axis for C-plots."> + <option value="correlation">Plot VIP versus correlation</option> + <option value="covariance">Plot VIP versus covariance</option> + </param> + <param + name="cplot_p" + type="boolean" + checked="TRUE" + truevalue="TRUE" + falsevalue="FALSE" + label="Produce predictor C-plot" + help="When this option is set to 'Yes', correlation will be plotted against vip4 for predictor loadings." /> + <param + name="cplot_o" + type="boolean" + checked="TRUE" + truevalue="TRUE" + falsevalue="FALSE" + label="Produce orthogonal C-plot" + help="When this option is set to 'Yes', correlation will be plotted against vip4 for orthogonal loadings." /> + </xml> </macros> <requirements> <requirement type="package">r-batch</requirement> <requirement type="package">bioconductor-ropls</requirement> </requirements> - + <stdio> <exit_code range="1:" level="fatal" /> </stdio> @@ -40,6 +62,15 @@ levCSV '$levCSV' matchingC '$matchingC' labelFeatures '$labelFeatures' + #if str( $xplots.expPlot ) == "none": + cplot_p "FALSE" + cplot_o "FALSE" + cplot_y "correlation" + #else if str( $xplots.expPlot ) == "cplot": + cplot_p "$xplots.cplot_p" + cplot_o "$xplots.cplot_o" + cplot_y "$xplots.cplot_y" + #end if contrast_detail '$contrast_detail' contrast_corcov '$contrast_corcov' contrast_salience '$contrast_salience' @@ -55,7 +86,7 @@ <option value="none">none - Display all features from variableMetadata (rather than choosing a subset based on significance in univariate testing)</option> <option value="ttest">ttest - Student's t-test (parametric test, qualitative factor with exactly 2 levels)</option> <option value="anova">anova - Analysis of variance (parametric test, qualitative factor with more than 2 levels)</option> - <option value="wilcoxon">wilcoxon - Wilcoxon rank test (nonparametric test, qualitative factor with exactly 2 levels)</option> + <option value="wilcoxon">wilcoxon - Wilcoxon rank test (nonparametric test, qualitative factor with exactly 2 levels)</option> <option value="kruskal">kruskal - Kruskal-Wallis rank test (nonparametric test, qualitative factor with more than 2 levels)</option> </param> <when value="none" /> @@ -104,6 +135,16 @@ <option value="regex">use regular expressions for matching level-names</option> </param> <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."/> + <conditional name="xplots"> + <param name="expPlot" label="Extra plots to include" type="select" help="Choosing 'none' hides further choices."> + <option value="none">none - Do not include additonal extra plots.</option> + <option value="cplot">cplot - Include C-plots (predictor-loading vs. 'vip4p' and orthogonal-loading versus 'vip4o')</option> + </param> + <when value="none" /> + <when value="cplot"> + <expand macro="cplots" /> + </when> + </conditional> </inputs> <outputs> @@ -137,6 +178,7 @@ </outputs> <tests> + <!-- test #1 --> <test> <param name="dataMatrix_in" value="input_dataMatrix.tsv"/> <param name="sampleMetadata_in" value="input_sampleMetadata.tsv"/> @@ -160,22 +202,22 @@ <has_text text="level1Level2Sig" /> <!-- first matched line --> <has_text text="M349.2383T700" /> - <has_text text="-0.05007" /> - <has_text text="-5.8455" /> - <has_text text="0.0961269" /> - <has_text text="0.1848301" /> + <has_text text="-0.04051509" /> + <has_text text="-0.001964912" /> + <has_text text="0.02106343" /> + <has_text text="0.2446366813" /> <!-- second matched line --> <has_text text="M207.9308T206" /> - <has_text text="-0.2967565" /> - <has_text text="-19.56942" /> - <has_text text="1.6023" /> - <has_text text="1.35368" /> + <has_text text="0.504885262" /> + <has_text text="0.020749097" /> + <has_text text="0.207196379" /> + <has_text text="0.04438632" /> <!-- third matched line --> <has_text text="M211.0607T263" /> - <has_text text="0.47052" /> - <has_text text="15.910087" /> - <has_text text="0.89838" /> - <has_text text="0.125372" /> + <has_text text="0.0680900" /> + <has_text text="0.0020163" /> + <has_text text="0.0201345" /> + <has_text text="0.0690773" /> </assert_contents> </output> <output name="contrast_salience"> @@ -200,6 +242,7 @@ </assert_contents> </output> </test> + <!-- test #2 --> <test> <param name="dataMatrix_in" value="input_dataMatrix.tsv"/> <param name="sampleMetadata_in" value="input_sampleMetadata.tsv"/> @@ -222,21 +265,11 @@ <has_text text="vip4o" /> <has_text text="level1Level2Sig" /> <!-- first matched line --> - <has_text text="M349.2383T700" /> - <has_text text="-0.99601577" /> - <has_text text="-947.55795176" /> - <!-- second matched line --> - <has_text text="M207.9308T206" /> - <has_text text="0.688549" /> - <has_text text="58.22352" /> - <has_text text="1.394687" /> - <has_text text="0.06049885" /> - <!-- third matched line --> - <has_text text="M211.0607T263" /> - <has_text text="-0.572018" /> - <has_text text="-14.57769" /> - <has_text text="0.7780899" /> - <has_text text="0.3678166776" /> + <has_text text="M200.005T296" /> + <has_text text="-0.2803571" /> + <has_text text="-0.0115899" /> + <has_text text="0.1157346" /> + <has_text text="0.0647860" /> </assert_contents> </output> <output name="contrast_salience"> @@ -261,6 +294,7 @@ </assert_contents> </output> </test> + <!-- test #3 --> <test> <param name="dataMatrix_in" value="input_dataMatrix.tsv"/> <param name="sampleMetadata_in" value="input_sampleMetadata.tsv"/> @@ -282,22 +316,16 @@ <has_text text="vip4o" /> <!-- first matched line --> <has_text text="M349.2383T700" /> - <has_text text="-0.64331257" /> - <has_text text="-161.82220" /> - <has_text text="1.862455" /> - <has_text text="0.2105143" /> + <has_text text="-0.4732226" /> + <has_text text="-0.0506172" /> + <has_text text="0.5246766" /> + <has_text text="0.0103341" /> <!-- second matched line --> <has_text text="M207.9308T206" /> - <has_text text="-0.313507" /> - <has_text text="-20.0476" /> - <has_text text="1.6956987" /> - <has_text text="1.19247" /> - <!-- third matched line --> - <has_text text="M211.0607T263" /> - <has_text text="-0.38986114" /> - <has_text text="-23.747718" /> - <has_text text="1.064296856" /> - <has_text text="1.16507455" /> + <has_text text="0.4927151" /> + <has_text text="0.0203715" /> + <has_text text="0.2111623" /> + <has_text text="0.0488654" /> </assert_contents> </output> <output name="contrast_salience"> @@ -322,6 +350,59 @@ </assert_contents> </output> </test> + <!-- test #4 --> + <test> + <param name="dataMatrix_in" value="issue1_input_dataMatrix.tsv"/> + <param name="sampleMetadata_in" value="issue1_input_sampleMetadata.tsv"/> + <param name="variableMetadata_in" value="issue1_input_variableMetadata.tsv"/> + <param name="tesC" value="none"/> + <param name="facC" value="tissue_flowering"/> + <param name="labelFeatures" value="3"/> + <param name="levCSV" value="*"/> + <param name="matchingC" value="wildcard"/> + <output name="contrast_corcov"> + <assert_contents> + <!-- column-labels line --> + <has_text text="featureID" /> + <has_text text="factorLevel1" /> + <has_text text="factorLevel2" /> + <has_text text="correlation" /> + <has_text text="covariance" /> + <has_text text="vip4p" /> + <has_text text="vip4o" /> + <!-- first matched line --> + <has_text text="NM516T251" /> + <has_text text="flower_yes" /> + <has_text text="other" /> + <has_text text="0.3490559" /> + <has_text text="0.0260147" /> + <has_text text="0.4377872" /> + <has_text text="0.5916089" /> + <has_text text="0.0260147" /> + <has_text text="0.0438942" /> + <has_text text="516.0845" /> + <has_text text="250.8762" /> + </assert_contents> + </output> + <output name="contrast_salience"> + <assert_contents> + <!-- column-labels line --> + <has_text text="featureID" /> + <has_text text="salientLevel" /> + <has_text text="salientRCV" /> + <has_text text="salience" /> + <has_text text="mz" /> + <has_text text="rt" /> + <!-- first matched line --> + <has_text text="PM518T369" /> + <has_text text="flower_yes" /> + <has_text text="0.58655260" /> + <has_text text="4.414469" /> + <has_text text="518.1656" /> + <has_text text="368.59817" /> + </assert_contents> + </output> + </test> </tests> <help><![CDATA[ @@ -356,11 +437,11 @@ - A column of variableMetadata would be labelled 'cluster_kruskal_sig' and would have values '1' and '0'; when the samples are grouped by 'cluster', '1' means that there is strong evidence against the hypothesis that there is no difference among the intensities for the feature across all sample-groups. - A column of variableMetadata would be labelled 'cluster_kruskal_k1.k2_sig' and would have values '1' and '0', where '1' means that there is significant evidence against the hypothesis that samples from sampleMetadata whose 'cluster' column contains 'k1' or 'k2' have the same intensity for that feature. -The 'PLS-DA Contrasts' tool produces graphics and data for OPLS-DA contrasts of feature-intensities between significantly different pairs of factor-levels. For each factor-level, the tool performs a contrast with all other factor-levels combined and then separately with each other factor-level. +The 'PLS-DA Contrasts' tool produces graphics and data for OPLS-DA contrasts of feature-intensities between significantly different pairs of factor-levels. For each factor-level, the tool performs a contrast with all other factor-levels combined and then separately with each other factor-level. **Along the left-to-right axis, the plots show the supervised projection of the variation explained by the predictor** (i.e., the factor specified when invoking the tool); **the top-to-bottom axis displays the variation that is orthogonal to the predictor level** (i.e., independent of it). -Although this tool can be used in a purely exploratory manner by supplying the variableMetadata file without the columns added by the W4M 'Univariate' tool, **the preferred workflow is to use univariate testing to exclude features that are not significantly different and use OPLS-DA to visualize the differences identified in univariate testing** (Th]]>é<![CDATA[venot *et al.*, 2015); an appropriate exception would be to visualize contrasts of a specific list of metabolites. +Although this tool can be used in a purely exploratory manner by supplying the variableMetadata file without the columns added by the W4M 'Univariate' tool, **the preferred workflow may be to use univariate testing to exclude features that are not significantly different and then to use OPLS-DA to visualize the differences identified in univariate testing** (Th]]>é<![CDATA[venot *et al.*, 2015); an appropriate exception would be to visualize contrasts of a specific list of metabolites. It must be stressed that there may be no *single* definitive computational approach to select features that are reliable biomarkers, especially from a small number of samples or experiments. A few possible choices are: @@ -369,7 +450,7 @@ - examining "variable importance in projection VIP for OPLS-DA" (Galindo-Prieto *et al.* 2014), and - examining a feature's "selectivity ratio" (Rajalahti *et al.*, 2009). -In this spirit, this tool reports the S-PLOT covariance and correlation (Wiklund *op. cit.*) and VIP metrics, and it introduces an informal "salience" metric to flag features that may merit attention without dimensional reduction; future versions may add selectivity ratio. +In this spirit, this tool reports the S-PLOT covariance and correlation (Wiklund *op. cit.*) and VIP metrics, and it introduces an informal "salience" metric to flag features that may merit attention without dimensional reduction; future versions may add selectivity ratio. For a more systematic approach to biomarker identification, please consider the W4M 'biosigner' tool (Rinuardo *et al.* 2016), which applies three different identification metrics to the selection process. @@ -453,10 +534,12 @@ [OUT] Contrast-detail output PDF | Several plots for each two-projection OPLS-DA analysis: -- (top-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) -- (bottom-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) -- (top-right) OPLS-DA **scores-plot** for the two projections (Th]]>é<![CDATA[venot *et al.*, 2015) -- (bottom-right) OPLS-DA **loadings-plot** for the two projections (*ibid.*) +- (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. +- (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) +- (first row, right) OPLS-DA **scores-plot** for the two projections (Th]]>é<![CDATA[venot *et al.*, 2015) +- (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). +- (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. +- (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. [OUT] Contrast Correlation-Covarinace data TABULAR | A tab-separated values file of metadata for each feature for each contrast in which it was included. @@ -475,7 +558,7 @@ - **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 [OUT] Feature "Salience" data TABULAR - | 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: + | 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: - **featureID** - feature identifier - **salientLevel** - salient level, i.e., for the feature, the class-level having the greatest median intensity @@ -651,6 +734,15 @@ Release notes ------------- +0.98.10 + +- new feature: C-plots of VIP versus correlation or relative covariance. +- bug fix: Handle issue 2 - features now are only pareto-scaled per Wikland *op cit.*. + +0.98.9 + +- bug fix: Handle issue 1 - handle features removed by ropls because variance is less than 2.2e-16. + 0.98.8 - 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
--- a/w4mcorcov_calc.R Fri Mar 30 14:59:19 2018 -0400 +++ b/w4mcorcov_calc.R Wed Jul 18 12:35:55 2018 -0400 @@ -7,9 +7,23 @@ #### 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, x_crossval_i) { +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 +) { off <- function(x) if (x_show_labels == "0") 0 else x - if ( x_is_match && ncol(x_dataMatrix) > 0 && length(unique(x_predictor))> 1 && x_crossval_i < nrow(x_dataMatrix) ) { + if ( x_is_match + && ncol(x_dataMatrix) > 0 + && length(unique(x_predictor))> 1 + && x_crossval_i < nrow(x_dataMatrix) + ) { my_oplsda <- opls( x = x_dataMatrix , y = x_predictor @@ -19,21 +33,47 @@ , printL = FALSE , plotL = FALSE , crossvalI = x_crossval_i + , scaleC = "none" # data have been pareto scaled outside this routine and must not be scaled again. This line fixes issue #2. ) 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] - do_s_plot <- function(parallel_x) { - my_cor_vs_cov <- cor_vs_cov( - matrix_x = x_dataMatrix - , ropls_x = my_oplsda - , parallel_x = parallel_x - ) + do_s_plot <- function( + x_env + , predictor_projection_x = TRUE + , cplot_x = FALSE + , cor_vs_cov_x = NULL + ) + { + #print(ls(x_env)) # "cplot_y" etc + #print(str(x_env$cplot_y)) # chr "covariance" + if (cplot_x) { + #print(x_env$cplot_y) # "covariance" + cplot_y_correlation <- (x_env$cplot_y == "correlation") + #print(cplot_y_correlation) # FALSE + } + if (is.null(cor_vs_cov_x)) { + my_cor_vs_cov <- cor_vs_cov( + matrix_x = x_dataMatrix + , ropls_x = my_oplsda + , predictor_projection_x = predictor_projection_x + ) + } else { + my_cor_vs_cov <- cor_vs_cov_x + } + # print("str(my_cor_vs_cov)") + # str(my_cor_vs_cov) + if (sum(!is.na(my_cor_vs_cov$tsv1$covariance)) < 2) { + x_progress("No cor_vs_cov data produced") + plot(x=1, y=1, xaxt="n", yaxt="n", xlab="", ylab="", type="n") + text(x=1, y=1, labels="too few covariance data") + return(my_cor_vs_cov) + } with( my_cor_vs_cov , { - min_x <- min(covariance) - max_x <- max(covariance) + min_x <- min(covariance, na.rm = TRUE) + max_x <- max(covariance, na.rm = TRUE) lim_x <- max(sapply(X=c(min_x, max_x), FUN=abs)) covariance <- covariance / lim_x lim_x <- 1.2 @@ -42,37 +82,78 @@ # (Mehmood 2012 doi:10.1186/1748-7188-6-27) plus_cor <- correlation plus_cov <- covariance - cex <- 0.75 + cex <- 0.65 which_projection <- if (projection == 1) "t1" else "o1" which_loading <- if (projection == 1) "parallel" else "orthogonal" - if (projection == 1) { - my_xlab <- "relative covariance(feature,t1)" - my_x <- plus_cov - my_ylab <- "correlation(feature,t1) [~ parallel loading]" - my_y <- plus_cor - my_xlim <- c( -lim_x - off(0.2), lim_x + off(0.2) ) + if (projection == 1) { # predictor-projection + vipcp <- pmax(0, pmin(1,(vip4p-0.83)/(1.21-0.83))) + if (!cplot_x) { # S-plot predictor-projection + my_xlab <- "relative covariance(feature,t1)" + my_x <- plus_cov + my_ylab <- "correlation(feature,t1)" + my_y <- plus_cor + } else { # C-plot predictor-projection + my_xlab <- "variable importance in predictor-projection" + my_x <- vip4p + if (cplot_y_correlation) { + my_ylab <- "correlation(feature,t1)" + my_y <- plus_cor + } else { + my_ylab <- "relative covariance(feature,t1)" + my_y <- plus_cov + } + } + if (cplot_x) { + lim_x <- max(my_x, na.rm = TRUE) * 1.1 + my_xlim <- c( 0, lim_x + off(0.2) ) + } else { + my_xlim <- c( -lim_x - off(0.2), lim_x + off(0.2) ) + } my_ylim <- c( -1.0 - off(0.2), 1.0 + off(0.2) ) my_load_distal <- loadp my_load_proximal <- loado - vipcp <- pmax(0, pmin(1,(vip4p-0.83)/(1.21-0.83))) red <- as.numeric(correlation > 0) * vipcp blue <- as.numeric(correlation < 0) * vipcp alpha <- 0.1 + 0.4 * vipcp - my_col = rgb(blue = blue, red = red, green = 0, alpha = alpha) - main_label <- sprintf("%s for level %s versus %s", x_prefix, fctr_lvl_1, fctr_lvl_2) - } else { - my_xlab <- "relative covariance(feature,to1)" - my_x <- -plus_cov - my_xlim <- c( -lim_x - off(0.2), lim_x + off(0.2) ) - my_ylab <- "correlation(feature,to1) [~ orthogonal loading]" - my_y <- plus_cor + red[is.na(red)] <- 0 + blue[is.na(blue)] <- 0 + alpha[is.na(alpha)] <- 0 + my_col <- rgb(blue = blue, red = red, green = 0, alpha = alpha) + main_label <- sprintf("%s for level %s versus %s" + , x_prefix, fctr_lvl_1, fctr_lvl_2) + } else { # orthogonal projection + vipco <- pmax(0, pmin(1,(vip4o-0.83)/(1.21-0.83))) + if (!cplot_x) { + my_xlab <- "relative covariance(feature,to1)" + my_x <- -plus_cov + } else { + my_xlab <- "variable importance in orthogonal projection" + my_x <- vip4o + } + if (!cplot_x) { # S-plot orthogonal projection + my_xlim <- c( -lim_x - off(0.2), lim_x + off(0.2) ) + my_ylab <- "correlation(feature,to1)" + my_y <- plus_cor + } else { # C-plot orthogonal projection + lim_x <- max(my_x, na.rm = TRUE) * 1.1 + my_xlim <- c( 0, lim_x + off(0.2) ) + if (cplot_y_correlation) { + my_ylab <- "correlation(feature,to1)" + my_y <- plus_cor + } else { + my_ylab <- "relative covariance(feature,to1)" + my_y <- plus_cov + } + } my_ylim <- c( -1.0 - off(0.2), 1.0 + off(0.2) ) my_load_distal <- loado my_load_proximal <- loadp - vipco <- pmax(0, pmin(1,(vip4o-0.83)/(1.21-0.83))) alpha <- 0.1 + 0.4 * vipco - my_col = rgb(blue = 0, red = 0, green = 0, alpha = alpha) - main_label <- sprintf("Features influencing orthogonal projection for level %s versus %s", fctr_lvl_1, fctr_lvl_2) + alpha[is.na(alpha)] <- 0 + my_col <- rgb(blue = 0, red = 0, green = 0, alpha = alpha) + main_label <- sprintf( + "Features influencing orthogonal projection for %s versus %s" + , fctr_lvl_1, fctr_lvl_2) } main_cex <- min(1.0, 46.0/nchar(main_label)) my_feature_label_slant <- -30 # slant feature labels 30 degrees downward @@ -92,7 +173,7 @@ ) low_x <- -0.7 * lim_x high_x <- 0.7 * lim_x - if (projection == 1) { + if (projection == 1 && !cplot_x) { text(x = low_x, y = -0.05, labels = fctr_lvl_1, col = "blue") text(x = high_x, y = 0.05, labels = fctr_lvl_2, col = "red") } @@ -109,54 +190,105 @@ names(head(sort(my_load_distal),n = n_labels)) , names(tail(sort(my_load_distal),n = n_labels)) ) - labels <- unname(sapply( X = tsv1$featureID, FUN = function(x) if( x %in% labels_to_show ) x else "" )) + labels <- unname( + sapply( + X = tsv1$featureID + , FUN = function(x) if( x %in% labels_to_show ) x else "" + ) + ) x_text_offset <- 0.024 - y_text_offset <- (if (projection == 1) 1 else -1) * -0.017 + y_text_off <- 0.017 + if (!cplot_x) { # S-plot + y_text_offset <- if (projection == 1) -y_text_off else y_text_off + } else { # C-plot + y_text_offset <- + sapply( + X = (my_y > 0) + , FUN = function(x) { if (x) y_text_off else -y_text_off } + ) + } label_features <- function(x_arg, y_arg, labels_arg, slant_arg) { - print("str(x_arg)") - print(str(x_arg)) - print("str(y_arg)") - print(str(y_arg)) - print("str(labels_arg)") - print(str(labels_arg)) - text( - y = y_arg - , x = x_arg + x_text_offset - , cex = 0.4 - , labels = labels_arg - , col = rgb(blue = 0, red = 0, green = 0, alpha = 0.5) # grey semi-transparent labels - , srt = slant_arg - , adj = 0 # left-justified - ) + # print("str(x_arg)") + # print(str(x_arg)) + # print("str(y_arg)") + # print(str(y_arg)) + # print("str(labels_arg)") + # print(str(labels_arg)) + if (length(labels_arg) > 0) { + unique_slant <- unique(slant_arg) + if (length(unique_slant) == 1) { + text( + y = y_arg + , x = x_arg + x_text_offset + , cex = 0.4 + , labels = labels_arg + , col = rgb(blue = 0, red = 0, green = 0, alpha = 0.5) # grey semi-transparent labels + , srt = slant_arg + , adj = 0 # left-justified + ) + } else { + for (slant in unique_slant) { + text( + y = y_arg[slant_arg == slant] + , x = x_arg[slant_arg == slant] + x_text_offset + , cex = 0.4 + , labels = labels_arg[slant_arg == slant] + , col = rgb(blue = 0, red = 0, green = 0, alpha = 0.5) # grey semi-transparent labels + , srt = slant + , adj = 0 # left-justified + ) + } + } + } } - my_slant <- (if (projection == 1) 1 else -1) * my_feature_label_slant + if (!cplot_x) { + my_slant <- (if (projection == 1) 1 else -1) * my_feature_label_slant + } else { + my_slant <- sapply( + X = (my_y > 0) + , FUN = function(x) if (x) 2 else -2 + ) * my_feature_label_slant + } if (length(my_x) > 1) { - label_features( + label_features( x_arg = my_x [my_x > 0] , y_arg = my_y [my_x > 0] - y_text_offset , labels_arg = labels[my_x > 0] - , slant_arg = -my_slant + , slant_arg = (if (!cplot_x) -my_slant else (my_slant)) ) - label_features( - x_arg = my_x [my_x < 0] - , y_arg = my_y [my_x < 0] + y_text_offset - , labels_arg = labels[my_x < 0] + if (!cplot_x) { + label_features( + x_arg = my_x [my_x < 0] + , y_arg = my_y [my_x < 0] + y_text_offset + , labels_arg = labels[my_x < 0] + , slant_arg = my_slant + ) + } + } else { + if (!cplot_x) { + my_slant <- (if (my_x > 1) -1 else 1) * my_slant + my_y_arg = my_y + (if (my_x > 1) -1 else 1) * y_text_offset + } else { + my_slant <- (if (my_y > 1) -1 else 1) * my_slant + my_y_arg = my_y + (if (my_y > 1) -1 else 1) * y_text_offset + } + label_features( + x_arg = my_x + , y_arg = my_y_arg + , labels_arg = labels , slant_arg = my_slant ) - } else { - label_features( - x_arg = my_x - , y_arg = my_y + (if (my_x > 1) -1 else 1) * y_text_offset - , labels_arg = labels - , slant_arg = (if (my_x > 1) -1 else 1) * my_slant - ) } } } ) return (my_cor_vs_cov) } - my_cor_vs_cov <- do_s_plot( parallel_x = TRUE ) + my_cor_vs_cov <- do_s_plot( + x_env = x_env + , predictor_projection_x = TRUE + , cplot_x = FALSE + ) typeVc <- c("correlation", # 1 "outlier", # 2 "overview", # 3 @@ -175,6 +307,7 @@ } else { my_typevc <- c("(dummy)","overview","(dummy)") } + my_ortho_cor_vs_cov_exists <- FALSE for (my_type in my_typevc) { if (my_type %in% typeVc) { tryCatch({ @@ -187,17 +320,65 @@ , parLayL = TRUE , parEllipsesL = TRUE ) + if (my_type == "overview") { + sub_label <- sprintf("%s versus %s", fctr_lvl_1, fctr_lvl_2) + title(sub = sub_label) + } } else { - do_s_plot( parallel_x = FALSE ) + my_ortho_cor_vs_cov <- do_s_plot( + x_env = x_env + , predictor_projection_x = FALSE + , cplot_x = FALSE + ) + my_ortho_cor_vs_cov_exists <- TRUE } }, error = function(e) { - x_progress( sprintf( "factor level %s or %s may have only one sample - %s", fctr_lvl_1, fctr_lvl_2, e$message ) ) + x_progress( + sprintf( + "factor level %s or %s may have only one sample - %s" + , fctr_lvl_1 + , fctr_lvl_2 + , e$message + ) + ) }) } else { plot(x=1, y=1, xaxt="n", yaxt="n", xlab="", ylab="", type="n") text(x=1, y=1, labels="no orthogonal projection is possible") } } + cplot_p <- x_env$cplot_p + cplot_o <- x_env$cplot_o + if (cplot_p || cplot_o) { + if (cplot_p) { + do_s_plot( + x_env = x_env + , predictor_projection_x = TRUE + , cplot_x = TRUE + , cor_vs_cov_x = my_cor_vs_cov + ) + did_plots <- 1 + } else { + did_plots <- 0 + } + if (cplot_o) { + if (my_ortho_cor_vs_cov_exists) { + do_s_plot( + x_env = x_env + , predictor_projection_x = FALSE + , cplot_x = TRUE + , cor_vs_cov_x = my_ortho_cor_vs_cov + ) + } else { + plot(x=1, y=1, xaxt="n", yaxt="n", xlab="", ylab="", type="n") + text(x=1, y=1, labels="no orthogonal projection is possible") + } + did_plots <- 1 + did_plots + } + if (did_plots == 1) { + plot(x=1, y=1, xaxt="n", yaxt="n", xlab="", ylab="", type="n", fg = "white") + } + } return (my_cor_vs_cov) } else { return (NULL) @@ -205,7 +386,14 @@ } # 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){}) { +corcov_calc <- function( + calc_env +, failure_action = stop +, progress_action = function(x) { } +, corcov_tsv_action = function(t) { } +, salience_tsv_action = function(t) { } +, extra_plots = c() +) { if ( ! is.environment(calc_env) ) { failure_action("corcov_calc: fatal error - 'calc_env' is not an environment") return ( FALSE ) @@ -235,18 +423,20 @@ # arg/env checking if (!(facC %in% names(smpl_metadata))) { - failure_action(sprintf("bad parameter! Factor name '%s' not found in sampleMetadata", facC)) + failure_action( + sprintf("bad parameter! Factor name '%s' not found in sampleMetadata" + , facC)) return ( FALSE ) } mz <- vrbl_metadata$mz names(mz) <- vrbl_metadata$variableMetadata mz_lookup <- function(feature) unname(mz[feature]) - + rt <- vrbl_metadata$rt names(rt) <- vrbl_metadata$variableMetadata rt_lookup <- function(feature) unname(rt[feature]) - + # 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 @@ -322,7 +512,10 @@ # 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)) + failure_action( + sprintf( + "bad parameter! variableMetadata must contain results of W4M Univariate test '%s'." + , tesC)) return ( FALSE ) } col_matches <- lapply( @@ -349,21 +542,36 @@ } } level_union <- unique(sort(level_union)) - overall_significant <- 1 == ( if (intersample_sig_col %in% colnames(vrbl_metadata)) vrbl_metadata[,intersample_sig_col] else TRUE ) + overall_significant <- 1 == ( + if (intersample_sig_col %in% colnames(vrbl_metadata)) { + vrbl_metadata[,intersample_sig_col] + } else { + 1 + } + ) if ( length(level_union) > 2 ) { 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 ] plot_action <- function(fctr_lvl_1, fctr_lvl_2) { - 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 ) + 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_prefix = if (pairSigFeatOnly) { + "Significantly contrasting features" + } else { + "Significant features" + } , x_show_labels = labelFeatures , x_progress = progress_action , x_crossval_i = min(7, length(chosen_samples)) @@ -375,7 +583,10 @@ my_tsv <- my_cor_cov$tsv1 my_tsv$mz <- mz_lookup(my_tsv$featureID) my_tsv$rt <- rt_lookup(my_tsv$featureID) - my_tsv["level1Level2Sig"] <- vrbl_metadata[ match(my_tsv$featureID, vrbl_metadata_names), vrbl_metadata_col ] + my_tsv["level1Level2Sig"] <- vrbl_metadata[ + match(my_tsv$featureID, vrbl_metadata_names) + , vrbl_metadata_col + ] tsv <<- my_tsv corcov_tsv_action(tsv) did_plot <<- TRUE @@ -404,20 +615,34 @@ 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)) + progress_action( + sprintf("calculating/plotting contrast of %s vs. %s" + , fctr_lvl_1, fctr_lvl_2)) # 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] * ( if (intersample_sig_col %in% colnames(vrbl_metadata)) vrbl_metadata[,intersample_sig_col] else TRUE ) - col_selector <- vrbl_metadata_names[ if ( pairSigFeatOnly ) fully_significant else overall_significant ] + fully_significant <- 1 == vrbl_metadata[,vrbl_metadata_col] * + ( if (intersample_sig_col %in% colnames(vrbl_metadata)) { + vrbl_metadata[,intersample_sig_col] + } else { + 1 + } + ) + 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_prefix = if (pairSigFeatOnly) { + "Significantly contrasting features" + } else { + "Significant features" + } , x_show_labels = labelFeatures , x_progress = progress_action , x_crossval_i = min(7, length(chosen_samples)) @@ -429,7 +654,10 @@ tsv <- my_cor_cov$tsv1 tsv$mz <- mz_lookup(tsv$featureID) tsv$rt <- rt_lookup(tsv$featureID) - tsv["level1Level2Sig"] <- vrbl_metadata[ match(tsv$featureID, vrbl_metadata_names), vrbl_metadata_col ] + tsv["level1Level2Sig"] <- vrbl_metadata[ + match(tsv$featureID, vrbl_metadata_names) + , vrbl_metadata_col + ] corcov_tsv_action(tsv) did_plot <- TRUE } @@ -444,7 +672,7 @@ completed <- c() lapply( X = level_union - , FUN = function(x) { + , FUN = function(x) { fctr_lvl_1 <- x[1] fctr_lvl_2 <- { if ( fctr_lvl_1 %in% completed ) @@ -455,12 +683,20 @@ } chosen_samples <- smpl_metadata_facC %in% c(fctr_lvl_1, fctr_lvl_2) fctr_lvl_2 <- "other" - progress_action(sprintf("calculating/plotting contrast of %s vs. %s", fctr_lvl_1, 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 ) + 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) @@ -494,11 +730,13 @@ utils::combn( x = level_union , m = 2 - , FUN = function(x) { + , 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) - progress_action(sprintf("calculating/plotting contrast of %s vs. %s", fctr_lvl_1, 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 { @@ -536,7 +774,10 @@ } } if (!did_plot) { - failure_action(sprintf("bad parameter! sampleMetadata must have at least two levels of factor '%s' matching '%s'", facC, originalLevCSV)) + failure_action( + sprintf( + "bad parameter! sampleMetadata must have at least two levels of factor '%s' matching '%s'" + , facC, originalLevCSV)) return ( FALSE ) } return ( TRUE ) @@ -547,31 +788,41 @@ # 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, parallel_x = TRUE) { +cor_vs_cov <- function( + matrix_x, ropls_x +, predictor_projection_x = TRUE +) { 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) ) + if ( !( as.character(x_class) == "opls" ) ) { # || !( attr(class(x_class),"package") == "ropls" ) ) + stop( + paste( + "cor_vs_cov: Expected ropls_x to be of class ropls::opls but instead it was of class " + , as.character(x_class) + ) + ) } result <- list() - result$projection <- projection <- if (parallel_x) 1 else 2 + result$projection <- projection <- if (predictor_projection_x) 1 else 2 # 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])) - if (parallel_x) + if (predictor_projection_x) score_matrix <- ropls_x@scoreMN else score_matrix <- ropls_x@orthoScoreMN 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 ) + 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 - if (parallel_x) + if (predictor_projection_x) score_matrix <- ropls_x@scoreMN else score_matrix <- ropls_x@orthoScoreMN @@ -612,9 +863,20 @@ result$level2 <- rep.int(x = fctr_lvl_2, times = feature_count) superresult <- list() if (length(result$vip4o) == 0) result$vip4o <- NA - greaterLevel <- sapply( X = result$correlation, FUN = function(my_corr) if ( my_corr < 0 ) fctr_lvl_1 else fctr_lvl_2 ) - superresult$tsv1 <- data.frame( - featureID = names(ropls_x@vipVn) + greaterLevel <- sapply( + X = result$correlation + , FUN = function(my_corr) if ( my_corr < 0 ) fctr_lvl_1 else fctr_lvl_2 + ) + + # begin fixes for https://github.com/HegemanLab/w4mcorcov_galaxy_wrapper/issues/1 + featureID <- names(ropls_x@vipVn) + greaterLevel <- greaterLevel[featureID] + result$correlation <- result$correlation[featureID] + result$covariance <- result$covariance[featureID] + # end fixes for https://github.com/HegemanLab/w4mcorcov_galaxy_wrapper/issues/1 + + tsv1 <- data.frame( + featureID = featureID , factorLevel1 = result$level1 , factorLevel2 = result$level2 , greaterLevel = greaterLevel @@ -627,7 +889,10 @@ , loado = result$loado , row.names = NULL ) - rownames(superresult$tsv1) <- superresult$tsv1$featureID + tsv1 <- tsv1[!is.na(tsv1$correlation),] + tsv1 <- tsv1[!is.na(tsv1$covariance),] + superresult$tsv1 <- tsv1 + rownames(superresult$tsv1) <- tsv1$featureID superresult$projection <- result$projection superresult$covariance <- result$covariance superresult$correlation <- result$correlation @@ -638,7 +903,7 @@ superresult$details <- result 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$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) }
--- 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 --------------------------") }