Mercurial > repos > george-weingart > maaslin
diff maaslin-4450aa4ecc84/src/lib/BoostGLM.R @ 1:a87d5a5f2776
Uploaded the version running on the prod server
| author | george-weingart |
|---|---|
| date | Sun, 08 Feb 2015 23:08:38 -0500 |
| parents | |
| children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/maaslin-4450aa4ecc84/src/lib/BoostGLM.R Sun Feb 08 23:08:38 2015 -0500 @@ -0,0 +1,887 @@ +##################################################################################### +#Copyright (C) <2012> +# +#Permission is hereby granted, free of charge, to any person obtaining a copy of +#this software and associated documentation files (the "Software"), to deal in the +#Software without restriction, including without limitation the rights to use, copy, +#modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, +#and to permit persons to whom the Software is furnished to do so, subject to +#the following conditions: +# +#The above copyright notice and this permission notice shall be included in all copies +#or substantial portions of the Software. +# +#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, +#INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A +#PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +#HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION +#OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE +#SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. +# +# This file is a component of the MaAsLin (Multivariate Associations Using Linear Models), +# authored by the Huttenhower lab at the Harvard School of Public Health +# (contact Timothy Tickle, ttickle@hsph.harvard.edu). +##################################################################################### + +inlinedocs <- function( +##author<< Curtis Huttenhower <chuttenh@hsph.harvard.edu> and Timothy Tickle <ttickle@hsph.harvard.edu> +##description<< Manages the quality control of data and the performance of analysis (univariate or multivariate), regularization, and data (response) transformation. +) { return( pArgs ) } + +### Load libraries quietly +suppressMessages(library( gam, warn.conflicts=FALSE, quietly=TRUE, verbose=FALSE)) +suppressMessages(library( gbm, warn.conflicts=FALSE, quietly=TRUE, verbose=FALSE)) +suppressMessages(library( logging, warn.conflicts=FALSE, quietly=TRUE, verbose=FALSE)) +suppressMessages(library( outliers, warn.conflicts=FALSE, quietly=TRUE, verbose=FALSE)) +suppressMessages(library( robustbase, warn.conflicts=FALSE, quietly=TRUE, verbose=FALSE)) +suppressMessages(library( pscl, warn.conflicts=FALSE, quietly=TRUE, verbose=FALSE)) + +### Get constants +#source(file.path("input","maaslin","src","Constants.R")) +#source("Constants.R") + +## Get logger +c_logrMaaslin <- getLogger( "maaslin" ) + +funcDoGrubbs <- function( +### Use the Grubbs Test to identify outliers +iData, +### Column index in the data frame to test +frmeData, +### The data frame holding the data +dPOutlier, +### P-value threshold to indicate an outlier is significant +lsQC +### List holding the QC info of the cleaning step. Which indices are outliers is added. +){ + adData <- frmeData[,iData] + + # Original number of NA + viNAOrig = which(is.na(adData)) + + while( TRUE ) + { + lsTest <- try( grubbs.test( adData ), silent = TRUE ) + if( ( class( lsTest ) == "try-error" ) || is.na( lsTest$p.value ) || ( lsTest$p.value > dPOutlier ) ) + {break} + viOutliers = outlier( adData, logical = TRUE ) + adData[viOutliers] <- NA + } + + # Record removed data + viNAAfter = which(is.na(adData)) + + # If all were set to NA then ignore the filtering + if(length(adData)==length(viNAAfter)) + { + viNAAfter = viNAOrig + adData = frmeData[,iData] + c_logrMaaslin$info( paste("Grubbs Test:: Identifed all data as outliers so was inactived for index=",iData," data=",paste(as.vector(frmeData[,iData]),collapse=","), "number zeros=", length(which(frmeData[,iData]==0)), sep = " " )) + } else if(mean(adData, na.rm=TRUE) == 0) { + viNAAfter = viNAOrig + adData = frmeData[,iData] + c_logrMaaslin$info( paste("Grubbs Test::Removed all values but 0, ignored. Index=",iData,".",sep=" " ) ) + } else { + # Document removal + if( sum( is.na( adData ) ) ) + { + c_logrMaaslin$info( "Grubbs Test::Removing %d outliers from %s", sum( is.na( adData ) ), colnames(frmeData)[iData] ) + c_logrMaaslin$info( format( rownames( frmeData )[is.na( adData )] )) + } + } + + return(list(data=adData,outliers=length(viNAAfter)-length(viNAOrig),indices=setdiff(viNAAfter,viNAOrig))) +} + +funcDoFenceTest <- function( +### Use a threshold based on the quartiles of the data to identify outliers +iData, +### Column index in the data frame to test +frmeData, +### The data frame holding the data +dFence +### The fence outside the first and third quartiles to use as a threshold for cutt off. +### This many times the interquartile range +/- to the 3rd/1st quartiles +){ + # Establish fence + adData <- frmeData[,iData] + adQ <- quantile( adData, c(0.25, 0.5, 0.75), na.rm = TRUE ) + + dIQR <- adQ[3] - adQ[1] + if(!dIQR) + { + dIQR = sd(adData,na.rm = TRUE) + } + dUF <- adQ[3] + ( dFence * dIQR ) + dLF <- adQ[1] - ( dFence * dIQR ) + + # Record indices of values outside of fence to remove and remove. + aiRemove <- c() + for( j in 1:length( adData ) ) + { + d <- adData[j] + if( !is.na( d ) && ( ( d < dLF ) || ( d > dUF ) ) ) + { + aiRemove <- c(aiRemove, j) + } + } + + if(length(aiRemove)==length(adData)) + { + aiRemove = c() + c_logrMaaslin$info( "OutliersByFence:: Identified all data as outlier so was inactivated for index=", iData,"data=", paste(as.vector(frmeData[,iData]),collapse=","), "number zeros=", length(which(frmeData[,iData]==0)), sep=" " ) + } else { + adData[aiRemove] <- NA + + # Document to screen + if( length( aiRemove ) ) + { + c_logrMaaslin$info( "OutliersByFence::Removing %d outliers from %s", length( aiRemove ), colnames(frmeData)[iData] ) + c_logrMaaslin$info( format( rownames( frmeData )[aiRemove] )) + } + } + + return(list(data=adData,outliers=length(aiRemove),indices=aiRemove)) +} + +funcZerosAreUneven = function( +### +vdRawData, +### Raw data to be checked during transformation +funcTransform, +### Data transform to perform +vsStratificationFeatures, +### Groupings to check for unevenness +dfData +### Data frame holding the features +){ + # Return indicator of unevenness + fUneven = FALSE + + # Transform the data to compare + vdTransformed = funcTransform( vdRawData ) + + # Go through each stratification of data + for( sStratification in vsStratificationFeatures ) + { + # Current stratification + vFactorStrats = dfData[[ sStratification ]] + + # If the metadata is not a factor then skip + # Only binned data can be evaluated this way. + if( !is.factor( vFactorStrats )){ next } + + viZerosCountsRaw = c() + for( sLevel in levels( vFactorStrats ) ) + { + vdTest = vdRawData[ which( vFactorStrats == sLevel ) ] + viZerosCountsRaw = c( viZerosCountsRaw, length(which(vdTest == 0))) + vdTest = vdTransformed[ which( vFactorStrats == sLevel ) ] + } + dExpectation = 1 / length( viZerosCountsRaw ) + dMin = dExpectation / 2 + dMax = dExpectation + dMin + viZerosCountsRaw = viZerosCountsRaw / sum( viZerosCountsRaw ) + if( ( length( which( viZerosCountsRaw <= dMin ) ) > 0 ) || ( length( which( viZerosCountsRaw >= dMax ) ) > 0 ) ) + { + return( TRUE ) + } + } + return( fUneven ) +} + +funcTransformIncreasesOutliers = function( +### Checks if a data transform increases outliers in a distribution +vdRawData, +### Raw data to check for outlier zeros +funcTransform +){ + iUnOutliers = length( boxplot( vdRawData, plot = FALSE )$out ) + iTransformedOutliers = length( boxplot( funcTransform( vdRawData ), plot = FALSE )$out ) + + return( iUnOutliers <= iTransformedOutliers ) +} + +funcClean <- function( +### Properly clean / get data ready for analysis +### Includes custom analysis from the custom R script if it exists +frmeData, +### Data frame, input data to be acted on +funcDataProcess, +### Custom script that can be given to perform specialized processing before MaAsLin does. +aiMetadata, +### Indices of columns in frmeData which are metadata for analysis. +aiData, +### Indices of column in frmeData which are (abundance) data for analysis. +lsQCCounts, +### List that will hold the quality control information which is written in the output directory. +astrNoImpute = c(), +### An array of column names of frmeData not to impute. +dMinSamp, +### Minimum number of samples +dMinAbd, +# Minimum sample abundance +dFence, +### How many quartile ranges defines the fence to define outliers. +funcTransform, +### The data transformation function or a dummy function that does not affect the data +dPOutlier = 0.05 +### The significance threshold for the grubbs test to identify an outlier. +){ + # Call the custom script and set current data and indicies to the processes data and indicies. + c_logrMaaslin$debug( "Start Clean") + if( !is.null( funcDataProcess ) ) + { + c_logrMaaslin$debug("Additional preprocess function attempted.") + + pTmp <- funcDataProcess( frmeData=frmeData, aiMetadata=aiMetadata, aiData=aiData) + frmeData = pTmp$frmeData + aiMetadata = pTmp$aiMetadata + aiData = pTmp$aiData + lsQCCounts$lsQCCustom = pTmp$lsQCCounts + } + # Set data indicies after custom QC process. + lsQCCounts$aiAfterPreprocess = aiData + + # Remove missing data, remove any sample that has less than dMinSamp * the number of data or low abundance + aiRemove = c() + aiRemoveLowAbundance = c() + for( iCol in aiData ) + { + adCol = frmeData[,iCol] + adCol[!is.finite( adCol )] <- NA + if( ( sum( !is.na( adCol ) ) < ( dMinSamp * length( adCol ) ) ) || + ( length( unique( na.omit( adCol ) ) ) < 2 ) ) + { + aiRemove = c(aiRemove, iCol) + } + if( sum(adCol > dMinAbd, na.rm=TRUE ) < (dMinSamp * length( adCol))) + { + aiRemoveLowAbundance = c(aiRemoveLowAbundance, iCol) + } + } + # Remove and document + aiData = setdiff( aiData, aiRemove ) + aiData = setdiff( aiData, aiRemoveLowAbundance ) + lsQCCounts$iMissingData = aiRemove + lsQCCounts$iLowAbundanceData = aiRemoveLowAbundance + if(length(aiRemove)) + { + c_logrMaaslin$info( "Removing the following for data lower bound.") + c_logrMaaslin$info( format( colnames( frmeData )[aiRemove] )) + } + if(length(aiRemoveLowAbundance)) + { + c_logrMaaslin$info( "Removing the following for too many low abundance bugs.") + c_logrMaaslin$info( format( colnames( frmeData )[aiRemoveLowAbundance] )) + } + + #Transform data + iTransformed = 0 + viNotTransformedData = c() + for(aiDatum in aiData) + { + adValues = frmeData[,aiDatum] +# if( ! funcTransformIncreasesOutliers( adValues, funcTransform ) ) +# { + frmeData[,aiDatum] = funcTransform( adValues ) +# iTransformed = iTransformed + 1 +# } else { +# viNotTransformedData = c( viNotTransformedData, aiDatum ) +# } + } + c_logrMaaslin$info(paste("Number of features transformed = ",iTransformed)) + + # Metadata: Properly factorize all logical data and integer and number data with less than iNonFactorLevelThreshold + # Also record which are numeric metadata + aiNumericMetadata = c() + for( i in aiMetadata ) + { + if( ( class( frmeData[,i] ) %in% c("integer", "numeric", "logical") ) && + ( length( unique( frmeData[,i] ) ) < c_iNonFactorLevelThreshold ) ) { + c_logrMaaslin$debug(paste("Changing metadatum from numeric/integer/logical to factor",colnames(frmeData)[i],sep="=")) + frmeData[,i] = factor( frmeData[,i] ) + } + if( class( frmeData[,i] ) %in% c("integer","numeric") ) + { + aiNumericMetadata = c(aiNumericMetadata,i) + } + } + + # Remove outliers + # If the dFence Value is set use the method of defining the outllier as + # dFence * the interquartile range + or - the 3rd and first quartile respectively. + # If not the gibbs test is used. + lsQCCounts$aiDataSumOutlierPerDatum = c() + lsQCCounts$aiMetadataSumOutlierPerDatum = c() + lsQCCounts$liOutliers = list() + + if( dFence > 0.0 ) + { + # For data + for( iData in aiData ) + { + lOutlierInfo <- funcDoFenceTest(iData=iData,frmeData=frmeData,dFence=dFence) + frmeData[,iData] <- lOutlierInfo[["data"]] + lsQCCounts$aiDataSumOutlierPerDatum <- c(lsQCCounts$aiDataSumOutlierPerDatum,lOutlierInfo[["outliers"]]) + if(lOutlierInfo[["outliers"]]>0) + { + lsQCCounts$liOutliers[[paste(iData,sep="")]] <- lOutlierInfo[["indices"]] + } + } + + # Remove outlier non-factor metadata + for( iMetadata in aiNumericMetadata ) + { + lOutlierInfo <- funcDoFenceTest(iData=iMetadata,frmeData=frmeData,dFence=dFence) + frmeData[,iMetadata] <- lOutlierInfo[["data"]] + lsQCCounts$aiMetadataSumOutlierPerDatum <- c(lsQCCounts$aiMetadataSumOutlierPerDatum,lOutlierInfo[["outliers"]]) + if(lOutlierInfo[["outliers"]]>0) + { + lsQCCounts$liOutliers[[paste(iMetadata,sep="")]] <- lOutlierInfo[["indices"]] + } + } + #Do not use the fence, use the Grubbs test + } else if(dPOutlier!=0.0){ + # For data + for( iData in aiData ) + { + lOutlierInfo <- funcDoGrubbs(iData=iData,frmeData=frmeData,dPOutlier=dPOutlier) + frmeData[,iData] <- lOutlierInfo[["data"]] + lsQCCounts$aiDataSumOutlierPerDatum <- c(lsQCCounts$aiDataSumOutlierPerDatum,lOutlierInfo[["outliers"]]) + if(lOutlierInfo[["outliers"]]>0) + { + lsQCCounts$liOutliers[[paste(iData,sep="")]] <- lOutlierInfo[["indices"]] + } + } + for( iMetadata in aiNumericMetadata ) + { + lOutlierInfo <- funcDoGrubbs(iData=iMetadata,frmeData=frmeData,dPOutlier=dPOutlier) + frmeData[,iMetadata] <- lOutlierInfo[["data"]] + lsQCCounts$aiMetadataSumOutlierPerDatum <- c(lsQCCounts$aiMetadataSumOutlierPerDatum,lOutlierInfo[["outliers"]]) + if(lOutlierInfo[["outliers"]]>0) + { + lsQCCounts$liOutliers[[paste(iMetadata,sep="")]] <- lOutlierInfo[["indices"]] + } + } + } + + # Metadata: Remove missing data + # This is defined as if there is only one non-NA value or + # if the number of NOT NA data is less than a percentage of the data defined by dMinSamp + aiRemove = c() + for( iCol in c(aiMetadata) ) + { + adCol = frmeData[,iCol] + if( ( sum( !is.na( adCol ) ) < ( dMinSamp * length( adCol ) ) ) || + ( length( unique( na.omit( adCol ) ) ) < 2 ) ) + { + aiRemove = c(aiRemove, iCol) + } + } + + # Remove metadata + aiMetadata = setdiff( aiMetadata, aiRemove ) + + # Update the data which was removed. + lsQCCounts$iMissingMetadata = aiRemove + if(length(aiRemove)) + { + c_logrMaaslin$info("Removing the following metadata for too much missing data or only one data value outside of NA.") + c_logrMaaslin$info(format(colnames( frmeData )[aiRemove])) + } + + # Keep track of factor levels in a list for later use + lslsFactors <- list() + for( iCol in c(aiMetadata) ) + { + aCol <- frmeData[,iCol] + if( class( aCol ) == "factor" ) + { + lslsFactors[[length( lslsFactors ) + 1]] <- list(iCol, levels( aCol )) + } + } + + # Replace missing data values by the mean of the data column. + # Remove samples that were all NA from the cleaning and so could not be imputed. + aiRemoveData = c() + for( iCol in aiData ) + { + adCol <- frmeData[,iCol] + adCol[is.infinite( adCol )] <- NA + adCol[is.na( adCol )] <- mean( adCol[which(adCol>0)], na.rm = TRUE ) + frmeData[,iCol] <- adCol + + if(length(which(is.na(frmeData[,iCol]))) == length(frmeData[,iCol])) + { + c_logrMaaslin$info( paste("Removing data", iCol, "for being all NA after QC")) + aiRemoveData = c(aiRemoveData,iCol) + } + } + + # Remove and document + aiData = setdiff( aiData, aiRemoveData ) + lsQCCounts$iMissingData = c(lsQCCounts$iMissingData,aiRemoveData) + if(length(aiRemoveData)) + { + c_logrMaaslin$info( "Removing the following for having only NAs after cleaning (maybe due to only having NA after outlier testing).") + c_logrMaaslin$info( format( colnames( frmeData )[aiRemoveData] )) + } + + #Use na.gam.replace to manage NA metadata + aiTmp <- setdiff( aiMetadata, which( colnames( frmeData ) %in% astrNoImpute ) ) + # Keep tack of NAs so the may not be plotted later. + liNaIndices = list() + lsNames = names(frmeData) + for( i in aiTmp) + { + liNaIndices[[lsNames[i]]] = which(is.na(frmeData[,i])) + } + frmeData[,aiTmp] <- na.gam.replace( frmeData[,aiTmp] ) + + #If NA is a value in factor data, set the NA as a level. + for( lsFactor in lslsFactors ) + { + iCol <- lsFactor[[1]] + aCol <- frmeData[,iCol] + if( "NA" %in% levels( aCol ) ) + { + if(! lsNames[iCol] %in% astrNoImpute) + { + liNaIndices[[lsNames[iCol]]] = union(which(is.na(frmeData[,iCol])),which(frmeData[,iCol]=="NA")) + } + frmeData[,iCol] <- factor( aCol, levels = c(lsFactor[[2]], "NA") ) + } + } + + # Make sure there is a minimum number of non-0 measurements + aiRemove = c() + for( iCol in aiData ) + { + adCol = frmeData[,iCol] + if(length( which(adCol!=0)) < ( dMinSamp * length( adCol ) ) ) + { + aiRemove = c(aiRemove, iCol) + } + } + + # Remove and document + aiData = setdiff( aiData, aiRemove) + lsQCCounts$iZeroDominantData = aiRemove + if(length(aiRemove)) + { + c_logrMaaslin$info( "Removing the following for having not enough non-zero measurments for analysis.") + c_logrMaaslin$info( format( colnames( frmeData )[aiRemove] )) + } + + c_logrMaaslin$debug("End FuncClean") + return( list(frmeData = frmeData, aiMetadata = aiMetadata, aiData = aiData, lsQCCounts = lsQCCounts, liNaIndices=liNaIndices, viNotTransformedData = viNotTransformedData) ) + ### Return list of + ### frmeData: The Data after cleaning + ### aiMetadata: The indices of the metadata still being used after filtering + ### aiData: The indices of the data still being used after filtering + ### lsQCCOunts: QC info +} + +funcBugs <- function( +### Run analysis of all data features against all metadata +frmeData, +### Cleaned data including metadata, and data +lsData, +### This list is a general container for data as the analysis occurs, think about it as a cache for the analysis +aiMetadata, +### Indices of metadata used in analysis +aiData, +### Indices of response data +aiNotTransformedData, +### Indicies of the data not transformed +strData, +### Log file name +dSig, +### Significance threshold for the qvalue cut off +fInvert=FALSE, +### Invert images to have a black background +strDirOut = NA, +### Output project directory +funcReg=NULL, +### Function for regularization +funcTransform=NULL, +### Function used to transform the data +funcUnTransform=NULL, +### If a transform is used the opposite of that transfor must be used on the residuals in the partial residual plots +lsNonPenalizedPredictors=NULL, +### These predictors will not be penalized in the feature (model) selection step +funcAnalysis=NULL, +### Function to perform association analysis +lsRandomCovariates=NULL, +### List of string names of metadata which will be treated as random covariates +funcGetResults=NULL, +### Function to unpack results from analysis +fDoRPlot=TRUE, +### Plot residuals +fOmitLogFile = FALSE, +### Stops the creation of the log file +fAllvAll=FALSE, +### Flag to turn on all against all comparisons +liNaIndices = list(), +### Indicies of imputed NA data +lxParameters=list(), +### List holds parameters for different variable selection techniques +strTestingCorrection = "BH", +### Correction for multiple testing +fIsUnivariate = FALSE, +### Indicates if the function is univariate +fZeroInflated = FALSE +### Indicates to use a zero infalted model +){ + c_logrMaaslin$debug("Start funcBugs") + # If no output directory is indicated + # Then make it the current directory + if( is.na( strDirOut ) || is.null( strDirOut ) ) + { + if( !is.na( strData ) ) + { + strDirOut <- paste( dirname( strData ), "/", sep = "" ) + } else { strDirOut = "" } + } + + # Make th log file and output file names based on the log file name + strLog = NA + strBase = "" + if(!is.na(strData)) + { + strBaseOut <- paste( strDirOut, sub( "\\.([^.]+)$", "", basename(strData) ), sep = "/" ) + strLog <- paste( strBaseOut,c_sLogFileSuffix, ".txt", sep = "" ) + } + + # If indicated, stop the creation of the log file + # Otherwise delete the log file if it exists and log + if(fOmitLogFile){ strLog = NA } + if(!is.na(strLog)) + { + c_logrMaaslin$info( "Outputting to: %s", strLog ) + unlink( strLog ) + } + + # Will contain pvalues + adP = c() + adPAdj = c() + + # List of lists with association information + lsSig <- list() + # Go through each data that was not previously removed and perform inference + for( iTaxon in aiData ) + { + # Log to screen progress per 10 associations. + # Can be thown off if iTaxon is missing a mod 10 value + # So the taxons may not be logged every 10 but not a big deal + if( !( iTaxon %% 10 ) ) + { + c_logrMaaslin$info( "Taxon %d/%d", iTaxon, max( aiData ) ) + } + + # Call analysis method + lsOne <- funcBugHybrid( iTaxon=iTaxon, frmeData=frmeData, lsData=lsData, aiMetadata=aiMetadata, dSig=dSig, adP=adP, lsSig=lsSig, funcTransform=funcTransform, funcUnTransform=funcUnTransform, strLog=strLog, funcReg=funcReg, lsNonPenalizedPredictors=lsNonPenalizedPredictors, funcAnalysis=funcAnalysis, lsRandomCovariates=lsRandomCovariates, funcGetResult=funcGetResults, fAllvAll=fAllvAll, fIsUnivariate=fIsUnivariate, lxParameters=lxParameters, fZeroInflated=fZeroInflated, fIsTransformed= ! iTaxon %in% aiNotTransformedData ) + + # If you get a NA (happens when the lmm gets all random covariates) move on + if( is.na( lsOne ) ){ next } + + # The updating of the following happens in the inference method call in the funcBugHybrid call + # New pvalue array + adP <- lsOne$adP + # New lsSig contains data about significant feature v metadata comparisons + lsSig <- lsOne$lsSig + # New qc data + lsData$lsQCCounts = lsOne$lsQCCounts + } + + # Log the QC info + c_logrMaaslin$debug("lsData$lsQCCounts") + c_logrMaaslin$debug(format(lsData$lsQCCounts)) + + if( is.null( adP ) ) { return( NULL ) } + + # Perform bonferonni corrections on factor data (for levels), calculate the number of tests performed, and FDR adjust for multiple hypotheses + # Perform Bonferonni adjustment on factor data + for( iADIndex in 1:length( adP ) ) + { + # Only perform on factor data + if( is.factor( lsSig[[ iADIndex ]]$metadata ) ) + { + adPAdj = c( adPAdj, funcBonferonniCorrectFactorData( dPvalue = adP[ iADIndex ], vsFactors = lsSig[[ iADIndex ]]$metadata, fIgnoreNAs = length(liNaIndices)>0) ) + } else { + adPAdj = c( adPAdj, adP[ iADIndex ] ) + } + } + + iTests = funcCalculateTestCounts(iDataCount = length(aiData), asMetadata = intersect( lsData$astrMetadata, colnames( frmeData )[aiMetadata] ), asForced = lsNonPenalizedPredictors, asRandom = lsRandomCovariates, fAllvAll = fAllvAll) + + #Get indices of sorted data after the factor correction but before the multiple hypothesis corrections. + aiSig <- sort.list( adPAdj ) + + # Perform FDR BH + adQ = p.adjust(adPAdj, method=strTestingCorrection, n=max(length(adPAdj), iTests)) + + # Find all covariates that had significant associations + astrNames <- c() + for( i in 1:length( lsSig ) ) + { + astrNames <- c(astrNames, lsSig[[i]]$name) + } + astrNames <- unique( astrNames ) + + # Sets up named label return for global plotting + lsReturnTaxa <- list() + for( j in aiSig ) + { + if( adQ[j] > dSig ) { next } + strTaxon <- lsSig[[j]]$taxon + if(strTaxon %in% names(lsReturnTaxa)) + { + lsReturnTaxa[[strTaxon]] = min(lsReturnTaxa[[strTaxon]],adQ[j]) + } else { lsReturnTaxa[[strTaxon]] = adQ[j]} + } + + # For each covariate with significant associations + # Write out a file with association information + for( strName in astrNames ) + { + strFileTXT <- NA + strFilePDF <- NA + for( j in aiSig ) + { + lsCur <- lsSig[[j]] + strCur <- lsCur$name + + if( strCur != strName ) { next } + + strTaxon <- lsCur$taxon + adData <- lsCur$data + astrFactors <- lsCur$factors + adCur <- lsCur$metadata + adY <- adData + + if( is.na( strData ) ) { next } + + ## If the text file output is not written to yet + ## make the file names, and delete any previous file output + if( is.na( strFileTXT ) ) + { + strFileTXT <- sprintf( "%s-%s.txt", strBaseOut, strName ) + unlink(strFileTXT) + funcWrite( c("Variable", "Feature", "Value", "Coefficient", "N", "N not 0", "P-value", "Q-value"), strFileTXT ) + } + + ## Write text output + funcWrite( c(strName, strTaxon, lsCur$orig, lsCur$value, length( adData ), sum( adData > 0 ), adP[j], adQ[j]), strFileTXT ) + + ## If the significance meets the threshold + ## Write PDF file output + if( adQ[j] > dSig ) { next } + + # Do not make residuals plots if univariate is selected + strFilePDF = funcPDF( frmeTmp=frmeData, lsCur=lsCur, curPValue=adP[j], curQValue=adQ[j], strFilePDF=strFilePDF, strBaseOut=strBaseOut, strName=strName, funcUnTransform= funcUnTransform, fDoResidualPlot=fDoRPlot, fInvert=fInvert, liNaIndices=liNaIndices ) + } + if( dev.cur( ) != 1 ) { dev.off( ) } + } + aiTmp <- aiData + + logdebug("End funcBugs", c_logMaaslin) + return(list(lsReturnBugs=lsReturnTaxa, lsQCCounts=lsData$lsQCCounts)) + ### List of data features successfully associated without error and quality control data +} + +#Lightly Tested +### Performs analysis for 1 feature +### iTaxon: integer Taxon index to be associated with data +### frmeData: Data frame The full data +### lsData: List of all associated data +### aiMetadata: Numeric vector of indices +### dSig: Numeric significance threshold for q-value cut off +### adP: List of pvalues from associations +### lsSig: List which serves as a cache of data about significant associations +### strLog: String file to log to +funcBugHybrid <- function( +### Performs analysis for 1 feature +iTaxon, +### integer Taxon index to be associated with data +frmeData, +### Data frame, the full data +lsData, +### List of all associated data +aiMetadata, +### Numeric vector of indices +dSig, +### Numeric significance threshold for q-value cut off +adP, +### List of pvalues from associations +lsSig, +### List which serves as a cache of data about significant associations +funcTransform, +### The tranform used on the data +funcUnTransform, +### The reverse transform on the data +strLog = NA, +### String, file to which to log +funcReg=NULL, +### Function to perform regularization +lsNonPenalizedPredictors=NULL, +### These predictors will not be penalized in the feature (model) selection step +funcAnalysis=NULL, +### Function to perform association analysis +lsRandomCovariates=NULL, +### List of string names of metadata which will be treated as random covariates +funcGetResult=NULL, +### Function to unpack results from analysis +fAllvAll=FALSE, +### Flag to turn on all against all comparisons +fIsUnivariate = FALSE, +### Indicates the analysis function is univariate +lxParameters=list(), +### List holds parameters for different variable selection techniques +fZeroInflated = FALSE, +### Indicates if to use a zero infalted model +fIsTransformed = TRUE +### Indicates that the bug is transformed +){ +#dTime00 <- proc.time()[3] + #Get metadata column names + astrMetadata = intersect( lsData$astrMetadata, colnames( frmeData )[aiMetadata] ) + + #Get data measurements that are not NA + aiRows <- which( !is.na( frmeData[,iTaxon] ) ) + + #Get the dataframe of non-na data measurements + frmeTmp <- frmeData[aiRows,] + + #Set the min boosting selection frequency to a default if not given + if( is.na( lxParameters$dFreq ) ) + { + lxParameters$dFreq <- 0.5 / length( c(astrMetadata) ) + } + + # Get the full data for the bug feature + adCur = frmeTmp[,iTaxon] + lxParameters$sBugName = names(frmeTmp[iTaxon]) + + # This can run multiple models so some of the results are held in lists and some are not + llmod = list() + liTaxon = list() + lastrTerms = list() + + # Build formula for simple mixed effects models + # Removes random covariates from variable selection + astrMetadata = setdiff(astrMetadata, lsRandomCovariates) + strFormula <- paste( "adCur ~", paste( sprintf( "`%s`", astrMetadata ), collapse = " + " ), sep = " " ) + + # Document the model + funcWrite( c("#taxon", colnames( frmeTmp )[iTaxon]), strLog ) + funcWrite( c("#metadata", astrMetadata), strLog ) + funcWrite( c("#samples", rownames( frmeTmp )), strLog ) + + #Model terms + astrTerms <- c() + + # Attempt feature (model) selection + if(!is.na(funcReg)) + { + #Count model selection method attempts + lsData$lsQCCounts$iBoosts = lsData$lsQCCounts$iBoosts + 1 + #Perform model selection + astrTerms <- funcReg(strFormula=strFormula, frmeTmp=frmeTmp, adCur=adCur, lsParameters=lxParameters, lsForcedParameters=lsNonPenalizedPredictors, strLog=strLog) + #If the feature selection function is set to None, set all terms of the model to all the metadata + } else { astrTerms = astrMetadata } + + # Get look through the boosting results to get a model + # Holds the predictors in the predictors in the model that were selected by the boosting + if(is.null( astrTerms )){lsData$lsQCCounts$iBoostErrors = lsData$lsQCCounts$iBoostErrors + 1} + + # Get the indices that are transformed + # Of those indices check for uneven metadata + # Untransform any of the metadata that failed + # Failed means true for uneven occurences of zeros +# if( fIsTransformed ) +# { +# vdUnevenZeroCheck = funcUnTransform( frmeData[[ iTaxon ]] ) +# if( funcZerosAreUneven( vdRawData=vdUnevenZeroCheck, funcTransform=funcTransform, vsStratificationFeatures=astrTerms, dfData=frmeData ) ) +# { +# frmeData[[ iTaxon ]] = vdUnevenZeroCheck +# c_logrMaaslin$debug( paste( "Taxon transformation reversed due to unevenness of zero distribution.", iTaxon ) ) +# } +# } + + # Run association analysis if predictors exist and an analysis function is specified + # Run analysis + if(!is.na(funcAnalysis) ) + { + #If there are selected and forced fixed covariates + if( length( astrTerms ) ) + { + #Count the association attempt + lsData$lsQCCounts$iLms = lsData$lsQCCounts$iLms + 1 + + #Make the lm formula + #Build formula for simple mixed effects models using random covariates + strRandomCovariatesFormula = NULL + #Random covariates are forced + if(length(lsRandomCovariates)>0) + { + #Format for lme + #Needed for changes to not allowing random covariates through the boosting process + strRandomCovariatesFormula <- paste( "adCur ~ ", paste( sprintf( "1|`%s`", lsRandomCovariates), collapse = " + " )) + } + + #Set up a list of formula containing selected fixed variables changing and the forced fixed covariates constant + vstrFormula = c() + #Set up suppressing forced covariates in a all v all scenario only + asSuppress = c() + #Enable all against all comparisons + if(fAllvAll && !fIsUnivariate) + { + lsVaryingCovariates = setdiff(astrTerms,lsNonPenalizedPredictors) + lsConstantCovariates = setdiff(lsNonPenalizedPredictors,lsRandomCovariates) + strConstantFormula = paste( sprintf( "`%s`", lsConstantCovariates ), collapse = " + " ) + asSuppress = lsConstantCovariates + + if(length(lsVaryingCovariates)==0L) + { + vstrFormula <- c( paste( "adCur ~ ", paste( sprintf( "`%s`", lsConstantCovariates ), collapse = " + " )) ) + } else { + for( sVarCov in lsVaryingCovariates ) + { + strTempFormula = paste( "adCur ~ `", sVarCov,"`",sep="") + if(length(lsConstantCovariates)>0){ strTempFormula = paste(strTempFormula,strConstantFormula,sep=" + ") } + vstrFormula <- c( vstrFormula, strTempFormula ) + } + } + } else { + #This is either the multivariate case formula for all covariates in an lm or fixed covariates in the lmm + vstrFormula <- c( paste( "adCur ~ ", paste( sprintf( "`%s`", astrTerms ), collapse = " + " )) ) + } + + #Run the association + for( strAnalysisFormula in vstrFormula ) + { + i = length(llmod)+1 + llmod[[i]] = funcAnalysis(strFormula=strAnalysisFormula, frmeTmp=frmeTmp, iTaxon=iTaxon, lsHistory=list(adP=adP, lsSig=lsSig, lsQCCounts=lsData$lsQCCounts), strRandomFormula=strRandomCovariatesFormula, fZeroInflated=fZeroInflated) + + liTaxon[[i]] = iTaxon + lastrTerms[[i]] = funcFormulaStrToList(strAnalysisFormula) + } + } else { + #If there are no selected or forced fixed covariates + lsData$lsQCCounts$iNoTerms = lsData$lsQCCounts$iNoTerms + 1 + return(list(adP=adP, lsSig=lsSig, lsQCCounts=lsData$lsQCCounts)) + } + } + + #Call funcBugResults and return it's return + if(!is.na(funcGetResult)) + { + #Format the results to a consistent expected result. + return( funcGetResult( llmod=llmod, frmeData=frmeData, liTaxon=liTaxon, dSig=dSig, adP=adP, lsSig=lsSig, strLog=strLog, lsQCCounts=lsData$lsQCCounts, lastrCols=lastrTerms, asSuppressCovariates=asSuppress ) ) + } else { + return(list(adP=adP, lsSig=lsSig, lsQCCounts=lsData$lsQCCounts)) + } + ### List containing a list of pvalues, a list of significant data per association, and a list of QC data +}