comparison MicroPITA.py @ 0:0de566f21448 draft default tip

v2
author sagun98
date Thu, 03 Jun 2021 18:13:32 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:0de566f21448
1 #!/usr/bin/env python
2 """
3 Author: Timothy Tickle
4 Description: Class to Run analysis for the microPITA paper
5 """
6
7 #####################################################################################
8 #Copyright (C) <2012>
9 #
10 #Permission is hereby granted, free of charge, to any person obtaining a copy of
11 #this software and associated documentation files (the "Software"), to deal in the
12 #Software without restriction, including without limitation the rights to use, copy,
13 #modify, merge, publish, distribute, sublicense, and/or sell copies of the Software,
14 #and to permit persons to whom the Software is furnished to do so, subject to
15 #the following conditions:
16 #
17 #The above copyright notice and this permission notice shall be included in all copies
18 #or substantial portions of the Software.
19 #
20 #THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
21 #INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
22 #PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
23 #HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
24 #OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
25 #SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
26 #####################################################################################
27
28 __author__ = "Timothy Tickle"
29 __copyright__ = "Copyright 2012"
30 __credits__ = ["Timothy Tickle"]
31 __license__ = "MIT"
32 __maintainer__ = "Timothy Tickle"
33 __email__ = "ttickle@sph.harvard.edu"
34 __status__ = "Development"
35
36 import sys
37 import argparse
38 from src.breadcrumbs.src.AbundanceTable import AbundanceTable
39 from src.breadcrumbs.src.ConstantsBreadCrumbs import ConstantsBreadCrumbs
40 from src.breadcrumbs.src.Metric import Metric
41 from src.breadcrumbs.src.KMedoids import Kmedoids
42 from src.breadcrumbs.src.MLPYDistanceAdaptor import MLPYDistanceAdaptor
43 from src.breadcrumbs.src.SVM import SVM
44 from src.breadcrumbs.src.UtilityMath import UtilityMath
45
46 from src.ConstantsMicropita import ConstantsMicropita
47 import csv
48 import logging
49 import math
50 import mlpy
51 import numpy as np
52 import operator
53 import os
54 import random
55 import scipy.cluster.hierarchy as hcluster
56 import scipy.spatial.distance
57 from types import *
58
59 class MicroPITA:
60 """
61 Selects samples from a first tier of a multi-tiered study to be used in a second tier.
62 Different methods can be used for selection.
63 The expected input is an abundance table (and potentially a text file of targeted features,
64 if using the targeted features option). Output is a list of samples exhibiting the
65 characteristics of interest.
66 """
67
68 #Constants
69 #Diversity metrics Alpha
70 c_strInverseSimpsonDiversity = Metric.c_strInvSimpsonDiversity
71 c_strChao1Diversity = Metric.c_strChao1Diversity
72
73 #Diversity metrics Beta
74 c_strBrayCurtisDissimilarity = Metric.c_strBrayCurtisDissimilarity
75
76 #Additive inverses of diversity metrics beta
77 c_strInvBrayCurtisDissimilarity = Metric.c_strInvBrayCurtisDissimilarity
78
79 #Technique Names
80 ConstantsMicropita.c_strDiversity2 = ConstantsMicropita.c_strDiversity+"_C"
81
82 #Targeted feature settings
83 c_strTargetedRanked = ConstantsMicropita.c_strTargetedRanked
84 c_strTargetedAbundance = ConstantsMicropita.c_strTargetedAbundance
85
86 #Technique groupings
87 # c_lsDiversityMethods = [ConstantsMicropita.c_strDiversity,ConstantsMicropita.c_strDiversity2]
88
89 #Converts ecology metrics into standardized method selection names
90 dictConvertAMetricDiversity = {c_strInverseSimpsonDiversity:ConstantsMicropita.c_strDiversity, c_strChao1Diversity:ConstantsMicropita.c_strDiversity2}
91 # dictConvertMicroPITAToAMetric = {ConstantsMicropita.c_strDiversity:c_strInverseSimpsonDiversity, ConstantsMicropita.c_strDiversity2:c_strChao1Diversity}
92 dictConvertBMetricToMethod = {c_strBrayCurtisDissimilarity:ConstantsMicropita.c_strRepresentative}
93 dictConvertInvBMetricToMethod = {c_strBrayCurtisDissimilarity:ConstantsMicropita.c_strExtreme}
94
95 #Linkage used in the Hierarchical clustering
96 c_strHierarchicalClusterMethod = 'average'
97
98 ####Group 1## Diversity
99 #Testing: Happy path Testing (8)
100 def funcGetTopRankedSamples(self, lldMatrix = None, lsSampleNames = None, iTopAmount = None):
101 """
102 Given a list of lists of measurements, for each list the indices of the highest values are returned. If lsSamplesNames is given
103 it is treated as a list of string names that is in the order of the measurements in each list. Indices are returned or the sample
104 names associated with the indices.
105
106 :param lldMatrix: List of lists [[value,value,value,value],[value,value,value,value]].
107 :type: List of lists List of measurements. Each list is a different measurement. Each measurement in positionally related to a sample.
108 :param lsSampleNames: List of sample names positionally related (the same) to each list (Optional).
109 :type: List of strings List of strings.
110 :param iTopAmount: The amount of top measured samples (assumes the higher measurements are better).
111 :type: integer Integer amount of sample names/ indices to return.
112 :return List: List of samples to be selected.
113 """
114 topRankListRet = []
115 for rowMetrics in lldMatrix:
116 #Create 2 d array to hold value and index and sort
117 liIndexX = [rowMetrics,range(len(rowMetrics))]
118 liIndexX[1].sort(key = liIndexX[0].__getitem__,reverse = True)
119
120 if lsSampleNames:
121 topRankListRet.append([lsSampleNames[iIndex] for iIndex in liIndexX[1][:iTopAmount]])
122 else:
123 topRankListRet.append(liIndexX[1][:iTopAmount])
124
125 return topRankListRet
126
127 ####Group 2## Representative Dissimilarity
128 #Testing: Happy path tested 1
129 def funcGetCentralSamplesByKMedoids(self, npaMatrix=None, sMetric=None, lsSampleNames=None, iNumberSamplesReturned=0, istmBetaMatrix=None, istrmTree=None, istrmEnvr=None):
130 """
131 Gets centroid samples by k-medoids clustering of a given matrix.
132
133 :param npaMatrix: Numpy array where row=features and columns=samples
134 :type: Numpy array Abundance Data.
135 :param sMetric: String name of beta metric used as the distance metric.
136 :type: String String name of beta metric.
137 :param lsSampleNames: The names of the sample
138 :type: List List of strings
139 :param iNumberSamplesReturned: Number of samples to return, each will be a centroid of a sample.
140 :type: Integer Number of samples to return
141 :return List: List of selected samples.
142 :param istmBetaMatrix: File with beta-diversity matrix
143 :type: File stream or file path string
144 """
145
146 #Count of how many rows
147 sampleCount = npaMatrix.shape[0]
148 if iNumberSamplesReturned > sampleCount:
149 logging.error("MicroPITA.funcGetCentralSamplesByKMedoids:: There are not enough samples to return the amount of samples specified. Return sample count = "+str(iNumberSamplesReturned)+". Sample number = "+str(sampleCount)+".")
150 return False
151
152 #If the cluster count is equal to the sample count return all samples
153 if sampleCount == iNumberSamplesReturned:
154 return list(lsSampleNames)
155
156 #Get distance matrix
157 distanceMatrix=scipy.spatial.distance.squareform(Metric.funcReadMatrixFile(istmMatrixFile=istmBetaMatrix,lsSampleOrder=lsSampleNames)[0]) if istmBetaMatrix else Metric.funcGetBetaMetric(npadAbundancies=npaMatrix, sMetric=sMetric, istrmTree=istrmTree, istrmEnvr=istrmEnvr, lsSampleOrder=lsSampleNames)
158 if type(distanceMatrix) is BooleanType:
159 logging.error("MicroPITA.funcGetCentralSamplesByKMedoids:: Could not read in the supplied distance matrix, returning false.")
160 return False
161
162 # Handle unifrac output
163 if sMetric in [Metric.c_strUnifracUnweighted,Metric.c_strUnifracWeighted]:
164 distanceMatrix = distanceMatrix[0]
165
166 #Log distance matrix
167 logging.debug("MicroPITA.funcGetCentralSamplesByKMedoids:: Distance matrix for representative selection using metric="+str(sMetric))
168
169 distance = MLPYDistanceAdaptor(npaDistanceMatrix=distanceMatrix, fIsCondensedMatrix=True)
170
171 #Create object to determine clusters/medoids
172 medoidsMaker = Kmedoids(k=iNumberSamplesReturned, dist=distance)
173 #medoidsData includes(1d numpy array, medoids indexes;
174 # 1d numpy array, non-medoids indexes;
175 # 1d numpy array, cluster membership for non-medoids;
176 # double, cost of configuration)
177 #npaMatrix is samples x rows
178 #Build a matrix of lists of indicies to pass to the distance matrix
179 lliIndicesMatrix = [[iIndexPosition] for iIndexPosition in xrange(0,len(npaMatrix))]
180 medoidsData = medoidsMaker.compute(np.array(lliIndicesMatrix))
181 logging.debug("MicroPITA.funcGetCentralSamplesByKMedoids:: Results from the kmedoid method in representative selection:")
182 logging.debug(str(medoidsData))
183
184 #If returning the same amount of clusters and samples
185 #Return centroids
186 selectedIndexes = medoidsData[0]
187 return [lsSampleNames[selectedIndexes[index]] for index in xrange(0,iNumberSamplesReturned)]
188
189 ####Group 3## Highest Dissimilarity
190 #Testing: Happy path tested
191 def funcSelectExtremeSamplesFromHClust(self, strBetaMetric, npaAbundanceMatrix, lsSampleNames, iSelectSampleCount, istmBetaMatrix=None, istrmTree=None, istrmEnvr=None):
192 """
193 Select extreme samples from HClustering.
194
195 :param strBetaMetric: The beta metric to use for distance matrix generation.
196 :type: String The name of the beta metric to use.
197 :param npaAbundanceMatrix: Numpy array where row=samples and columns=features.
198 :type: Numpy Array Abundance data.
199 :param lsSampleNames: The names of the sample.
200 :type: List List of strings.
201 :param iSelectSampleCount: Number of samples to select (return).
202 :type: Integer Integer number of samples returned.
203 :return Samples: List of samples.
204 :param istmBetaMatrix: File with beta-diversity matrix
205 :type: File stream or file path string
206 """
207
208 #If they want all the sample count, return all sample names
209 iSampleCount=len(npaAbundanceMatrix[:,0])
210 if iSelectSampleCount==iSampleCount:
211 return lsSampleNames
212
213 #Holds the samples to be returned
214 lsReturnSamplesRet = []
215
216 #Generate beta matrix
217 #Returns condensed matrix
218 tempDistanceMatrix = scipy.spatial.distance.squareform(Metric.funcReadMatrixFile(istmMatrixFile=istmBetaMatrix,lsSampleOrder=lsSampleNames)[0]) if istmBetaMatrix else Metric.funcGetBetaMetric(npadAbundancies=npaAbundanceMatrix, sMetric=strBetaMetric, istrmTree=istrmTree, istrmEnvr=istrmEnvr, lsSampleOrder=lsSampleNames, fAdditiveInverse = True)
219
220 if strBetaMetric in [Metric.c_strUnifracUnweighted,Metric.c_strUnifracWeighted]:
221 tempDistanceMatrix = tempDistanceMatrix[0]
222
223 if type(tempDistanceMatrix) is BooleanType:
224 logging.error("MicroPITA.funcSelectExtremeSamplesFromHClust:: Could not read in the supplied distance matrix, returning false.")
225 return False
226
227 if istmBetaMatrix:
228 tempDistanceMatrix = 1-tempDistanceMatrix
229
230 #Feed beta matrix to linkage to cluster
231 #Send condensed matrix
232 linkageMatrix = hcluster.linkage(tempDistanceMatrix, method=self.c_strHierarchicalClusterMethod)
233
234 #Extract cluster information from dendrogram
235 #The linakge matrix is of the form
236 #[[int1 int2 doube int3],...]
237 #int1 and int1 are the paired samples indexed at 0 and up.
238 #each list is an entry for a branch that is number starting with the first
239 #list being sample count index + 1
240 #each list is then named by an increment as they appear
241 #this means that if a number is in the list and is = sample count or greater it is not
242 #terminal and is instead a branch.
243 #This method just takes the lowest metric measurement (highest distance pairs/clusters)
244 #Works much better than the original technique
245 #get total number of samples
246
247 iCurrentSelectCount = 0
248 for row in linkageMatrix:
249 #Get nodes ofthe lowest pairing (so the furthest apart pair)
250 iNode1 = int(row[0])
251 iNode2 = int(row[1])
252 #Make sure the nodes are a terminal node (sample) and not a branch in the dendrogram
253 #The branching in the dendrogram will start at the number of samples and increment higher.
254 #Add each of the pair one at a time breaking when enough samples are selected.
255 if iNode1<iSampleCount:
256 lsReturnSamplesRet.append(lsSampleNames[iNode1])
257 iCurrentSelectCount = iCurrentSelectCount + 1
258 if iCurrentSelectCount == iSelectSampleCount:
259 break
260 if iNode2<iSampleCount:
261 lsReturnSamplesRet.append(lsSampleNames[iNode2])
262 iCurrentSelectCount = iCurrentSelectCount + 1
263 if iCurrentSelectCount == iSelectSampleCount:
264 break
265
266 #Return selected samples
267 return lsReturnSamplesRet
268
269 ####Group 4## Rank Average of user Defined Taxa
270 #Testing: Happy Path Tested
271 def funcGetAverageAbundanceSamples(self, abndTable, lsTargetedFeature, fRank=False):
272 """
273 Averages feature abundance or ranked abundance. Expects a column 0 of taxa id that is skipped.
274
275 :param abndTable: Abundance Table to analyse
276 :type: AbundanceTable Abundance Table
277 :param lsTargetedFeature: String names
278 :type: list list of string names of features (bugs) which are measured after ranking against the full sample
279 :param fRank: Indicates to rank the abundance before getting the average abundance of the features (default false)
280 :type: boolean Flag indicating ranking abundance before calculating average feature measurement (false= no ranking)
281 :return List of lists or boolean: List of lists or False on error. One internal list per sample indicating the sample,
282 feature average abundance or ranked abundance. Lists will already be sorted.
283 For not Ranked [[sample,average abundance of selected feature,1]]
284 For Ranked [[sample,average ranked abundance, average abundance of selected feature]]
285 Error Returns false
286 """
287
288 llAbundance = abndTable.funcGetAverageAbundancePerSample(lsTargetedFeature)
289 if not llAbundance:
290 logging.error("MicroPITA.funcGetAverageAbundanceSamples:: Could not get average abundance, returned false. Make sure the features (bugs) are spelled correctly and in the abundance table.")
291 return False
292 #Add a space for ranking if needed
293 #Not ranked will be [[sSample,average abundance,1]]
294 #(where 1 will not discriminant ties if used in later functions, so this generalizes)
295 #Ranked will be [[sSample, average rank, average abundance]]
296 llRetAbundance = [[llist[0],-1,llist[1]] for llist in llAbundance]
297 #Rank if needed
298 if fRank:
299 abndRanked = abndTable.funcRankAbundance()
300 if abndRanked == None:
301 logging.error("MicroPITA.funcGetAverageAbundanceSamples:: Could not rank the abundance table, returned false.")
302 return False
303 llRetRank = abndRanked.funcGetAverageAbundancePerSample(lsTargetedFeature)
304 if not llRetRank:
305 logging.error("MicroPITA.funcGetAverageAbundanceSamples:: Could not get average ranked abundance, returned false. Make sure the features (bugs) are spelled correctly and in the abundance table.")
306 return False
307 dictRanks = dict(llRetRank)
308 llRetAbundance = [[a[0],dictRanks[a[0]],a[2]] for a in llRetAbundance]
309
310 #Sort first for ties and then for the main feature
311 if not fRank or ConstantsMicropita.c_fBreakRankTiesByDiversity:
312 llRetAbundance = sorted(llRetAbundance, key = lambda sampleData: sampleData[2], reverse = not fRank)
313 if fRank:
314 llRetAbundance = sorted(llRetAbundance, key = lambda sampleData: sampleData[1], reverse = not fRank)
315 return llRetAbundance
316
317 #Testing: Happy Path Tested
318 def funcSelectTargetedTaxaSamples(self, abndMatrix, lsTargetedTaxa, iSampleSelectionCount, sMethod = ConstantsMicropita.lsTargetedFeatureMethodValues[0]):
319 """
320 Selects samples with the highest ranks or abundance of targeted features.
321 If ranked, select the highest abundance for tie breaking
322
323 :param abndMatrix: Abundance table to analyse
324 :type: AbundanceTable Abundance table
325 :param lsTargetedTaxa: List of features
326 :type: list list of strings
327 :param iSampleSelectionCount: Number of samples to select
328 :type: integer integer
329 :param sMethod: Method to select targeted features
330 :type: string String (Can be values found in ConstantsMicropita.lsTargetedFeatureMethodValues)
331 :return List of strings: List of sample names which were selected
332 List of strings Empty list is returned on an error.
333 """
334
335 #Check data
336 if(len(lsTargetedTaxa) < 1):
337 logging.error("MicroPITA.funcSelectTargetedTaxaSamples. Taxa defined selection was requested but no features were given.")
338 return []
339
340 lsTargetedSamples = self.funcGetAverageAbundanceSamples(abndTable=abndMatrix, lsTargetedFeature=lsTargetedTaxa,
341 fRank=sMethod.lower() == self.c_strTargetedRanked.lower())
342 #If an error occured or the key word for the method was not recognized
343 if lsTargetedSamples == False:
344 logging.error("MicroPITA.funcSelectTargetedTaxaSamples:: Was not able to select for the features given. So targeted feature selection was performed. Check to make sure the features are spelled correctly and exist in the abundance file.")
345 return []
346
347 #Select from results
348 return [sSample[0] for sSample in lsTargetedSamples[:iSampleSelectionCount]]
349
350 ####Group 5## Random
351 #Testing: Happy path Tested
352 def funcGetRandomSamples(self, lsSamples=None, iNumberOfSamplesToReturn=0):
353 """
354 Returns random sample names of the number given. No replacement.
355
356 :param lsSamples: List of sample names
357 :type: list list of strings
358 :param iNumberOfSamplesToReturn: Number of samples to select
359 :type: integer integer.
360 :return List: List of selected samples (strings).
361 """
362
363 #Input matrix sample count
364 sampleCount = len(lsSamples)
365
366 #Return the full matrix if they ask for a return matrix where length == original
367 if(iNumberOfSamplesToReturn >= sampleCount):
368 return lsSamples
369
370 #Get the random indices for the sample (without replacement)
371 liRandomIndices = random.sample(range(sampleCount), iNumberOfSamplesToReturn)
372
373 #Create a boolean array of if indexes are to be included in the reduced array
374 return [sSample for iIndex, sSample in enumerate(lsSamples) if iIndex in liRandomIndices]
375
376 #Happy path tested (case 3)
377 def funcGetAveragePopulation(self, abndTable, lfCompress):
378 """
379 Get the average row per column in the abndtable.
380
381 :param abndTable: AbundanceTable of data to be averaged
382 :type: AbudanceTable
383 :param lfCompress: List of boolean flags (false means to remove sample before averaging
384 :type: List of floats
385 :return List of doubles:
386 """
387 if sum(lfCompress) == 0:
388 return []
389
390 #Get the average populations
391 lAverageRet = []
392
393 for sFeature in abndTable.funcGetAbundanceCopy():
394 sFeature = list(sFeature)[1:]
395 sFeature=np.compress(lfCompress,sFeature,axis=0)
396 lAverageRet.append(sum(sFeature)/float(len(sFeature)))
397 return lAverageRet
398
399 #Happy path tested (2 cases)
400 def funcGetDistanceFromAverage(self, abndTable,ldAverage,lsSamples,lfSelected):
401 """
402 Given an abundance table and an average sample, this returns the distance of each sample
403 (measured using brays-curtis dissimilarity) from the average.
404 The distances are reduced by needing to be in the lsSamples and being a true in the lfSelected
405 (which is associated with the samples in the order of the samples in the abundance table;
406 use abundancetable.funcGetSampleNames() to see the order if needed).
407
408 :param abndTable: Abundance table holding the data to be analyzed.
409 :type: AbundanceTable
410 :param ldAverage: Average population (Average features of the abundance table of samples)
411 :type: List of doubles which represent the average population
412 :param lsSamples: These are the only samples used in the analysis
413 :type: List of strings (sample ids)
414 :param lfSelected: Samples to be included in the analysis
415 :type: List of boolean (true means include)
416 :return: List of distances (doubles)
417 """
418 #Get the distance from label 1 of all samples in label0 splitting into selected and not selected lists
419 ldSelectedDistances = []
420
421 for sSampleName in [sSample for iindex, sSample in enumerate(lsSamples) if lfSelected[iindex]]:
422 #Get the sample measurements
423 ldSelectedDistances.append(Metric.funcGetBrayCurtisDissimilarity(np.array([abndTable.funcGetSample(sSampleName),ldAverage]))[0])
424 return ldSelectedDistances
425
426 #Happy path tested (1 case)
427 def funcMeasureDistanceFromLabelToAverageOtherLabel(self, abndTable, lfGroupOfInterest, lfGroupOther):
428 """
429 Get the distance of samples from one label from the average sample of not the label.
430 Note: This assumes 2 classes.
431
432 :param abndTable: Table of data to work out of.
433 :type: Abundace Table
434 :param lfGroupOfInterest: Boolean indicator of the sample being in the first group.
435 :type: List of floats, true indicating an individual in the group of interest.
436 :param lfGroupOther: Boolean indicator of the sample being in the other group.
437 :type: List of floats, true indicating an individual in the
438 :return List of List of doubles: [list of tuples (string sample name,double distance) for the selected population, list of tuples for the not selected population]
439 """
440 #Get all sample names
441 lsAllSamples = abndTable.funcGetSampleNames()
442
443 #Get average populations
444 lAverageOther = self.funcGetAveragePopulation(abndTable=abndTable, lfCompress=lfGroupOther)
445
446 #Get the distance from the average of the other label (label 1)
447 ldSelectedDistances = self.funcGetDistanceFromAverage(abndTable=abndTable, ldAverage=lAverageOther,
448 lsSamples=lsAllSamples, lfSelected=lfGroupOfInterest)
449
450 return zip([lsAllSamples[iindex] for iindex, fGroup in enumerate(lfGroupOfInterest) if fGroup],ldSelectedDistances)
451
452 #Happy path tested (1 test case)
453 def funcPerformDistanceSelection(self, abndTable, iSelectionCount, sLabel, sValueOfInterest):
454 """
455 Given metadata, metadata of one value (sValueOfInterest) is measured from the average (centroid) value of another label group.
456 An iSelectionCount of samples is selected from the group of interest closest to and furthest from the centroid of the other group.
457
458 :params abndTable: Abundance of measurements
459 :type: AbundanceTable
460 :params iSelectionCount: The number of samples selected per sample.
461 :type: Integer Integer greater than 0
462 :params sLabel: ID of the metadata which is the supervised label
463 :type: String
464 :params sValueOfInterest: Metadata value in the sLabel metadta row of the abundance table which defines the group of interest.
465 :type: String found in the abundance table metadata row indicated by sLabel.
466 :return list list of tuples (samplename, distance) [[iSelectionCount of tuples closest to the other centroid], [iSelectionCount of tuples farthest from the other centroid], [all tuples of samples not selected]]
467 """
468
469 lsMetadata = abndTable.funcGetMetadata(sLabel)
470 #Other metadata values
471 lsUniqueOtherValues = list(set(lsMetadata)-set(sValueOfInterest))
472
473 #Get boolean indicator of values of interest
474 lfLabelsInterested = [sValueOfInterest == sValue for sValue in lsMetadata]
475
476 #Get the distances of the items of interest from the other metadata values
477 dictDistanceAverages = {}
478 for sOtherLabel in lsUniqueOtherValues:
479 #Get boolean indicator of labels not of interest
480 lfLabelsOther = [sOtherLabel == sValue for sValue in lsMetadata]
481
482 #Get the distances of data from two different groups to the average of the other
483 ldValueDistances = dict(self.funcMeasureDistanceFromLabelToAverageOtherLabel(abndTable, lfLabelsInterested, lfLabelsOther))
484
485 for sKey in ldValueDistances:
486 dictDistanceAverages[sKey] = ldValueDistances[sKey] + dictDistanceAverages[sKey] if sKey in dictDistanceAverages else ldValueDistances[sKey]
487
488 #Finish average by dividing by length of lsUniqueOtherValues
489 ltpleAverageDistances = [(sKey, dictDistanceAverages[sKey]/float(len(lsUniqueOtherValues))) for sKey in dictDistanceAverages]
490
491 #Sort to extract extremes
492 ltpleAverageDistances = sorted(ltpleAverageDistances,key=operator.itemgetter(1))
493
494 #Get the closest and farthest distances
495 ltupleDiscriminantSamples = ltpleAverageDistances[:iSelectionCount]
496 ltupleDistinctSamples = ltpleAverageDistances[iSelectionCount*-1:]
497
498 #Remove the selected samples from the larger population of distances (better visualization)
499 ldSelected = [tpleSelected[0] for tpleSelected in ltupleDiscriminantSamples+ltupleDistinctSamples]
500
501 #Return discriminant tuples, distinct tuples, other tuples
502 return [ltupleDiscriminantSamples, ltupleDistinctSamples,
503 [tplData for tplData in ltpleAverageDistances if tplData[0] not in ldSelected]]
504
505 #Run the supervised method surrounding distance from centroids
506 #Happy path tested (3 test cases)
507 def funcRunSupervisedDistancesFromCentroids(self, abundanceTable, fRunDistinct, fRunDiscriminant,
508 xOutputSupFile, xPredictSupFile, strSupervisedMetadata,
509 iSampleSupSelectionCount, lsOriginalSampleNames, lsOriginalLabels, fAppendFiles = False):
510 """
511 Runs supervised methods based on measuring distances of one label from the centroid of another. NAs are evaluated as theirown group.
512
513 :param abundanceTable: AbundanceTable
514 :type: AbudanceTable Data to analyze
515 :param fRunDistinct: Run distinct selection method
516 :type: Boolean boolean (true runs method)
517 :param fRunDiscriminant: Run discriminant method
518 :type: Boolean boolean (true runs method)
519 :param xOutputSupFile: File output from supervised methods detailing data going into the method.
520 :type: String or FileStream
521 :param xPredictSupFile: File output from supervised methods distance results from supervised methods.
522 :type: String or FileStream
523 :param strSupervisedMetadata: The metadata that will be used to group samples.
524 :type: String
525 :param iSampleSupSelectionCount: Number of samples to select
526 :type: Integer int sample selection count
527 :param lsOriginalSampleNames: List of the sample names, order is important and should be preserved from the abundanceTable.
528 :type: List of samples
529 :param fAppendFiles: Indicates that output files already exist and appending is occuring.
530 :type: Boolean
531 :return Selected Samples: A dictionary of selected samples by selection ID
532 Dictionary {"Selection Method":["SampleID","SampleID"...]}
533 """
534 #Get labels and run one label against many
535 lstrMetadata = abundanceTable.funcGetMetadata(strSupervisedMetadata)
536 dictlltpleDistanceMeasurements = {}
537 for sMetadataValue in set(lstrMetadata):
538
539 #For now perform the selection here for the label of interest against the other labels
540 dictlltpleDistanceMeasurements.setdefault(sMetadataValue,[]).extend(self.funcPerformDistanceSelection(abndTable=abundanceTable,
541 iSelectionCount=iSampleSupSelectionCount, sLabel=strSupervisedMetadata, sValueOfInterest=sMetadataValue))
542
543 #Make expected output files for supervised methods
544 #1. Output file which is similar to an input file for SVMs
545 #2. Output file that is similar to the probabilitic output of a SVM (LibSVM)
546 #Manly for making output of supervised methods (Distance from Centroid) similar
547 #MicropitaVis needs some of these files
548 if xOutputSupFile:
549 if fAppendFiles:
550 SVM.funcUpdateSVMFileWithAbundanceTable(abndAbundanceTable=abundanceTable, xOutputSVMFile=xOutputSupFile,
551 lsOriginalLabels=lsOriginalLabels, lsSampleOrdering=lsOriginalSampleNames)
552 else:
553 SVM.funcConvertAbundanceTableToSVMFile(abndAbundanceTable=abundanceTable, xOutputSVMFile=xOutputSupFile,
554 sMetadataLabel=strSupervisedMetadata, lsOriginalLabels=lsOriginalLabels, lsSampleOrdering=lsOriginalSampleNames)
555
556 #Will contain the samples selected to return
557 #One or more of the methods may be active so this is why I am extending instead of just returning the result of each method type
558 dictSelectedSamplesRet = dict()
559 for sKey, ltplDistances in dictlltpleDistanceMeasurements.items():
560 if fRunDistinct:
561 dictSelectedSamplesRet.setdefault(ConstantsMicropita.c_strDistinct,[]).extend([ltple[0] for ltple in ltplDistances[1]])
562 if fRunDiscriminant:
563 dictSelectedSamplesRet.setdefault(ConstantsMicropita.c_strDiscriminant,[]).extend([ltple[0] for ltple in ltplDistances[0]])
564
565 if xPredictSupFile:
566 dictFlattenedDistances = dict()
567 [dictFlattenedDistances.setdefault(sKey, []).append(tple)
568 for sKey, lltple in dictlltpleDistanceMeasurements.items()
569 for ltple in lltple for tple in ltple]
570 if fAppendFiles:
571 self._updatePredictFile(xPredictSupFile=xPredictSupFile, xInputLabelsFile=xOutputSupFile,
572 dictltpleDistanceMeasurements=dictFlattenedDistances, abundanceTable=abundanceTable, lsOriginalSampleNames=lsOriginalSampleNames)
573 else:
574 self._writeToPredictFile(xPredictSupFile=xPredictSupFile, xInputLabelsFile=xOutputSupFile,
575 dictltpleDistanceMeasurements=dictFlattenedDistances, abundanceTable=abundanceTable, lsOriginalSampleNames=lsOriginalSampleNames)
576 return dictSelectedSamplesRet
577
578 #Two happy path test cases
579 def _updatePredictFile(self, xPredictSupFile, xInputLabelsFile, dictltpleDistanceMeasurements, abundanceTable, lsOriginalSampleNames):
580 """
581 Manages updating the predict file.
582
583 :param xPredictSupFile: File that has predictions (distances) from the supervised method.
584 :type: FileStream or String file path
585 :param xInputLabelsFile: File that as input to the supervised methods.
586 :type: FileStream or String file path
587 :param dictltpleDistanceMeasurements:
588 :type: Dictionary of lists of tuples {"labelgroup":[("SampleName",dDistance)...], "labelgroup":[("SampleName",dDistance)...]}
589 """
590
591 if not isinstance(xPredictSupFile, str):
592 xPredictSupFile.close()
593 xPredictSupFile = xPredictSupFile.name
594 csvr = open(xPredictSupFile,'r')
595
596 f = csv.reader(csvr,delimiter=ConstantsBreadCrumbs.c_strBreadCrumbsSVMSpace)
597 lsHeader = f.next()[1:]
598 dictlltpleRead = dict([(sHeader,[]) for sHeader in lsHeader])
599
600 #Read data in
601 iSampleIndex = 0
602 for sRow in f:
603 sLabel = sRow[0]
604 [dictlltpleRead[lsHeader[iDistanceIndex]].append((lsOriginalSampleNames[iSampleIndex],dDistance)) for iDistanceIndex, dDistance in enumerate(sRow[1:])
605 if not dDistance == ConstantsMicropita.c_sEmptyPredictFileValue]
606 iSampleIndex += 1
607
608 #Combine dictltpleDistanceMeasurements with new data
609 #If they share a key then merge keeping parameter data
610 #If they do not share the key, keep the full data
611 dictNew = {}
612 for sKey in dictltpleDistanceMeasurements.keys():
613 lsSamples = [tple[0] for tple in dictltpleDistanceMeasurements[sKey]]
614 dictNew[sKey] = dictltpleDistanceMeasurements[sKey]+[tple for tple in dictlltpleRead[sKey] if tple[0] not in lsSamples] if sKey in dictlltpleRead.keys() else dictltpleDistanceMeasurements[sKey]
615 for sKey in dictlltpleRead:
616 if sKey not in dictltpleDistanceMeasurements.keys():
617 dictNew[sKey] = dictlltpleRead[sKey]
618
619 #Call writer
620 self._writeToPredictFile(xPredictSupFile=xPredictSupFile, xInputLabelsFile=xInputLabelsFile,
621 dictltpleDistanceMeasurements=dictNew, abundanceTable=abundanceTable,
622 lsOriginalSampleNames=lsOriginalSampleNames, fFromUpdate=True)
623
624 #2 happy path test cases
625 def _writeToPredictFile(self, xPredictSupFile, xInputLabelsFile, dictltpleDistanceMeasurements, abundanceTable, lsOriginalSampleNames, fFromUpdate=False):
626 """
627 Write to the predict file.
628
629 :param xPredictSupFile: File that has predictions (distances) from the supervised method.
630 :type: FileStream or String file path
631 :param xInputLabelsFile: File that as input to the supervised methods.
632 :type: FileStream or String file path
633 :param dictltpleDistanceMeasurements:
634 :type: Dictionary of lists of tuples {"labelgroup":[("SampleName",dDistance)...], "labelgroup":[("SampleName",dDistance)...]}
635 :param abundanceTable: An abundance table of the sample data.
636 :type: AbundanceTable
637 :param lsOriginalSampleNames: Used if the file is being updated as the sample names so that it may be passed in and consistent with other writing.
638 Otherwise will use the sample names from the abundance table.
639 :type: List of strings
640 :param fFromUpdate: Indicates if this is part of an update to the file or not.
641 :type: Boolean
642 """
643
644 xInputLabelsFileName = xInputLabelsFile
645 if not isinstance(xInputLabelsFile,str):
646 xInputLabelsFileName = xInputLabelsFile.name
647 f = csv.writer(open(xPredictSupFile,"w") if isinstance(xPredictSupFile, str) else xPredictSupFile,delimiter=ConstantsBreadCrumbs.c_strBreadCrumbsSVMSpace)
648
649 lsAllSampleNames = abundanceTable.funcGetSampleNames()
650 lsLabels = SVM.funcReadLabelsFromFile(xSVMFile=xInputLabelsFileName, lsAllSampleNames= lsOriginalSampleNames if fFromUpdate else lsAllSampleNames,
651 isPredictFile=False)
652 dictLabels = dict([(sSample,sLabel) for sLabel in lsLabels.keys() for sSample in lsLabels[sLabel]])
653
654 #Dictionay keys will be used to order the predict file
655 lsMeasurementKeys = dictltpleDistanceMeasurements.keys()
656 #Make header
657 f.writerow(["labels"]+lsMeasurementKeys)
658
659 #Reformat dictionary to make it easier to use
660 for sKey in dictltpleDistanceMeasurements:
661 dictltpleDistanceMeasurements[sKey] = dict([ltpl for ltpl in dictltpleDistanceMeasurements[sKey]])
662
663 for sSample in lsOriginalSampleNames:
664 #Make body of file
665 f.writerow([dictLabels.get(sSample,ConstantsMicropita.c_sEmptyPredictFileValue)]+
666 [str(dictltpleDistanceMeasurements[sKey].get(sSample,ConstantsMicropita.c_sEmptyPredictFileValue))
667 for sKey in lsMeasurementKeys])
668
669 def _funcRunNormalizeSensitiveMethods(self, abndData, iSampleSelectionCount, dictSelectedSamples, lsAlphaMetrics, lsBetaMetrics, lsInverseBetaMetrics,
670 fRunDiversity, fRunRepresentative, fRunExtreme, strAlphaMetadata=None,
671 istmBetaMatrix=None, istrmTree=None, istrmEnvr=None, fInvertDiversity=False):
672 """
673 Manages running methods that are sensitive to normalization. This is called twice, once for the set of methods which should not be normalized and the other
674 for the set that should be normalized.
675
676 :param abndData: Abundance table object holding the samples to be measured.
677 :type: AbundanceTable
678 :param iSampleSelectionCount The number of samples to select per method.
679 :type: Integer
680 :param dictSelectedSamples Will be added to as samples are selected {"Method:["strSelectedSampleID","strSelectedSampleID"...]}.
681 :type: Dictionary
682 :param lsAlphaMetrics: List of alpha metrics to use on alpha metric dependent assays (like highest diversity).
683 :type: List of strings
684 :param lsBetaMetrics: List of beta metrics to use on beta metric dependent assays (like most representative).
685 :type: List of strings
686 :param lsInverseBetaMetrics: List of inverse beta metrics to use on inverse beta metric dependent assays (like most dissimilar).
687 :type: List of strings
688 :param fRunDiversity: Run Diversity based methods (true indicates run).
689 :type: Boolean
690 :param fRunRepresentative: Run Representative based methods (true indicates run).
691 :type: Boolean
692 :param fRunExtreme: Run Extreme based methods (true indicates run).
693 :type: Boolean
694 :param istmBetaMatrix: File that has a precalculated beta matrix
695 :type: File stream or File path string
696 :return Selected Samples: Samples selected by methods.
697 Dictionary {"Selection Method":["SampleID","SampleID","SampleID",...]}
698 """
699
700 #Sample ids/names
701 lsSampleNames = abndData.funcGetSampleNames()
702
703 #Generate alpha metrics and get most diverse
704 if fRunDiversity:
705
706 #Get Alpha metrics matrix
707 internalAlphaMatrix = None
708 #Name of technique
709 strMethod = [strAlphaMetadata] if strAlphaMetadata else lsAlphaMetrics
710
711 #If given an alpha-diversity metadata
712 if strAlphaMetadata:
713 internalAlphaMatrix = [[float(strNum) for strNum in abndData.funcGetMetadata(strAlphaMetadata)]]
714 else:
715 #Expects Observations (Taxa (row) x sample (column))
716 #Returns [[metric1-sample1, metric1-sample2, metric1-sample3],[metric1-sample1, metric1-sample2, metric1-sample3]]
717 internalAlphaMatrix = Metric.funcBuildAlphaMetricsMatrix(npaSampleAbundance = abndData.funcGetAbundanceCopy()
718 if not abndData.funcIsSummed()
719 else abndData.funcGetFeatureAbundanceTable(abndData.funcGetTerminalNodes()).funcGetAbundanceCopy(),
720 lsSampleNames = lsSampleNames, lsDiversityMetricAlpha = lsAlphaMetrics)
721
722 if internalAlphaMatrix:
723 #Invert measurments
724 if fInvertDiversity:
725 lldNewDiversity = []
726 for lsLine in internalAlphaMatrix:
727 lldNewDiversity.append([1/max(dValue,ConstantsMicropita.c_smallNumber) for dValue in lsLine])
728 internalAlphaMatrix = lldNewDiversity
729 #Get top ranked alpha diversity by most diverse
730 #Expects [[sample1,sample2,sample3...],[sample1,sample2,sample3..],...]
731 #Returns [[sampleName1, sampleName2, sampleNameN],[sampleName1, sampleName2, sampleNameN]]
732 mostDiverseAlphaSamplesIndexes = self.funcGetTopRankedSamples(lldMatrix=internalAlphaMatrix, lsSampleNames=lsSampleNames, iTopAmount=iSampleSelectionCount)
733
734 #Add to results
735 for index in xrange(0,len(strMethod)):
736 strSelectionMethod = self.dictConvertAMetricDiversity.get(strMethod[index],ConstantsMicropita.c_strDiversity+"="+strMethod[index])
737 dictSelectedSamples.setdefault(strSelectionMethod,[]).extend(mostDiverseAlphaSamplesIndexes[index])
738
739 logging.info("MicroPITA.funcRunNormalizeSensitiveMethods:: Selected Samples 1b")
740 logging.info(dictSelectedSamples)
741
742 #Generate beta metrics and
743 if fRunRepresentative or fRunExtreme:
744
745 #Abundance matrix transposed
746 npaTransposedAbundance = UtilityMath.funcTransposeDataMatrix(abndData.funcGetAbundanceCopy(), fRemoveAdornments=True)
747
748 #Get center selection using clusters/tiling
749 #This will be for beta metrics in normalized space
750 if fRunRepresentative:
751
752 if istmBetaMatrix:
753 #Get representative dissimilarity samples
754 medoidSamples=self.funcGetCentralSamplesByKMedoids(npaMatrix=npaTransposedAbundance, sMetric=ConstantsMicropita.c_custom, lsSampleNames=lsSampleNames, iNumberSamplesReturned=iSampleSelectionCount, istmBetaMatrix=istmBetaMatrix, istrmTree=istrmTree, istrmEnvr=istrmEnvr)
755
756 if medoidSamples:
757 dictSelectedSamples.setdefault(ConstantsMicropita.c_strRepresentative+"="+ConstantsMicropita.c_custom,[]).extend(medoidSamples)
758 else:
759 logging.info("MicroPITA.funcRunNormalizeSensitiveMethods:: Performing representative selection on normalized data.")
760 for bMetric in lsBetaMetrics:
761
762 #Get representative dissimilarity samples
763 medoidSamples=self.funcGetCentralSamplesByKMedoids(npaMatrix=npaTransposedAbundance, sMetric=bMetric, lsSampleNames=lsSampleNames, iNumberSamplesReturned=iSampleSelectionCount, istmBetaMatrix=istmBetaMatrix, istrmTree=istrmTree, istrmEnvr=istrmEnvr)
764
765 if medoidSamples:
766 dictSelectedSamples.setdefault(self.dictConvertBMetricToMethod.get(bMetric,ConstantsMicropita.c_strRepresentative+"="+bMetric),[]).extend(medoidSamples)
767
768 #Get extreme selection using clusters, tiling
769 if fRunExtreme:
770 logging.info("MicroPITA.funcRunNormalizeSensitiveMethods:: Performing extreme selection on normalized data.")
771 if istmBetaMatrix:
772
773 #Samples for representative dissimilarity
774 #This involves inverting the distance metric,
775 #Taking the dendrogram level of where the number cluster == the number of samples to select
776 #Returning a repersentative sample from each cluster
777 extremeSamples = self.funcSelectExtremeSamplesFromHClust(strBetaMetric=ConstantsMicropita.c_custom, npaAbundanceMatrix=npaTransposedAbundance, lsSampleNames=lsSampleNames, iSelectSampleCount=iSampleSelectionCount, istmBetaMatrix=istmBetaMatrix, istrmTree=istrmTree, istrmEnvr=istrmEnvr)
778
779 #Add selected samples
780 if extremeSamples:
781 dictSelectedSamples.setdefault(ConstantsMicropita.c_strExtreme+"="+ConstantsMicropita.c_custom,[]).extend(extremeSamples)
782
783 else:
784 #Run KMedoids with inverse custom distance metric in normalized space
785 for bMetric in lsInverseBetaMetrics:
786
787 #Samples for representative dissimilarity
788 #This involves inverting the distance metric,
789 #Taking the dendrogram level of where the number cluster == the number of samples to select
790 #Returning a repersentative sample from each cluster
791 extremeSamples = self.funcSelectExtremeSamplesFromHClust(strBetaMetric=bMetric, npaAbundanceMatrix=npaTransposedAbundance, lsSampleNames=lsSampleNames, iSelectSampleCount=iSampleSelectionCount, istmBetaMatrix=istmBetaMatrix, istrmTree=istrmTree, istrmEnvr=istrmEnvr)
792
793 #Add selected samples
794 if extremeSamples:
795 dictSelectedSamples.setdefault(self.dictConvertInvBMetricToMethod.get(bMetric,ConstantsMicropita.c_strExtreme+"="+bMetric),[]).extend(extremeSamples)
796
797 logging.info("MicroPITA.funcRunNormalizeSensitiveMethods:: Selected Samples 2,3b")
798 logging.info(dictSelectedSamples)
799 return dictSelectedSamples
800
801 def funcRun(self, strIDName, strLastMetadataName, istmInput,
802 ostmInputPredictFile, ostmPredictFile, ostmCheckedFile, ostmOutput,
803 cDelimiter, cFeatureNameDelimiter, strFeatureSelection,
804 istmFeatures, iCount, lstrMethods, strLastRowMetadata = None, strLabel = None, strStratify = None,
805 strCustomAlpha = None, strCustomBeta = None, strAlphaMetadata = None, istmBetaMatrix = None, istrmTree = None, istrmEnvr = None,
806 iMinSeqs = ConstantsMicropita.c_liOccurenceFilter[0], iMinSamples = ConstantsMicropita.c_liOccurenceFilter[1], fInvertDiversity = False):
807 """
808 Manages the selection of samples given different metrics.
809
810 :param strIDName: Sample Id metadata row
811 :type: String
812 :param strLastMetadataName: The id of the metadata positioned last in the abundance table.
813 :type: String String metadata id.
814 :param istmInput: File to store input data to supervised methods.
815 :type: FileStream of String file path
816 :param ostmInputPredictFile: File to store distances from supervised methods.
817 :type: FileStream or String file path
818 :param ostmCheckedFile: File to store the AbundanceTable data after it is being checked.
819 :type: FileStream or String file path
820 :param ostmOutPut: File to store sample selection by methods of interest.
821 :type: FileStream or String file path
822 :param cDelimiter: Delimiter of abundance table.
823 :type: Character Char (default TAB).
824 :param cFeatureNameDelimiter: Delimiter of the name of features (for instance if they contain consensus lineages indicating clades).
825 :type: Character (default |).
826 :param stFeatureSelectionMethod: Which method to use to select features in a targeted manner (Using average ranked abundance or average abundance).
827 :type: String (specific values indicated in ConstantsMicropita.lsTargetedFeatureMethodValues).
828 :param istmFeatures: File which holds the features of interest if using targeted feature methodology.
829 :type: FileStream or String file path
830 :param iCount: Number of samples to select in each methods, supervised methods select this amount per label if possible.
831 :type: Integer integer.
832 :param lstrMethods: List of strings indicating selection techniques.
833 :type: List of string method names
834 :param strLabel: The metadata used for supervised labels.
835 :type: String
836 :param strStratify: The metadata used to stratify unsupervised data.
837 :type: String
838 :param strCustomAlpha: Custom alpha diversity metric
839 :type: String
840 :param strCustomBeta: Custom beta diversity metric
841 :type: String
842 :param strAlphaMetadata: Metadata id which is a diveristy metric to use in highest diversity sampling
843 :type: String
844 :param istmBetaMatrix: File containing precalculated beta-diversity matrix for representative sampling
845 :type: FileStream or String file path
846 :param istrmTree: File containing tree for phylogentic beta-diversity analysis
847 :type: FileStream or String file path
848 :param istrmEnvr: File containing environment for phylogentic beta-diversity analysis
849 :type: FileStream or String file path
850 :param iMinSeqs: Minimum sequence in the occurence filter which filters all features not with a minimum number of sequences in each of a minimum number of samples.
851 :type: Integer
852 :param iMinSamples: Minimum sample count for the occurence filter.
853 :type: Integer
854 :param fInvertDiversity: When true will invert diversity measurements before using.
855 :type: boolean
856 :return Selected Samples: Samples selected by methods.
857 Dictionary {"Selection Method":["SampleID","SampleID","SampleID",...]}
858 """
859
860 #Holds the top ranked samples from different metrics
861 #dict[metric name] = [samplename,samplename...]
862 selectedSamples = dict()
863
864 #If a target feature file is given make sure that targeted feature is in the selection methods, if not add
865 if ConstantsMicropita.c_strFeature in lstrMethods:
866 if not istmFeatures:
867 logging.error("MicroPITA.funcRun:: Did not receive both the Targeted feature file and the feature selection method. MicroPITA did not run.")
868 return False
869
870 #Diversity metrics to run
871 #Use custom metrics if specified
872 #Custom beta metrics set to normalized only, custom alpha metrics set to count only
873 diversityMetricsAlpha = [] if strCustomAlpha or strAlphaMetadata else [MicroPITA.c_strInverseSimpsonDiversity]
874 diversityMetricsBeta = [] if istmBetaMatrix else [strCustomBeta] if strCustomBeta else [MicroPITA.c_strBrayCurtisDissimilarity]
875 # inverseDiversityMetricsBeta = [MicroPITA.c_strInvBrayCurtisDissimilarity]
876 diversityMetricsAlphaNoNormalize = [strAlphaMetadata] if strAlphaMetadata else [strCustomAlpha] if strCustomAlpha else []
877 diversityMetricsBetaNoNormalize = []
878 # inverseDiversityMetricsBetaNoNormalize = []
879
880 #Targeted taxa
881 userDefinedTaxa = []
882
883 #Perform different flows flags
884 c_RUN_MAX_DIVERSITY_1 = ConstantsMicropita.c_strDiversity in lstrMethods
885 c_RUN_REPRESENTIVE_DISSIMILARITY_2 = ConstantsMicropita.c_strRepresentative in lstrMethods
886 c_RUN_MAX_DISSIMILARITY_3 = ConstantsMicropita.c_strExtreme in lstrMethods
887 c_RUN_RANK_AVERAGE_USER_4 = False
888 if ConstantsMicropita.c_strFeature in lstrMethods:
889 c_RUN_RANK_AVERAGE_USER_4 = True
890 if not istmFeatures:
891 logging.error("MicroPITA.funcRun:: No taxa file was given for taxa selection.")
892 return False
893 #Read in taxa list, break down to lines and filter out empty strings
894 userDefinedTaxa = filter(None,(s.strip( ) for s in istmFeatures.readlines()))
895 c_RUN_RANDOM_5 = ConstantsMicropita.c_strRandom in lstrMethods
896 c_RUN_DISTINCT = ConstantsMicropita.c_strDistinct in lstrMethods
897 c_RUN_DISCRIMINANT = ConstantsMicropita.c_strDiscriminant in lstrMethods
898
899 #Read in abundance data
900 #Abundance is a structured array. Samples (column) by Taxa (rows) with the taxa id row included as the column index=0
901 #Abundance table object to read in and manage data
902 totalAbundanceTable = AbundanceTable.funcMakeFromFile(xInputFile=istmInput, lOccurenceFilter = [iMinSeqs, iMinSamples],
903 cDelimiter=cDelimiter, sMetadataID=strIDName, sLastMetadataRow=strLastRowMetadata,
904 sLastMetadata=strLastMetadataName, cFeatureNameDelimiter=cFeatureNameDelimiter, xOutputFile=ostmCheckedFile)
905 if not totalAbundanceTable:
906 logging.error("MicroPITA.funcRun:: Could not read in the abundance table. Analysis was not performed."+
907 " This often occurs when the Last Metadata is not specified correctly."+
908 " Please check to make sure the Last Metadata selection is the row of the last metadata,"+
909 " all values after this selection should be microbial measurements and should be numeric.")
910 return False
911
912 lsOriginalLabels = SVM.funcMakeLabels(totalAbundanceTable.funcGetMetadata(strLabel)) if strLabel else strLabel
913
914 dictTotalMetadata = totalAbundanceTable.funcGetMetadataCopy()
915 logging.debug("MicroPITA.funcRun:: Received metadata=" + str(dictTotalMetadata))
916 #If there is only 1 unique value for the labels, do not run the Supervised methods
917 if strLabel and ( len(set(dictTotalMetadata.get(strLabel,[]))) < 2 ):
918 logging.error("The label " + strLabel + " did not have 2 or more values. Labels found=" + str(dictTotalMetadata.get(strLabel,[])))
919 return False
920
921 #Run unsupervised methods###
922 #Stratify the data if need be and drop the old data
923 lStratifiedAbundanceTables = totalAbundanceTable.funcStratifyByMetadata(strStratify) if strStratify else [totalAbundanceTable]
924
925 #For each stratified abundance block or for the unstratfified abundance
926 #Run the unsupervised blocks
927 fAppendSupFiles = False
928 for stratAbundanceTable in lStratifiedAbundanceTables:
929 logging.info("MicroPITA.funcRun:: Running abundance block:"+stratAbundanceTable.funcGetName())
930
931 ###NOT SUMMED, NOT NORMALIZED
932 #Only perform if the data is not yet normalized
933 if not stratAbundanceTable.funcIsNormalized( ):
934 #Need to first work with unnormalized data
935 if c_RUN_MAX_DIVERSITY_1 or c_RUN_REPRESENTIVE_DISSIMILARITY_2 or c_RUN_MAX_DISSIMILARITY_3:
936
937 self._funcRunNormalizeSensitiveMethods(abndData=stratAbundanceTable, iSampleSelectionCount=iCount,
938 dictSelectedSamples=selectedSamples, lsAlphaMetrics=diversityMetricsAlphaNoNormalize,
939 lsBetaMetrics=diversityMetricsBetaNoNormalize,
940 lsInverseBetaMetrics=diversityMetricsBetaNoNormalize,
941 fRunDiversity=c_RUN_MAX_DIVERSITY_1,fRunRepresentative=c_RUN_REPRESENTIVE_DISSIMILARITY_2,
942 fRunExtreme=c_RUN_MAX_DISSIMILARITY_3, strAlphaMetadata=strAlphaMetadata,
943 istrmTree=istrmTree, istrmEnvr=istrmEnvr, fInvertDiversity=fInvertDiversity)
944
945
946 #Generate selection by the rank average of user defined taxa
947 #Expects (Taxa (row) by Samples (column))
948 #Expects a column 0 of taxa id that is skipped
949 #Returns [(sample name,average,rank)]
950 #SUMMED AND NORMALIZED
951 stratAbundanceTable.funcSumClades()
952 #Normalize data at this point
953 stratAbundanceTable.funcNormalize()
954 if c_RUN_RANK_AVERAGE_USER_4:
955 selectedSamples[ConstantsMicropita.c_strFeature] = self.funcSelectTargetedTaxaSamples(abndMatrix=stratAbundanceTable,
956 lsTargetedTaxa=userDefinedTaxa, iSampleSelectionCount=iCount, sMethod=strFeatureSelection)
957 logging.info("MicroPITA.funcRun:: Selected Samples Rank")
958 logging.info(selectedSamples)
959
960 ###SUMMED AND NORMALIZED analysis block
961 #Diversity based metric will move reduce to terminal taxa as needed
962 if c_RUN_MAX_DIVERSITY_1 or c_RUN_REPRESENTIVE_DISSIMILARITY_2 or c_RUN_MAX_DISSIMILARITY_3:
963
964 self._funcRunNormalizeSensitiveMethods(abndData=stratAbundanceTable, iSampleSelectionCount=iCount,
965 dictSelectedSamples=selectedSamples, lsAlphaMetrics=diversityMetricsAlpha,
966 lsBetaMetrics=diversityMetricsBeta,
967 lsInverseBetaMetrics=diversityMetricsBeta,
968 fRunDiversity=c_RUN_MAX_DIVERSITY_1,fRunRepresentative=c_RUN_REPRESENTIVE_DISSIMILARITY_2,
969 fRunExtreme=c_RUN_MAX_DISSIMILARITY_3,
970 istmBetaMatrix=istmBetaMatrix, istrmTree=istrmTree, istrmEnvr=istrmEnvr, fInvertDiversity=fInvertDiversity)
971
972 #5::Select randomly
973 #Expects sampleNames = List of sample names [name, name, name...]
974 if(c_RUN_RANDOM_5):
975 #Select randomly from sample names
976 selectedSamples[ConstantsMicropita.c_strRandom] = self.funcGetRandomSamples(lsSamples=stratAbundanceTable.funcGetSampleNames(), iNumberOfSamplesToReturn=iCount)
977 logging.info("MicroPITA.funcRun:: Selected Samples Random")
978 logging.info(selectedSamples)
979
980 #Perform supervised selection
981 if c_RUN_DISTINCT or c_RUN_DISCRIMINANT:
982 if strLabel:
983 dictSelectionRet = self.funcRunSupervisedDistancesFromCentroids(abundanceTable=stratAbundanceTable,
984 fRunDistinct=c_RUN_DISTINCT, fRunDiscriminant=c_RUN_DISCRIMINANT,
985 xOutputSupFile=ostmInputPredictFile,xPredictSupFile=ostmPredictFile,
986 strSupervisedMetadata=strLabel, iSampleSupSelectionCount=iCount,
987 lsOriginalSampleNames = totalAbundanceTable.funcGetSampleNames(),
988 lsOriginalLabels = lsOriginalLabels,
989 fAppendFiles=fAppendSupFiles)
990
991 [selectedSamples.setdefault(sKey,[]).extend(lValue) for sKey,lValue in dictSelectionRet.items()]
992
993 if not fAppendSupFiles:
994 fAppendSupFiles = True
995 logging.info("MicroPITA.funcRun:: Selected Samples Unsupervised")
996 logging.info(selectedSamples)
997 return selectedSamples
998
999 #Testing: Happy path tested
1000 @staticmethod
1001 def funcWriteSelectionToFile(dictSelection,xOutputFilePath):
1002 """
1003 Writes the selection of samples by method to an output file.
1004
1005 :param dictSelection: The dictionary of selections by method to be written to a file.
1006 :type: Dictionary The dictionary of selections by method {"method":["sample selected","sample selected"...]}
1007 :param xOutputFilePath: FileStream or String path to file inwhich the dictionary is written.
1008 :type: String FileStream or String path to file
1009 """
1010
1011 if not dictSelection:
1012 return
1013
1014 #Open file
1015 f = csv.writer(open(xOutputFilePath,"w") if isinstance(xOutputFilePath, str) else xOutputFilePath, delimiter=ConstantsMicropita.c_outputFileDelim )
1016
1017 #Create output content from dictionary
1018 for sKey in dictSelection:
1019 f.writerow([sKey]+dictSelection[sKey])
1020 logging.debug("MicroPITA.funcRun:: Selected samples output to file:"+str(dictSelection[sKey]))
1021
1022 #Testing: Happy Path tested
1023 @staticmethod
1024 def funcReadSelectionFileToDictionary(xInputFile):
1025 """
1026 Reads in an output selection file from micropita and formats it into a dictionary.
1027
1028 :param xInputFile: String path to file or file stream to read and translate into a dictionary.
1029 {"method":["sample selected","sample selected"...]}
1030 :type: FileStream or String Path to file
1031 :return Dictionary: Samples selected by methods.
1032 Dictionary {"Selection Method":["SampleID","SampleID","SampleID",...]}
1033 """
1034
1035 #Open file
1036 istmReader = csv.reader(open(xInputFile,'r') if isinstance(xInputFile, str) else xInputFile, delimiter = ConstantsMicropita.c_outputFileDelim)
1037
1038 #Dictionary to hold selection data
1039 return dict([(lsLine[0], lsLine[1:]) for lsLine in istmReader])
1040
1041 #Set up arguments reader
1042 argp = argparse.ArgumentParser( prog = "MicroPITA.py",
1043 description = """Selects samples from abundance tables based on various selection schemes.""" )
1044
1045 args = argp.add_argument_group( "Common", "Commonly modified options" )
1046 args.add_argument(ConstantsMicropita.c_strCountArgument,"--num", dest="iCount", metavar = "samples", default = 10, type = int, help = ConstantsMicropita.c_strCountHelp)
1047 args.add_argument("-m","--method", dest = "lstrMethods", metavar = "method", default = [], help = ConstantsMicropita.c_strSelectionTechniquesHelp,
1048 choices = ConstantsMicropita.c_lsAllMethods, action = "append")
1049
1050 args = argp.add_argument_group( "Custom", "Selecting and inputing custom metrics" )
1051 args.add_argument("-a","--alpha", dest = "strAlphaDiversity", metavar = "AlphaDiversity", default = None, help = ConstantsMicropita.c_strCustomAlphaDiversityHelp, choices = Metric.setAlphaDiversities)
1052 args.add_argument("-b","--beta", dest = "strBetaDiversity", metavar = "BetaDiversity", default = None, help = ConstantsMicropita.c_strCustomBetaDiversityHelp, choices = list(Metric.setBetaDiversities)+[Metric.c_strUnifracUnweighted,Metric.c_strUnifracWeighted])
1053 args.add_argument("-q","--alphameta", dest = "strAlphaMetadata", metavar = "AlphaDiversityMetadata", default = None, help = ConstantsMicropita.c_strCustomAlphaDiversityMetadataHelp)
1054 args.add_argument("-x","--betamatrix", dest = "istmBetaMatrix", metavar = "BetaDiversityMatrix", default = None, help = ConstantsMicropita.c_strCustomBetaDiversityMatrixHelp)
1055 args.add_argument("-o","--tree", dest = "istrmTree", metavar = "PhylogeneticTree", default = None, help = ConstantsMicropita.c_strCustomPhylogeneticTreeHelp)
1056 args.add_argument("-i","--envr", dest = "istrmEnvr", metavar = "EnvironmentFile", default = None, help = ConstantsMicropita.c_strCustomEnvironmentFileHelp)
1057 args.add_argument("-f","--invertDiversity", dest = "fInvertDiversity", action="store_true", default = False, help = ConstantsMicropita.c_strInvertDiversityHelp)
1058
1059 args = argp.add_argument_group( "Miscellaneous", "Row/column identifiers and feature targeting options" )
1060 args.add_argument("-d",ConstantsMicropita.c_strIDNameArgument, dest="strIDName", metavar="sample_id", help= ConstantsMicropita.c_strIDNameHelp)
1061 args.add_argument("-l",ConstantsMicropita.c_strLastMetadataNameArgument, dest="strLastMetadataName", metavar = "metadata_id", default = None,
1062 help= ConstantsMicropita.c_strLastMetadataNameHelp)
1063 args.add_argument("-r",ConstantsMicropita.c_strTargetedFeatureMethodArgument, dest="strFeatureSelection", metavar="targeting_method", default=ConstantsMicropita.lsTargetedFeatureMethodValues[0],
1064 choices=ConstantsMicropita.lsTargetedFeatureMethodValues, help= ConstantsMicropita.c_strTargetedFeatureMethodHelp)
1065 args.add_argument("-t",ConstantsMicropita.c_strTargetedSelectionFileArgument, dest="istmFeatures", metavar="feature_file", type=argparse.FileType("rU"), help=ConstantsMicropita.c_strTargetedSelectionFileHelp)
1066 args.add_argument("-w",ConstantsMicropita.c_strFeatureMetadataArgument, dest="strLastFeatureMetadata", metavar="Last_Feature_Metadata", default=None, help=ConstantsMicropita.c_strFeatureMetadataHelp)
1067
1068 args = argp.add_argument_group( "Data labeling", "Metadata IDs for strata and supervised label values" )
1069 args.add_argument("-e",ConstantsMicropita.c_strSupervisedLabelArgument, dest="strLabel", metavar= "supervised_id", help=ConstantsMicropita.c_strSupervisedLabelHelp)
1070 args.add_argument("-s",ConstantsMicropita.c_strUnsupervisedStratifyMetadataArgument, dest="strUnsupervisedStratify", metavar="stratify_id",
1071 help= ConstantsMicropita.c_strUnsupervisedStratifyMetadataHelp)
1072
1073 args = argp.add_argument_group( "File formatting", "Rarely modified file formatting options" )
1074 args.add_argument("-j",ConstantsMicropita.c_strFileDelimiterArgument, dest="cFileDelimiter", metavar="column_delimiter", default="\t", help=ConstantsMicropita.c_strFileDelimiterHelp)
1075 args.add_argument("-k",ConstantsMicropita.c_strFeatureNameDelimiterArgument, dest="cFeatureNameDelimiter", metavar="taxonomy_delimiter", default="|", help=ConstantsMicropita.c_strFeatureNameDelimiterHelp)
1076
1077 args = argp.add_argument_group( "Debugging", "Debugging options - modify at your own risk!" )
1078 args.add_argument("-v",ConstantsMicropita.c_strLoggingArgument, dest="strLogLevel", metavar = "log_level", default="WARNING",
1079 choices=ConstantsMicropita.c_lsLoggingChoices, help= ConstantsMicropita.c_strLoggingHelp)
1080 args.add_argument("-c",ConstantsMicropita.c_strCheckedAbundanceFileArgument, dest="ostmCheckedFile", metavar = "output_qc", type = argparse.FileType("w"), help = ConstantsMicropita.c_strCheckedAbundanceFileHelp)
1081 args.add_argument("-g",ConstantsMicropita.c_strLoggingFileArgument, dest="ostmLoggingFile", metavar = "output_log", type = argparse.FileType("w"), help = ConstantsMicropita.c_strLoggingFileHelp)
1082 args.add_argument("-u",ConstantsMicropita.c_strSupervisedInputFile, dest="ostmInputPredictFile", metavar = "output_scaled", type = argparse.FileType("w"), help = ConstantsMicropita.c_strSupervisedInputFileHelp)
1083 args.add_argument("-p",ConstantsMicropita.c_strSupervisedPredictedFile, dest="ostmPredictFile", metavar = "output_labels", type = argparse.FileType("w"), help = ConstantsMicropita.c_strSupervisedPredictedFileHelp)
1084
1085 argp.add_argument("istmInput", metavar = "input.pcl/biome", type = argparse.FileType("rU"), help = ConstantsMicropita.c_strAbundanceFileHelp,
1086 default = sys.stdin)
1087 argp.add_argument("ostmOutput", metavar = "output.txt", type = argparse.FileType("w"), help = ConstantsMicropita.c_strGenericOutputDataFileHelp,
1088 default = sys.stdout)
1089
1090 __doc__ = "::\n\n\t" + argp.format_help( ).replace( "\n", "\n\t" ) + __doc__
1091
1092 def _main( ):
1093 args = argp.parse_args( )
1094
1095 #Set up logger
1096 iLogLevel = getattr(logging, args.strLogLevel.upper(), None)
1097 logging.basicConfig(stream = args.ostmLoggingFile if args.ostmLoggingFile else sys.stderr, filemode = 'w', level=iLogLevel)
1098
1099 #Run micropita
1100 logging.info("MicroPITA:: Start microPITA")
1101 microPITA = MicroPITA()
1102
1103 #Argparse will append to the default but will not remove the default so I do this here
1104 if not len(args.lstrMethods):
1105 args.lstrMethods = [ConstantsMicropita.c_strRepresentative]
1106
1107 dictSelectedSamples = microPITA.funcRun(
1108 strIDName = args.strIDName,
1109 strLastMetadataName = args.strLastMetadataName,
1110 istmInput = args.istmInput,
1111 ostmInputPredictFile = args.ostmInputPredictFile,
1112 ostmPredictFile = args.ostmPredictFile,
1113 ostmCheckedFile = args.ostmCheckedFile,
1114 ostmOutput = args.ostmOutput,
1115 cDelimiter = args.cFileDelimiter,
1116 cFeatureNameDelimiter = args.cFeatureNameDelimiter,
1117 istmFeatures = args.istmFeatures,
1118 strFeatureSelection = args.strFeatureSelection,
1119 iCount = args.iCount,
1120 strLastRowMetadata = args.strLastFeatureMetadata,
1121 strLabel = args.strLabel,
1122 strStratify = args.strUnsupervisedStratify,
1123 strCustomAlpha = args.strAlphaDiversity,
1124 strCustomBeta = args.strBetaDiversity,
1125 strAlphaMetadata = args.strAlphaMetadata,
1126 istmBetaMatrix = args.istmBetaMatrix,
1127 istrmTree = args.istrmTree,
1128 istrmEnvr = args.istrmEnvr,
1129 lstrMethods = args.lstrMethods,
1130 fInvertDiversity = args.fInvertDiversity
1131 )
1132
1133 if not dictSelectedSamples:
1134 logging.error("MicroPITA:: Error, did not get a result from analysis.")
1135 return -1
1136 logging.info("End microPITA")
1137
1138 #Log output for debugging
1139 logging.debug("MicroPITA:: Returned the following samples:"+str(dictSelectedSamples))
1140
1141 #Write selection to file
1142 microPITA.funcWriteSelectionToFile(dictSelection=dictSelectedSamples, xOutputFilePath=args.ostmOutput)
1143
1144 if __name__ == "__main__":
1145 _main( )