Mercurial > repos > george-weingart > maaslin

comparison src/lib/Utility.R @ 8:e9677425c6c3 default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Updated the structure of the libraries
author george.weingart@gmail.com
date Mon, 09 Feb 2015 12:17:40 -0500
parents e0b5980139d9
children
comparison
equal deleted inserted replaced
7:c72e14eabb08 8:e9677425c6c3
1 #####################################################################################
2 #Copyright (C) <2012>
3 #
4 #Permission is hereby granted, free of charge, to any person obtaining a copy of
5 #this software and associated documentation files (the "Software"), to deal in the
6 #Software without restriction, including without limitation the rights to use, copy,
7 #modify, merge, publish, distribute, sublicense, and/or sell copies of the Software,
8 #and to permit persons to whom the Software is furnished to do so, subject to
9 #the following conditions:
10 #
11 #The above copyright notice and this permission notice shall be included in all copies
12 #or substantial portions of the Software.
13 #
14 #THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
15 #INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
16 #PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
17 #HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
18 #OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
19 #SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
20 #
21 # This file is a component of the MaAsLin (Multivariate Associations Using Linear Models),
22 # authored by the Huttenhower lab at the Harvard School of Public Health
23 # (contact Timothy Tickle, ttickle@hsph.harvard.edu).
24 #####################################################################################
25
26 inlinedocs <- function(
27 ##author<< Curtis Huttenhower <chuttenh@hsph.harvard.edu> and Timothy Tickle <ttickle@hsph.harvard.edu>
28 ##description<< Collection of minor utility scripts
29 ) { return( pArgs ) }
30
31 #source("Constants.R")
32
33 funcRename <- function(
34 ### Modifies labels for plotting
35 ### If the name is not an otu collapse to the last two clades
36 ### Otherwise use the most terminal clade
37 astrNames
38 ### Names to modify for plotting
39 ){
40 astrRet <- c()
41 for( strName in astrNames )
42 {
43 astrName <- strsplit( strName, c_cFeatureDelimRex )[[1]]
44 i <- length( astrName )
45 if( ( astrName[i] == c_strUnclassified ) || !is.na( as.numeric( astrName[i] ) ) )
46 {
47 strRet <- paste( astrName[( i - 1 ):i], collapse = c_cFeatureDelim )
48 } else {
49 strRet <- astrName[i]
50 }
51 astrRet <- c(astrRet, strRet)
52 }
53 return( astrRet )
54 ### List of modified names
55 }
56
57 funcBonferonniCorrectFactorData <- function
58 ### Bonferroni correct for factor data
59 (dPvalue,
60 ### P-value to correct
61 vsFactors,
62 ### Factors of the data to correct
63 fIgnoreNAs = TRUE
64 ){
65 vsUniqueFactors = unique( vsFactors )
66 if( fIgnoreNAs ){ vsUniqueFactors = setdiff( vsUniqueFactors, c("NA","na","Na","nA") ) }
67 return( dPvalue * max( 1, ( length( vsUniqueFactors ) - 1 ) ) )
68 ### Numeric p-value that is correct for levels (excluding NA levels)
69 }
70
71 funcCalculateTestCounts <- function(
72 ### Calculates the number of tests used in inference
73 iDataCount,
74 asMetadata,
75 asForced,
76 asRandom,
77 fAllvAll
78 ){
79 iMetadata = length(asMetadata)
80 iForced = length(setdiff(intersect( asForced, asMetadata ), asRandom))
81 iRandom = length(intersect( asRandom, asMetadata ))
82 if(fAllvAll)
83 {
84 #AllvAll flow formula
85 return((iMetadata-iForced-iRandom) * iDataCount)
86 }
87
88 #Normal flow formula
89 return((iMetadata-iRandom) * iDataCount)
90 }
91
92 funcGetRandomColors=function(
93 #Generates a given number of random colors
94 tempNumberColors = 1
95 ### Number of colors to generate
96 ){
97 adRet = c()
98 return(sapply(1:tempNumberColors, function(x){
99 adRGB <- ( runif( 3 ) * 0.66 ) + 0.33
100 adRet <- c(adRet, rgb( adRGB[1], adRGB[2], adRGB[3] ))
101 }))
102 }
103
104 funcCoef2Col <- function(
105 ### Searches through a dataframe and looks for a column that would match the coefficient
106 ### by the name of the column or the column name and level appended together.
107 strCoef,
108 ### String coefficient name
109 frmeData,
110 ### Data frame of data
111 astrCols = c()
112 ### Column names of interest (if NULL is given, all column names are inspected).
113 ){
114 #If the coefficient is the intercept there is no data column to return so return null
115 if( strCoef %in% c("(Intercept)", "Intercept") ) { return( NULL ) }
116 #Remove ` from coefficient
117 strCoef <- gsub( "`", "", strCoef )
118
119 #If the coefficient name is not in the data frame
120 if( !( strCoef %in% colnames( frmeData ) ) )
121 {
122 fHit <- FALSE
123 #If the column names are not provided, use the column names of the dataframe.
124 if( is.null( astrCols ) ){astrCols <- colnames( frmeData )}
125
126 #Search through the different column names (factors)
127 for( strFactor in astrCols )
128 {
129 #Select a column, if it is not a factor or does not begin with the factor's name then skip
130 adCur <- frmeData[,strFactor]
131 if( ( class( adCur ) != "factor" ) ||
132 ( substr( strCoef, 1, nchar( strFactor ) ) != strFactor ) ) { next }
133
134 #For the factors, create factor-level name combinations to read in factors
135 #Then check to see the factor-level combination is the coefficient of interest
136 #If it is then store that factor as the coefficient of interest
137 #And break
138 for( strValue in levels( adCur ) )
139 {
140 strCur <- paste( strFactor, strValue, sep = c_sFactorNameSep )
141 if( strCur == strCoef )
142 {
143 strCoef <- strFactor
144 fHit <- TRUE
145 break
146 }
147 }
148
149 #If the factor was found, return
150 if( fHit ){break }
151 }
152 }
153
154 #If the original coefficient or the coefficient factor combination name are in the
155 #data frame, return the name. Otherwise return NA.
156 return( ifelse( ( strCoef %in% colnames( frmeData ) ), strCoef, NA ) )
157 ### Coefficient name
158 }
159
160 funcColToMFAValue = function(
161 ### Given a column name, return the MFA values that could be associated with the name
162 lsColNames,
163 ### String list of column names (as you would get from names(dataframe))
164 dfData
165 ### Data frame of data the column names refer to
166 ){
167 lsMFAValues = c()
168
169 for(sColName in lsColNames)
170 {
171 axCur = dfData[[sColName]]
172
173 if(is.logical(axCur)){axCur=as.factor(axCur)}
174 if(is.factor(axCur))
175 {
176 lsLevels = levels(axCur)
177 if((length(lsLevels)==2) && (!is.na(as.numeric(lsLevels[1]))) && (!is.na(as.numeric(lsLevels[2]))))
178 {
179 lsMFAValues = c(lsMFAValues,paste(sColName,lsLevels[1],sep=c_sMFANameSep1),paste(sColName,lsLevels[2],sep=c_sMFANameSep1))
180 }else{
181 for(sLevel in levels(axCur))
182 {
183 lsMFAValues = c(lsMFAValues,sLevel)
184 }
185 }
186 } else {
187 lsMFAValues = c(lsMFAValues,sColName)
188 }
189 }
190 return(setdiff(lsMFAValues,c("NA",NA)))
191 }
192
193 funcMFAValue2Col = function(
194 ### Given a value in a column, the column name is returned.
195 xValue,
196 dfData,
197 aiColumnIndicesToSearch = NULL
198 ){
199 lsColumnNames = names(dfData)
200
201 if(is.null(aiColumnIndicesToSearch))
202 {
203 aiColumnIndicesToSearch = c(1:dim(dfData)[2])
204 }
205
206 # Could be the column name
207 if(xValue %in% lsColumnNames){return(xValue)}
208
209 # Could be the column name and value
210 iValueLength = length(xValue)
211 for( iColIndex in c(1:length(lsColumnNames) ))
212 {
213 adCur = dfData[[lsColumnNames[iColIndex]]]
214 if(is.factor(adCur))
215 {
216 for(strValue in levels(adCur))
217 {
218 strCurVersion1 <- paste( lsColumnNames[iColIndex], strValue, sep = c_sMFANameSep1 )
219 strCurVersion2 <- paste( lsColumnNames[iColIndex], strValue, sep = c_sMFANameSep2 )
220 if((xValue == strCurVersion1) || (xValue == strCurVersion2)){return(lsColumnNames[iColIndex])}
221 }
222 }
223 }
224
225 # Could be the value
226 for(iColIndex in aiColumnIndicesToSearch)
227 {
228 if(xValue %in% dfData[[lsColumnNames[iColIndex]]]){return(lsColumnNames[iColIndex])}
229 }
230 return(NULL)
231 }
232
233 funcColorHelper <- function(
234 ### Makes sure the max is max and the min is min, and dmed is average
235 dMax = 1,
236 ### Max number
237 dMin = -1,
238 ### Min number
239 dMed = NA
240 ### Average value
241 ){
242 #Make sure max is max and min is min
243 vSort = sort(c(dMin,dMax))
244 return( list( dMin = vSort[1], dMax = vSort[2], dMed = ifelse((is.na(dMed)), (dMin+dMax)/2.0, dMed ) ))
245 ### List of min, max and med numbers
246 }
247
248 funcColor <- function(
249 ### Generate a color based on a number that is forced to be between a min and max range.
250 ### The color is based on how far the number is from the center of the given range
251 ### From red to green (high) are produced with default settings
252 dX,
253 ### Number from which to generate the color
254 dMax = 1,
255 ### Max possible value
256 dMin = -1,
257 ### Min possible value
258 dMed = NA,
259 ### Central value if you don't want to be the average
260 adMax = c(1, 1, 0),
261 ### Is used to generate the color for the higher values in the range, this can be changed to give different colors set to green
262 adMin = c(0, 0, 1),
263 ### Is used to generate the color for the lower values in the range, this can be changed to give different colors set to red
264 adMed = c(0, 0, 0)
265 ### Is used to generate the color for the central values in the range, this can be changed to give different colors set to black
266 ){
267 lsTmp <- funcColorHelper( dMax, dMin, dMed )
268 dMax <- lsTmp$dMax
269 dMin <- lsTmp$dMin
270 dMed <- lsTmp$dMed
271 if( is.na( dX ) )
272 {
273 dX <- dMed
274 }
275 if( dX > dMax )
276 {
277 dX <- dMax
278 } else if( dX < dMin )
279 {
280 dX <- dMin }
281 if( dX < dMed )
282 {
283 d <- ( dMed - dX ) / ( dMed - dMin )
284 adCur <- ( adMed * ( 1 - d ) ) + ( adMin * d )
285 } else {
286 d <- ( dMax - dX ) / ( dMax - dMed )
287 adCur <- ( adMed * d ) + ( adMax * ( 1 - d ) )
288 }
289 return( rgb( adCur[1], adCur[2], adCur[3] ) )
290 ### RGB object
291 }
292
293 funcColors <- function(
294 ### Generate a range of colors
295 dMax = 1,
296 ### Max possible value
297 dMin = -1,
298 ### Min possible value
299 dMed = NA,
300 ### Central value if you don't want to be the average
301 adMax = c(1, 1, 0),
302 ### Is used to generate the color for the higher values in the range, this can be changed to give different colors set to green
303 adMin = c(0, 0, 1),
304 ### Is used to generate the color for the lower values in the range, this can be changed to give different colors set to red
305 adMed = c(0, 0, 0),
306 ### Is used to generate the color for the central values in the range, this can be changed to give different colors set to black
307 iSteps = 64
308 ### Number of intermediary colors made in the range of colors
309 ){
310 lsTmp <- funcColorHelper( dMax, dMin, dMed )
311 dMax <- lsTmp$dMax
312 dMin <- lsTmp$dMin
313 dMed <- lsTmp$dMed
314 aRet <- c ()
315 for( dCur in seq( dMin, dMax, ( dMax - dMin ) / ( iSteps - 1 ) ) )
316 {
317 aRet <- c(aRet, funcColor( dCur, dMax, dMin, dMed, adMax, adMin, adMed ))
318 }
319 return( aRet )
320 ### List of colors
321 }
322
323 funcGetColor <- function(
324 ### Get a color based on col parameter
325 ) {
326 adCol <- col2rgb( par( "col" ) )
327 return( sprintf( "#%02X%02X%02X", adCol[1], adCol[2], adCol[3] ) )
328 ### Return hexadecimal color
329 }
330
331 funcTrim=function(
332 ### Remove whitespace at the beginning or the end of a string
333 tempString
334 ### tempString String to be trimmed.
335 ){
336 return(gsub("^\\s+|\\s+$","",tempString))
337 }
338
339 funcWrite <- function(
340 ### Write a string or a table of data
341 ### This transposes a table before it is written
342 pOut,
343 ### String or table to write
344 strFile
345 ### File to which to write
346 ){
347 if(!is.na(strFile))
348 {
349 if( length( intersect( class( pOut ), c("character", "numeric") ) ) )
350 {
351 write.table( t(pOut), strFile, quote = FALSE, sep = c_cTableDelimiter, col.names = FALSE, row.names = FALSE, na = "", append = TRUE )
352 } else {
353 capture.output( print( pOut ), file = strFile, append = TRUE )
354 }
355 }
356 }
357
358 funcWriteTable <- function(
359 ### Log a table to a file
360 frmeTable,
361 ### Table to write
362 strFile,
363 ### File to which to write
364 fAppend = FALSE
365 ### Append when writing
366 ){
367 if(!is.na(strFile))
368 {
369 write.table( frmeTable, strFile, quote = FALSE, sep = c_cTableDelimiter, na = "", col.names = NA, append = fAppend )
370 }
371 }
372
373 funcWriteQCReport <- function(
374 ### Write out the quality control report
375 strProcessFileName,
376 ### File name
377 lsQCData,
378 ### List of QC data generated by maaslin to be written
379 liDataDim,
380 ### Dimensions of the data matrix
381 liMetadataDim
382 ### Dimensions of the metadata matrix
383 ){
384 unlink(strProcessFileName)
385 funcWrite( paste("Initial Metadata Matrix Size: Rows ",liMetadataDim[1]," Columns ",liMetadataDim[2],sep=""), strProcessFileName )
386 funcWrite( paste("Initial Data Matrix Size: Rows ",liDataDim[1]," Columns ",liDataDim[2],sep=""), strProcessFileName )
387 funcWrite( paste("\nInitial Data Count: ",length(lsQCData$aiDataInitial),sep=""), strProcessFileName )
388 funcWrite( paste("Initial Metadata Count: ",length(lsQCData$aiMetadataInitial),sep=""), strProcessFileName )
389 funcWrite( paste("Data Count after preprocess: ",length(lsQCData$aiAfterPreprocess),sep=""), strProcessFileName )
390 funcWrite( paste("Removed for missing metadata: ",length(lsQCData$iMissingMetadata),sep=""), strProcessFileName )
391 funcWrite( paste("Removed for missing data: ",length(lsQCData$iMissingData),sep=""), strProcessFileName )
392 funcWrite( paste("Number of data with outliers: ",length(which(lsQCData$aiDataSumOutlierPerDatum>0)),sep=""), strProcessFileName )
393 funcWrite( paste("Number of metadata with outliers: ",length(which(lsQCData$aiMetadataSumOutlierPerDatum>0)),sep=""), strProcessFileName )
394 funcWrite( paste("Metadata count which survived clean: ",length(lsQCData$aiMetadataCleaned),sep=""), strProcessFileName )
395 funcWrite( paste("Data count which survived clean: ",length(lsQCData$aiDataCleaned),sep=""), strProcessFileName )
396 funcWrite( paste("\nBoostings: ",lsQCData$iBoosts,sep=""), strProcessFileName )
397 funcWrite( paste("Boosting Errors: ",lsQCData$iBoostErrors,sep=""), strProcessFileName )
398 funcWrite( paste("LMs with no terms suriving boosting: ",lsQCData$iNoTerms,sep=""), strProcessFileName )
399 funcWrite( paste("LMs performed: ",lsQCData$iLms,sep=""), strProcessFileName )
400 if(!is.null(lsQCData$lsQCCustom))
401 {
402 funcWrite("Custom preprocess QC data: ", strProcessFileName )
403 funcWrite(lsQCData$lsQCCustom, strProcessFileName )
404 } else {
405 funcWrite("No custom preprocess QC data.", strProcessFileName )
406 }
407 funcWrite( "\n#Details###########################", strProcessFileName )
408 funcWrite("\nInitial Data Count: ", strProcessFileName )
409 funcWrite(lsQCData$aiDataInitial, strProcessFileName )
410 funcWrite("\nInitial Metadata Count: ", strProcessFileName )
411 funcWrite(lsQCData$aiMetadataInitial, strProcessFileName )
412 funcWrite("\nData Count after preprocess: ", strProcessFileName )
413 funcWrite(lsQCData$aiAfterPreprocess, strProcessFileName )
414 funcWrite("\nRemoved for missing metadata: ", strProcessFileName )
415 funcWrite(lsQCData$iMissingMetadata, strProcessFileName )
416 funcWrite("\nRemoved for missing data: ", strProcessFileName )
417 funcWrite(lsQCData$iMissingData, strProcessFileName )
418 funcWrite("\nDetailed outlier indices: ", strProcessFileName )
419 for(sFeature in names(lsQCData$liOutliers))
420 {
421 funcWrite(paste("Feature",sFeature,"Outlier indice(s):", paste(lsQCData$liOutliers[[sFeature]],collapse=",")), strProcessFileName )
422 }
423 funcWrite("\nMetadata which survived clean: ", strProcessFileName )
424 funcWrite(lsQCData$aiMetadataCleaned, strProcessFileName )
425 funcWrite("\nData which survived clean: ", strProcessFileName )
426 funcWrite(lsQCData$aiDataCleaned, strProcessFileName )
427 }
428
429 funcLMToNoNAFormula <-function(
430 lMod,
431 frmeTmp,
432 adCur
433 ){
434 dfCoef = coef(lMod)
435 astrCoefNames = setdiff(names(dfCoef[as.vector(!is.na(dfCoef))==TRUE]),"(Intercept)")
436 astrPredictors = unique(as.vector(sapply(astrCoefNames,funcCoef2Col, frmeData=frmeTmp)))
437 strFormula = paste( "adCur ~", paste( sprintf( "`%s`", astrPredictors ), collapse = " + " ), sep = " " )
438 return(try( lm(as.formula( strFormula ), data=frmeTmp )))
439 }
440
441 funcFormulaStrToList <- function(
442 #Takes a lm or mixed model formula and returns a list of covariate names in the formula
443 strFormula
444 #Formula to extract covariates from
445 ){
446 #Return list
447 lsRetComparisons = c()
448
449 #If you get a null or na just return
450 if(is.null(strFormula)||is.na(strFormula)){return(lsRetComparisons)}
451
452 #Get test comparisons (predictor names from formula string)
453 asComparisons = gsub("`","",setdiff(unlist(strsplit(unlist(strsplit(strFormula,"~"))[2]," ")),c("","+")))
454
455 #Change metadata in formula to univariate comparisons
456 for(sComparison in asComparisons)
457 {
458 #Removed random covariate formating
459 lsParse = unlist(strsplit(sComparison, "[\\(\\|\\)]", perl=FALSE))
460 lsRetComparisons = c(lsRetComparisons,lsParse[length(lsParse)])
461 }
462 return(lsRetComparisons)
463 }
464
465 funcFormulaListToString <- function(
466 # Using covariate and random covariate names, creates a lm or mixed model formula
467 # returns a vector of c(strLM, strMixedModel), one will be NA given the existance of random covariates.
468 # On error c(NA,NA) is given
469 astrTerms,
470 #Fixed covariates or all covariates if using an lm
471 astrRandomCovariates = NULL
472 #Random covariates for a mixed model
473 ){
474 strRetLMFormula = NA
475 strRetMMFormula = NA
476
477 #If no covariates return NA
478 if(is.null(astrTerms)){return(c(strRetLMFormula, strRetMMFormula))}
479
480 #Get fixed covariates
481 astrFixedCovariates = setdiff(astrTerms,astrRandomCovariates)
482
483 #If no fixed coavariates return NA
484 # Can not run a model with no fixed covariate, restriction of lmm
485 if(length(astrFixedCovariates)==0){return(c(strRetLMFormula, strRetMMFormula))}
486
487 # Fixed Covariates
488 strFixedCovariates = paste( sprintf( "`%s`", astrFixedCovariates ), collapse = " + " )
489
490 #If random covariates, set up a formula for mixed models
491 if(length(astrRandomCovariates)>0)
492 {
493 #Format for lmer
494 #strRetFormula <- paste( "adCur ~ ", paste( sprintf( "(1|`%s`))", intersect(astrRandomCovariates, astrTerms)), collapse = " + " ))
495 #Format for glmmpql
496 strRandomCovariates = paste( sprintf( "1|`%s`", setdiff(astrRandomCovariates, astrTerms)), collapse = " + " )
497 strRetMMFormula <- paste( "adCur ~ ", strFixedCovariates, " + ", strRandomCovariates, sep="")
498 } else {
499 #This is either the formula for all covariates in an lm or fixed covariates in the lmm
500 strRetLMFormula <- paste( "adCur ~ ", strFixedCovariates, sep="")
501 }
502 return(c(strRetLMFormula, strRetMMFormula))
503 }