11# Module script for: biolockj.module.report.r.R_PlotMds
22
33# Import vegan library for distance plot support
4- # Main function generates 3 MDS plots for each attribute at each level in taxaLevels()
4+ # Main function generates numAxis MDS plots for each field at each level in taxaLevels()
55main <- function () {
66 importLibs( c( " vegan" ) )
7- numAxis = getProperty(" r_PlotMds.numAxis" )
87 mdsFields = getProperty( " r_PlotMds.reportFields" , c( getBinaryFields(), getNominalFields() ) )
98
109 for ( level in taxaLevels() ) {
@@ -13,8 +12,10 @@ main <- function() {
1312 metaTable = getMetaData( level )
1413 if ( is.null(countTable ) || is.null(metaTable ) ) { next }
1514 if ( doDebug() ) sink( getLogFile( level ) )
15+ logInfo( " mdsFields" , mdsFields )
1616
1717 myMDS = capscale( countTable ~ 1 , distance = getProperty(" r_PlotMds.distance" ) )
18+ numAxis = min( c( getProperty(" r_PlotMds.numAxis" ), ncol(myMDS $ CA $ u ) ) )
1819 metaColColors = getColorsByCategory( metaTable )
1920
2021 pcoaFileName = paste0( getPath( file.path(getModuleDir(), " temp" ), paste0(level , " _pcoa" ) ), " .tsv" )
@@ -25,57 +26,42 @@ main <- function() {
2526 write.table( data.frame (mds = names(myMDS $ CA $ eig ), eig = myMDS $ CA $ eig ), file = eigenFileName , col.names = FALSE , row.names = FALSE , sep = " \t " )
2627 logInfo( " Save Eigen value table" , pcoaFileName )
2728
28- # Make plots
2929 outputFile = paste0( getPath( getOutputDir(), paste0(level , " _MDS.pdf" ) ) )
3030 pdf( outputFile , paper = " letter" , width = 7.5 , height = 10.5 )
31- par(mfrow = c(3 , 2 ), las = 1 , oma = c(1 ,0 ,2 ,1 ), mar = c(5 , 4 , 2 , 2 ), cex = .95 )
32- percentVariance = as.numeric(eigenvals(myMDS )/ sum( eigenvals(myMDS ) ) ) * 100
31+ par( mfrow = c(3 , 2 ), las = 1 , oma = c(1 ,0 ,2 ,1 ), mar = c(5 , 4 , 2 , 2 ), cex = 0 .95 )
32+ perVariance = as.numeric(eigenvals(myMDS )/ sum( eigenvals(myMDS ) ) ) * 100
3333 pageNum = 0
3434
3535 for ( field in mdsFields ){
36- logInfo( " mdsFields" , mdsFields )
37- pageNum = pageNum + 1
38- metaColVals = as.character(metaTable [,field ])
39- logInfo( " metaColVals" , metaColVals )
40- par(mfrow = par(" mfrow" ), cex = par(" cex" ))
41- att = as.factor(metaColVals )
42- colorKey = metaColColors [[field ]]
43- logInfo( c( " Using colors: " , paste(colorKey , " for" , names(colorKey ), collapse = " , " )) )
44- position = 1
45- pageNum = 1
46- numAxis = min(c(numAxis , ncol(myMDS $ CA $ u )))
47- for (x in 1 : (numAxis - 1 )) {
48- for (y in (x + 1 ): numAxis ) {
49- if (position > prod(par(" mfrow" ) ) ) {
50- position = 1
51- pageNum = pageNum + 1
52- }
53- pch = getProperty(" r.pch" , 20 )
36+ par(mfrow = par(" mfrow" ), cex = par(" cex" ))
37+ position = 1
38+
39+ metaColVals = as.character(metaTable [,field ])
40+ colorKey = metaColColors [[field ]]
41+
42+ logInfo( " metaColVals" , metaColVals )
43+ logInfo( c( " Using colors: " , paste(colorKey , " for" , names(colorKey ), collapse = " , " )) )
44+
45+ for ( x in 1 : (numAxis - 1 ) ) {
46+ for ( y in (x + 1 ): numAxis ) {
47+
5448 plot( myMDS $ CA $ u [,x ], myMDS $ CA $ u [,y ], main = paste(" Axes" , x , " vs" , y ),
55- xlab = getMdsLabel( x , percentVariance [x ] ),
56- ylab = getMdsLabel( y , percentVariance [y ] ),
57- cex = 1.2 , pch = pch , col = colorKey [metaColVals ] )
49+ xlab = getMdsLabel( x , perVariance [x ] ),
50+ ylab = getMdsLabel( y , perVariance [y ] ),
51+ cex = 1.2 , pch = getProperty(" r.pch" , 20 ), col = colorKey [metaColVals ] )
52+
53+
54+ if ( position == 1 || position > prod( par(" mfrow" ) ) ) {
55+ position = 1
56+ pageNum = pageNum + 1
57+ addHeaderFooter( field , level , pageNum )
58+ }
59+
5860 position = position + 1
59- if ( position == 2 ){
60- addPageTitle( field , line = 1 )
61- addPageNumber( pageNum )
62- addPageFooter( " Multidimensional Scaling" )
63- # put this plot at the upper right position
64- # that puts the legend in a nice white space, and it makes axis 1 in line with itself in two plots (same for axis3)
65- plotRelativeVariance(percentVariance , numAxis )
66- position = position + 1
67- title( displayLevel( level ) )
68- # Add legend
69- legendKey = colorKey
70- legendLabels = paste0(names(legendKey ), " (n=" , table(metaColVals )[names(legendKey )], " )" )
71- legendKey = legendKey [ order(table(metaColVals )[names(colorKey )]) ]
72- maxInLegend = 6
73- if (length(colorKey ) > (maxInLegend + 1 )){
74- legendKey = c( colorKey [ 1 : maxInLegend ], NA )
75- numDropped = length(colorKey ) - length(legendKey ) + 1
76- legendLabels = c(legendLabels [1 : maxInLegend ], paste(" (" , numDropped , " other labels )" ))
77- }
78- legend(x = " topright" , title = field , legend = legendLabels , col = legendKey , pch = pch , bty = " n" )
61+
62+ if ( position == 2 ) {
63+ plotRelativeVariance( field , metaColVals , perVariance , level , numAxis , colorKey )
64+ position = position + 1
7965 }
8066 }
8167 }
@@ -85,19 +71,42 @@ main <- function() {
8571 }
8672}
8773
74+ # Add page title + footer with page number
75+ addHeaderFooter <- function ( field , level , pageNum ) {
76+ addPageTitle( field )
77+ addPageNumber( pageNum )
78+ addPageFooter( paste( str_to_title( level ), " Multidimensional Scaling" ) )
79+ }
80+
81+ # Get variance plot label as percentage
8882getMdsLabel <- function ( axisNum , variance ) {
89- return ( paste0(" Axis " , axisNum , " (" , paste0( round( variance ), " %)" ) ) )
83+ return ( paste0(" Axis " , axisNum , " ( " , paste0( round( variance ), " % )" ) ) )
9084}
9185
92- plotRelativeVariance <- function (percentVariance , numAxis ){
93- numBars = min(c(length(percentVariance ), 6 )) # arbitrary choice, don't show more than 6
86+ # This plot is always put in the upper right corner of the page
87+ plotRelativeVariance <- function ( field , metaColVals , perVariance , level , numAxis , colorKey ){
88+ numBars = min( c(length(perVariance ), maxInLegend ) )
9489 numBars = max(numBars , numAxis )
95- heights = percentVariance [1 : numBars ]
90+ heights = perVariance [1 : numBars ]
9691 bp = barplot(heights , col = " dodgerblue1" , ylim = c(0 ,100 ), names = 1 : numBars ,
9792 xlab = " Axis" , ylab = " Variance" )
9893 labels = round(heights )
9994 near0 = which(labels < 1 )
10095 labels [near0 ] = " <1"
101- if ( numBars < = 6 ){ labels = paste(labels , " %" ) }
96+ if ( numBars < = maxInLegend ) labels = paste( labels , " %" )
10297 text(x = bp , y = heights , labels = labels , pos = 3 , xpd = TRUE )
98+ title( str_to_title( level ) )
99+
100+ legendKey = colorKey
101+ legendLabels = paste0(names(legendKey ), " (n=" , table(metaColVals )[names(legendKey )], " )" )
102+ legendKey = legendKey [ order(table(metaColVals )[names(colorKey )]) ]
103+
104+ if (length(colorKey ) > (maxInLegend + 1 ) ){
105+ legendKey = c( colorKey [ 1 : maxInLegend ], NA )
106+ numDropped = length(colorKey ) - length(legendKey ) + 1
107+ legendLabels = c(legendLabels [1 : maxInLegend ], paste(" (" , numDropped , " other labels )" ))
108+ }
109+ legend( " topright" , title = field , legend = legendLabels , cex = 0.8 , col = legendKey , pch = getProperty(" r.pch" , 20 ), bty = " n" )
103110}
111+
112+ maxInLegend = 6
0 commit comments