0
+ − 1 #!/usr/bin/env python
+ − 2 """
+ − 3 Author: Timothy Tickle
+ − 4 Description: Performs common manipulations on tables
+ − 5 """
+ − 6
+ − 7 __author__ = "Timothy Tickle"
+ − 8 __copyright__ = "Copyright 2012"
+ − 9 __credits__ = ["Timothy Tickle"]
+ − 10 __license__ = ""
+ − 11 __version__ = ""
+ − 12 __maintainer__ = "Timothy Tickle"
+ − 13 __email__ = "ttickle@sph.harvard.edu"
+ − 14 __status__ = "Development"
+ − 15
+ − 16 import argparse
+ − 17 import csv
+ − 18 import sys
+ − 19 import re
+ − 20 import os
+ − 21 import numpy as np
+ − 22 from src.AbundanceTable import AbundanceTable
+ − 23 #from src.PCA import PCA
+ − 24 from src.ValidateData import ValidateData
+ − 25
+ − 26 #Set up arguments reader
+ − 27 argp = argparse.ArgumentParser( prog = "scriptManipulateTable.py",
+ − 28 description = """Performs common manipulations on tables.\nExample: python scriptManipulateTable.py -i TID -l STSite Test.pcl""" )
+ − 29
+ − 30 #Arguments
+ − 31 #Describe table
+ − 32 argp.add_argument("-i","--id", dest="sIDName", default="ID", help="Abundance Table ID")
+ − 33 argp.add_argument("-l","--meta", dest="sLastMetadataName", help="Last metadata name")
+ − 34 argp.add_argument("-d","--fileDelim", dest= "cFileDelimiter", action= "store", default="\t", help="File delimiter, default tab")
+ − 35 argp.add_argument("-f","--featureDelim", dest= "cFeatureDelimiter", action= "store", default="|", help="Feature (eg. bug or function) delimiter, default '|'")
+ − 36
+ − 37 #Checked x 2
+ − 38 argp.add_argument("-n","--doNorm", dest="fNormalize", action="store_true", default=False, help="Flag to turn on normalization")
+ − 39 argp.add_argument("-s","--doSum", dest="fSum", action="store_true", default=False, help="Flag to turn on summation")
+ − 40
+ − 41 #Unsupervised filtering
+ − 42 argp.add_argument("-A","--doFilterAbundance", dest="strFilterAbundance", action="store", default=None, help="Turns on filtering by abundance (remove features that do not have a minimum abundance in a minimum number of samples); Should be a real number and an integer in the form 'minAbundance,minSamples', (should be performed on a normalized file).")
+ − 43 argp.add_argument("-P","--doFilterPercentile", dest="strFilterPercentile", action="store", default=None, help="Turns on filtering by percentile Should be two numbers between 0 and 1 in the form 'percentile,percentage'. (should be performed on a normalized file).")
+ − 44 argp.add_argument("-O","--doFilterOccurrence", dest="strFilterOccurence", action="store", default=None, help="Turns on filtering by occurrence. Should be two integers in the form 'minSequence,minSample' (should NOT be performed on a normalized file).")
+ − 45 #argp.add_argument("-D","--doFilterDeviation", dest="dCuttOff", action="store", type=float, default=None, help="Flag to turn on filtering by standard deviation (should NOT be performed on a normalized file).")
+ − 46
+ − 47 #Change bug membership
+ − 48 argp.add_argument("-t","--makeTerminal", dest="fMakeTerminal", action="store_true", default=False, help="Works reduces the file to teminal features in the original file.")
+ − 49 argp.add_argument("-u","--reduceOTUs", dest="fRemoveOTUs", action="store_true", default=False, help="Remove otu entries from file.")
+ − 50 argp.add_argument("-c","--reduceToClade", dest="iClade", action="store", type=int, default=None, help="Specify a level of clade to reduce to [].")
+ − 51 argp.add_argument("-b","--reduceToFeatures", dest="strFeatures", action="store", default=None, help="Reduce measurements to certain features (bugs or functions). This can be a comma delimited string (of atleast 2 bugs) or a file.")
+ − 52
+ − 53 #Manipulate based on metadata
+ − 54 #Checked
+ − 55 argp.add_argument("-y","--stratifyBy", dest="strStratifyBy", action="store", default=None, help="Metadata to stratify tables by.")
+ − 56 argp.add_argument("-r","--removeMetadata", dest="strRemoveMetadata", action="store", default=None, help="Remove samples of this metadata and value (format comma delimited string with metadata id first and the values to remove after 'id,lvalue1,value2').")
+ − 57
+ − 58 #Manipulate lineage
+ − 59 #Checked
+ − 60 argp.add_argument("-x","--doPrefixClades", dest="fPrefixClades", action="store_true", default=False, help="Flag to turn on adding prefixes to clades to better identify them, for example s__ will be placed infront of each species.")
+ − 61
+ − 62 #Combine tables
+ − 63 #argp.add_argument("-m","--combineIntersect", dest="fCombineIntersect", action="store_true", default=False, help="Combine two tables including only common features/metadata (intersection).")
+ − 64 #argp.add_argument("-e","--combineUnion", dest="fCombineUnion", action="store_true", default=False, help="Combine two tables (union).")
+ − 65
+ − 66 #Dimensionality Reduction
+ − 67 #argp.add_argument("-p","--doPCA", dest="fDoPCA",action="store_true", default=False, help="Flag to turn on adding metabugs and metametadata by performing PCA on each of bug relative abundance and continuous metadata and add the resulting components")
+ − 68
+ − 69 #Checked
+ − 70 argp.add_argument("-o","--output", dest="strOutFile", action="store", default=None, help="Indicate output pcl file.")
+ − 71 argp.add_argument("strFileAbund", help ="Input data file")
+ − 72
+ − 73 args = argp.parse_args( )
+ − 74
+ − 75 # Creat output file if needed.
+ − 76 if not args.strOutFile:
+ − 77 args.strOutFile = os.path.splitext(args.strFileAbund)[0]+"-mod.pcl"
+ − 78 lsPieces = os.path.splitext(args.strOutFile)
+ − 79
+ − 80 #List of abundance tables
+ − 81 lsTables = []
+ − 82
+ − 83 #Read in abundance table
+ − 84 abndTable = AbundanceTable.funcMakeFromFile(xInputFile=args.strFileAbund,
+ − 85 cDelimiter = args.cFileDelimiter,
+ − 86 sMetadataID = args.sIDName,
+ − 87 sLastMetadata = args.sLastMetadataName,
+ − 88 lOccurenceFilter = None,
+ − 89 cFeatureNameDelimiter=args.cFeatureDelimiter,
+ − 90 xOutputFile = args.strOutFile)
+ − 91
+ − 92 #TODO Check filtering, can not have some filtering together
+ − 93
+ − 94 # Make feature list
+ − 95 lsFeatures = []
+ − 96 if args.strFeatures:
+ − 97 print "Get features not completed"
+ − 98 # if "," in args.strFeatures:
+ − 99 # lsFeatures = args.strFeatures.split(",")
+ − 100 # print "ManipulateTable::Reading in feature list "+str(len(lsFeatures))+"."
+ − 101 # else:
+ − 102 # csvr = csv.reader(open(args.strFeatures, "rU"))
+ − 103 # print "ManipulateTable::Reading in feature file "+args.strFeatures+"."
+ − 104 # for lsLine in csvr:
+ − 105 # lsFeatures.extend(lsLine)
+ − 106
+ − 107 lsTables.append(abndTable)
+ − 108
+ − 109 # Do summing
+ − 110 #Sum if need
+ − 111 if args.fSum:
+ − 112 for abndTable in lsTables:
+ − 113 print "ManipulateTable::"+abndTable.funcGetName()+" had "+str(len(abndTable.funcGetFeatureNames()))+" features before summing."
+ − 114 fResult = abndTable.funcSumClades()
+ − 115 if fResult:
+ − 116 print "ManipulateTable::"+abndTable.funcGetName()+" was summed."
+ − 117 print "ManipulateTable::"+abndTable.funcGetName()+" has "+str(len(abndTable.funcGetFeatureNames()))+" features after summing."
+ − 118 else:
+ − 119 print "ManipulateTable::ERROR. "+abndTable.funcGetName()+" was NOT summed."
+ − 120
+ − 121 # Filter on counts
+ − 122 if args.strFilterOccurence:
+ − 123 iMinimumSequence,iMinimumSample = args.strFilterOccurence.split(",")
+ − 124 for abndTable in lsTables:
+ − 125 if abndTable.funcIsNormalized():
+ − 126 print "ManipulateTable::ERROR. "+abndTable.funcGetName()+" is normalized and can not be filtered by occurence. This filter needs counts."
+ − 127 else:
+ − 128 fResult = abndTable.funcFilterAbundanceBySequenceOccurence(iMinSequence = int(iMinimumSequence), iMinSamples = int(iMinimumSample))
+ − 129 if fResult:
+ − 130 print "ManipulateTable::"+abndTable.funcGetName()+" was filtered by occurence and now has "+str(len(abndTable.funcGetFeatureNames()))+" features."
+ − 131 else:
+ − 132 print "ManipulateTable::ERROR. "+abndTable.funcGetName()+" was NOT filtered by occurence."
+ − 133
+ − 134 # Change bug membership
+ − 135 if args.fMakeTerminal:
+ − 136 lsTerminalTables = []
+ − 137 for abndTable in lsTables:
+ − 138 print "ManipulateTable::"+abndTable.funcGetName()+" had "+str(len(abndTable.funcGetFeatureNames()))+" features before making terminal."
+ − 139 abndTable = abndTable.funcGetFeatureAbundanceTable(abndTable.funcGetTerminalNodes())
+ − 140 if abndTable:
+ − 141 print "ManipulateTable::"+abndTable.funcGetName()+" has "+str(len(abndTable.funcGetFeatureNames()))+" terminal features."
+ − 142 lsTerminalTables.append(abndTable)
+ − 143 else:
+ − 144 print "ManipulateTable::ERROR. "+abndTable.funcGetName()+" was not made terminal."
+ − 145 lsTables = lsTerminalTables
+ − 146
+ − 147 if args.fRemoveOTUs:
+ − 148 lsNotOTUs = []
+ − 149 for abndTable in lsTables:
+ − 150 print "ManipulateTable::"+abndTable.funcGetName()+" had "+str(len(abndTable.funcGetFeatureNames()))+" features before removing OTUs."
+ − 151 abndTable = abndTable.funcGetWithoutOTUs()
+ − 152 if abndTable:
+ − 153 print "ManipulateTable::"+abndTable.funcGetName()+" had OTUs removed and now has "+str(len(abndTable.funcGetFeatureNames()))+" features."
+ − 154 lsNotOTUs.append(abndTable)
+ − 155 else:
+ − 156 print "ManipulateTable::ERROR. "+abndTable.funcGetName()+" OTUs were not removed."
+ − 157 lsTables = lsNotOTUs
+ − 158
+ − 159 if args.iClade:
+ − 160 for abndTable in lsTables:
+ − 161 fResult = abndTable.funcReduceFeaturesToCladeLevel(args.iClade)
+ − 162 if fResult:
+ − 163 print "ManipulateTable::"+abndTable.funcGetName()+" was reduced to clade level "+str(args.iClade)+"."
+ − 164 else:
+ − 165 print "ManipulateTable::ERROR. "+abndTable.funcGetName()+" was NOT reduced in clade levels."
+ − 166
+ − 167 if args.strFeatures:
+ − 168 for abndTable in lsTables:
+ − 169 fResult = abndTable.funcGetFeatureAbundanceTable(lsFeatures)
+ − 170 if fResult:
+ − 171 print "ManipulateTable::"+abndTable.funcGetName()+" has been reduced to given features and now has "+str(len(abndTable.funcGetFeatureNames()))+" features."
+ − 172 else:
+ − 173 print "ManipulateTable::ERROR. "+abndTable.funcGetName()+" could not be reduced to the given list."
+ − 174
+ − 175 if args.strRemoveMetadata:
+ − 176 lsMetadata = args.strRemoveMetadata.split(",")
+ − 177 for abndTable in lsTables:
+ − 178 fResult = abndTable.funcRemoveSamplesByMetadata(sMetadata=lsMetadata[0], lValuesToRemove=lsMetadata[1:])
+ − 179 if fResult:
+ − 180 print "ManipulateTable::"+abndTable.funcGetName()+" has had samples removed and now has "+str(len(abndTable.funcGetSampleNames()))+" samples."
+ − 181 else:
+ − 182 print "ManipulateTable::ERROR. Could not remove samples from "+abndTable.funcGetName()+"."
+ − 183
+ − 184 # Normalize if needed
+ − 185 if args.fNormalize:
+ − 186 for abndTable in lsTables:
+ − 187 fResult = abndTable.funcNormalize()
+ − 188 if fResult:
+ − 189 print "ManipulateTable::"+abndTable.funcGetName()+" was normalized."
+ − 190 else:
+ − 191 print "ManipulateTable::"+abndTable.funcGetName()+" was NOT normalized."
+ − 192
+ − 193 # Filter on percentile
+ − 194 if args.strFilterPercentile:
+ − 195 dPercentile,dPercentage = args.strFilterPercentile.split(",")
+ − 196 for abndTable in lsTables:
+ − 197 if abndTable.funcIsNormalized():
+ − 198 fResult = abndTable.funcFilterAbundanceByPercentile(dPercentileCutOff = float(dPercentile), dPercentageAbovePercentile = float(dPercentage))
+ − 199 if fResult:
+ − 200 print "ManipulateTable::"+abndTable.funcGetName()+" has been reduced by percentile and now has "+str(len(abndTable.funcGetFeatureNames()))+" features."
+ − 201 else:
+ − 202 print "ManipulateTable::ERROR. "+abndTable.funcGetName()+" could not be reduced by percentile."
+ − 203 else:
+ − 204 print "ManipulateTable::"+abndTable.funcGetName()+" was NOT normalized and so the percentile filter is invalid, please indicate to normalize the table."
+ − 205
+ − 206 # Filter on abundance (should go after filter on percentile because the filter on percentile
+ − 207 # needs the full distribution of features in a sample
+ − 208 if args.strFilterAbundance:
+ − 209 dAbundance,iMinSamples = args.strFilterAbundance.split(",")
+ − 210 dAbundance = float(dAbundance)
+ − 211 iMinSamples = int(iMinSamples)
+ − 212 for abndTable in lsTables:
+ − 213 if abndTable.funcIsNormalized():
+ − 214 fResult = abndTable.funcFilterAbundanceByMinValue(dMinAbundance=dAbundance,iMinSamples=iMinSamples)
+ − 215 if fResult:
+ − 216 print "ManipulateTable::"+abndTable.funcGetName()+" has been reduced by minimum relative abundance value and now has "+str(len(abndTable.funcGetFeatureNames()))+" features."
+ − 217 else:
+ − 218 print "ManipulateTable::ERROR. "+abndTable.funcGetName()+" could not be reduced by percentile."
+ − 219 else:
+ − 220 print "ManipulateTable::"+abndTable.funcGetName()+" was NOT normalized and so the abundance filter is invalid, please indicate to normalize the table."
+ − 221
+ − 222 #if args.dCuttOff:
+ − 223 # print "Standard deviation filtering not completed"
+ − 224 # for abndTable in lsTables:
+ − 225 # abndTable.funcFilterFeatureBySD(dMinSDCuttOff=args.dCuttOff)
+ − 226 # if fResult:
+ − 227 # print "ManipulateTable::"+abndTable.funcGetName()+" has been reduced by standard deviation and now has "+str(len(abndTable.funcGetFeatureNames()))+" features."
+ − 228 # else:
+ − 229 # print "ManipulateTable::ERROR. "+abndTable.funcGetName()+" could not be reduced by standard devation."
+ − 230
+ − 231 # Need to normalize again after abundance data filtering given removing features breaks the normalization
+ − 232 # This happends twice because normalization is required to make the abundance data to filter on ;-)
+ − 233 # Normalize if needed
+ − 234 if args.fNormalize:
+ − 235 for abndTable in lsTables:
+ − 236 fResult = abndTable.funcNormalize()
+ − 237 if fResult:
+ − 238 print "ManipulateTable::"+abndTable.funcGetName()+" was normalized after filtering on abundance data."
+ − 239
+ − 240 #Manipulate lineage
+ − 241 if args.fPrefixClades:
+ − 242 for abndTable in lsTables:
+ − 243 fResult = abndTable.funcAddCladePrefixToFeatures()
+ − 244 if fResult:
+ − 245 print "ManipulateTable::Clade Prefix was added to "+abndTable.funcGetName()
+ − 246 else:
+ − 247 print "ManipulateTable::ERROR. Clade Prefix was NOT added to "+abndTable.funcGetName()
+ − 248
+ − 249 # Under development
+ − 250 # Reduce dimensionality
+ − 251 #if args.fDoPCA:
+ − 252 # pcaCur = PCA()
+ − 253 # for abndTable in lsTables:
+ − 254 #
+ − 255 # # Add data features
+ − 256 # # Make data components and add to abundance table
+ − 257 # pcaCur.loadData(abndTable,True)
+ − 258 # pcaCur.run(fASTransform=True)
+ − 259 # ldVariance = pcaCur.getVariance()
+ − 260 # lldComponents = pcaCur.getComponents()
+ − 261 # # Make Names
+ − 262 # lsNamesData = ["Data_PC"+str((tpleVariance[0]+1))+"_"+re.sub("[\.|-]","_",str(tpleVariance[1])) for tpleVariance in enumerate(ldVariance)]
+ − 263 # abndTable.funcAddDataFeature(lsNamesData,lldComponents)
+ − 264 #
+ − 265 # # Add metadata features
+ − 266 # # Convert metadata to an input for PCA
+ − 267 # pcaCur.loadData(pcaCur.convertMetadataForPCA(abndTable),False)
+ − 268 # fSuccessful = pcaCur.run(fASTransform=False)
+ − 269 # if(fSuccessful):
+ − 270 # ldVariance = pcaCur.getVariance()
+ − 271 # lldComponents = pcaCur.getComponents()
+ − 272 # # Make Names
+ − 273 # lsNamesMetadata = ["Metadata_PC"+str((tpleVariance[0]+1))+"_"+re.sub("[\.|-]","_",str(tpleVariance[1])) for tpleVariance in enumerate(ldVariance)]
+ − 274 # # Make metadata components and add to abundance
+ − 275 # llsMetadata = [list(npdRow) for npdRow in lldComponents]
+ − 276 # abndTable.funcAddMetadataFeature(lsNamesMetadata, llsMetadata)
+ − 277 # else:
+ − 278 # print "ManipulateTable::No metadata to PCA, no PCA components added to file based on metadata"
+ − 279
+ − 280 #Manipulate based on metadata
+ − 281 if args.strStratifyBy:
+ − 282 labndStratifiedTables = []
+ − 283 for abndTable in lsTables:
+ − 284 labndResult = abndTable.funcStratifyByMetadata(strMetadata=args.strStratifyBy)
+ − 285 print "ManipulateTable::"+abndTable.funcGetName()+" was stratified by "+args.strStratifyBy+" in to "+str(len(labndResult))+" tables."
+ − 286 labndStratifiedTables.extend(labndResult)
+ − 287 lsTables = labndStratifiedTables
+ − 288
+ − 289 if len(lsTables) == 1:
+ − 290 lsTables[0].funcWriteToFile(args.strOutFile)
+ − 291 else:
+ − 292 iIndex = 1
+ − 293 for abndManTable in lsTables:
+ − 294 abndManTable.funcWriteToFile(lsPieces[0]+str(iIndex)+lsPieces[1])
+ − 295 iIndex = iIndex + 1
