← Back to branch summary

~ubuntu-branches/ubuntu/trusty/qiime/trusty

« back to all changes in this revision

Viewing changes to .pc/fix_t_test_bug/qiime/compare_alpha_diversity.py

  • Committer: Package Import Robot
  • Author(s): Andreas Tille
  • Date: 2013-06-17 18:28:26 UTC
  • mfrom: (9.1.2 sid)
  • Revision ID: package-import@ubuntu.com-20130617182826-376az5ad080a0sfe
Tags: 1.7.0+dfsg-1
Upload preparations done for BioLinux to Debian

Show diffs side-by-side

added added

removed removed

Lines of Context:
.pc/fix_t_test_bug/qiime/compare_alpha_diversity.py
1
 
#!/usr/bin/env python
2
 
# File created on 19 May 2011
3
 
from __future__ import division
4
 
 
5
 
__author__ = "William Van Treuren"
6
 
__copyright__ = "Copyright 2011, The QIIME project"
7
 
__credits__ = ["William Van Treuren","Greg Caporaso"]
8
 
__license__ = "GPL"
9
 
__version__ = "1.5.0"
10
 
__maintainer__ = "William Van Treuren"
11
 
__email__ = "vantreur@colorado.edu"
12
 
__status__ = "Release"
13
 
 
14
 
from numpy import array, isnan
15
 
from qiime.parse import parse_mapping_file_to_dict, parse_rarefaction
16
 
from cogent.maths.stats.test import t_two_sample
17
 
 
18
 
def make_value_pairs_from_category(mapping_data, category):
19
 
    """creates all pairs of unique category values from mapping data
20
 
    
21
 
    input:
22
 
        mapping_data - a nested dictionary that maps SampleIds to
23
 
        descriptor categories (e.g. {'Id1': {'Weight':'Fat'}}
24
 
    
25
 
        category - a string which specifies a category found in the 
26
 
        mapping data (e.g. 'Obese')
27
 
    
28
 
    output:
29
 
        unique_pairs - a list of unique pairs of the values specified by
30
 
        category
31
 
        (e.g. [('Obese','Fat'),('Obese','notFat'),('Fat','notFat')]
32
 
    
33
 
    """
34
 
    
35
 
    #gets the keys of the mapping file dictionary. corresponds to the 
36
 
    #names of the individuals in the mapping data file
37
 
    
38
 
    keys = mapping_data.keys()
39
 
       
40
 
    categories = []
41
 
    for key in keys:
42
 
        try:
43
 
            categories.append(mapping_data[key][category])
44
 
        except KeyError:
45
 
            raise ValueError(('the specified category ({0}) was '+\
46
 
                    'not found in the mapping file.').format(category))
47
 
    #strip duplicate values from this list
48
 
    
49
 
    unique_vals = []
50
 
    for val in categories:
51
 
        if unique_vals.count(val)==0:
52
 
            unique_vals.append(val)
53
 
            
54
 
    #create and populate list of unique pairs of category values
55
 
    
56
 
    unique_pairs = []
57
 
    for i in range(len(unique_vals)):
58
 
        for j in unique_vals[i+1:]:
59
 
            unique_pairs.append((unique_vals[i], j))
60
 
    
61
 
    return unique_pairs
62
 
            
63
 
def make_category_values_Id_dict(mapping_data, category):
64
 
    """makes dict lists all SampleIds that have given category value
65
 
    
66
 
    input:
67
 
        mapping_data - a nested dictionary that maps SampleID to
68
 
        descriptor categories (e.g. {'Id1': {'Weight':'Fat}}
69
 
    
70
 
        category - a string which specifies a category found in the 
71
 
        mapping data (e.g. 'Weight')
72
 
    
73
 
    output:
74
 
        cat_val_Ids - a dictionary with category values as keys
75
 
        and a list of SampleIds which have the specific category value
76
 
        as values 
77
 
        (e.g. {'Fat':['Id1','Id2'],'notFat':['Id3'],'Obese':['Id4']}).
78
 
    
79
 
    """
80
 
    
81
 
    keys = mapping_data.keys()
82
 
    
83
 
    # create and populate list of all the different values for the 
84
 
    # category that was specified 
85
 
    
86
 
    categories = []
87
 
    for key in keys:
88
 
        categories.append(mapping_data[key][category])
89
 
    
90
 
    #strip duplicate values from this list
91
 
    
92
 
    unique_vals = []
93
 
    for val in categories:
94
 
        if unique_vals.count(val)==0:
95
 
            unique_vals.append(val)
96
 
    
97
 
    #make a dictionary with keys that are the possible values of the 
98
 
    #category that was specified.
99
 
    
100
 
    cat_val_Ids = {}
101
 
    for val in unique_vals:
102
 
        cat_val_Ids[val] = []
103
 
    
104
 
    #populate the cat_val dict with Id's which have proper category 
105
 
    
106
 
    for key in keys:
107
 
        for val in cat_val_Ids.keys():
108
 
            if mapping_data[key][category] == val:
109
 
                cat_val_Ids[val].append(key)
110
 
    
111
 
    return cat_val_Ids
112
 
        
113
 
def map_category_value_pairs_to_Ids(value_pairs, cat_val_Ids):
114
 
    """maps category value pairs to Id's which have that category value
115
 
    
116
 
    input:
117
 
        value_pairs - a list of pairs of categories (e.g. 
118
 
        [('Obese','Fat'),('Fat','notFat'),('Obese','notFat')]
119
 
        
120
 
        cat_val_Ids - a dictionary with category values as keys
121
 
        and a list of SampleId's which have the specific category value 
122
 
        as values
123
 
        (e.g. {'Fat':['Id1','Id2'],'notFat':['Id3'],'Obese':['Id4']}).
124
 
    
125
 
    output:
126
 
        mapped_pairs - the list of value_pairs with the values replaced
127
 
        by the SampleIds which have the values specified in the pair e.g
128
 
        [(['Id4'],['Id1','Id2']),(['Id1','Id2'],\
129
 
        ['Id3]),(['Id4'],['Id3])] 
130
 
    """
131
 
    
132
 
    mapped_pairs = []
133
 
    
134
 
    for pair in value_pairs:
135
 
        mapped_pairs.append((cat_val_Ids[pair[0]],cat_val_Ids[pair[1]]))
136
 
    
137
 
    return mapped_pairs
138
 
 
139
 
 
140
 
def make_SampleIds_rarefaction_columns_dict(rarefaction_list):
141
 
    """maps SampleId to column in parsed rarefaction file output
142
 
    
143
 
    input:
144
 
        rarefaction_list - ouput of parse_rarefaction.py. a nested list
145
 
        of scores and SampleIds and other fields. 
146
 
    
147
 
    output:
148
 
        map_from_Id_to_col - a dict which has as keys the SampleIds, and
149
 
        as values the col they are in the in the parsed rarefaction list
150
 
    """
151
 
    
152
 
    map_from_Id_to_col = {}
153
 
    
154
 
    # the first 3 entries in the rarefaction_list are not SampleIDs
155
 
    # so we dump them
156
 
    
157
 
    Ids = rarefaction_list[0][3:]
158
 
    
159
 
    for Id in Ids:
160
 
        map_from_Id_to_col[Id] = Ids.index(Id)
161
 
    
162
 
    return map_from_Id_to_col
163
 
        
164
 
    
165
 
def extract_rarefaction_scores_at_depth(depth, rarefaction_list):
166
 
    """makes rarefaction matrix with row=iteration and col=SampleId
167
 
    
168
 
    input:
169
 
        depth - an integer which corresponds to the depth of the
170
 
        rarefaction. also called the "sequences per sample" in the 
171
 
        rarefaction file
172
 
        
173
 
        rarefaction_list - ouput of parse_rarefaction.py. a nested list
174
 
        of scores and SampleIds and other fields. 
175
 
    
176
 
    output:
177
 
        result - a matrix with rows=rarefaction scores at a given depth
178
 
        and iteration, and cols=SampleIds.
179
 
    
180
 
    """
181
 
    # make and populate an array that has as rows rarefaction values
182
 
    # at the same depth and iteration and as cols SampleIds.
183
 
    
184
 
    score_matrix = []
185
 
    
186
 
    # the 4th element of rarefaction_list is a list of scores for each
187
 
    # different SampleId
188
 
    
189
 
    for line in rarefaction_list[3]:
190
 
        if line[0] == depth:
191
 
            # the first two elements are just rarefaction depth and 
192
 
            # iteration, throw these away
193
 
            score_matrix.append(line[2:])
194
 
    
195
 
    # raise error if rarefaction depth not found in rarefaction file
196
 
    if score_matrix == []:
197
 
        raise ValueError(('Specified depth ({0}) was not found in '+\
198
 
                    'the rarefaction file.').format(depth))
199
 
                        
200
 
    
201
 
    score_matrix_elements = []
202
 
    
203
 
    for line in score_matrix:
204
 
        score_matrix_elements.append(line)
205
 
    
206
 
    result = array(score_matrix_elements) 
207
 
    
208
 
    # raise error if any rarefaction score at spec. depth is Nan
209
 
    if isnan(result).any():
210
 
        raise ValueError(('Specified depth ({0}) has NaNs for some '+\
211
 
                            'rarefaction scores.').format(depth))
212
 
    
213
 
    return result
214
 
 
215
 
 
216
 
def convert_SampleIds_to_rarefaction_mtx(chosen_SampleIds,score_matrix,\
217
 
                                        map_from_SampleIds_to_cols):
218
 
    """converts list of SampleIDs to score mtx from rarefaction file
219
 
    
220
 
    input:
221
 
        chosen_SampleIds - a list of SampleIds
222
 
        
223
 
        score_matrix - a matrix created by
224
 
        extract_rarefaction_scores_at_depth which represents the
225
 
        rarefaction scores for a given depth
226
 
        
227
 
        map_from_SampleIds_to_cols - a dict which maps a SampleId to 
228
 
        the column its scores are in in the score matrix
229
 
    
230
 
    output:
231
 
        reduced_scores_matrix - a matrix which is the input scores mtx
232
 
        with only the cols that correspond to the chosen_SampleIds
233
 
    
234
 
    """
235
 
    #create and populate a list that specifies the r_array columns which
236
 
    #correspond to the name_list
237
 
    
238
 
    cols=[]
239
 
    
240
 
    for Id in chosen_SampleIds:
241
 
        cols.append(map_from_SampleIds_to_cols[Id])
242
 
        
243
 
    # grab only the columns we need based on a passed list of names and
244
 
    # a dictionary to convert between those names and the proper cols
245
 
    reduced_scores_matrix = score_matrix.take(cols, axis=1)
246
 
    
247
 
    return reduced_scores_matrix
248
 
    
249
 
    
250
 
 
251
 
def compare_alpha_diversities(rarefaction_lines, mapping_lines, 
252
 
                              category, depth):
253
 
    """compares alpha diversities
254
 
    
255
 
    inputs:
256
 
        rarefaction_file - rarefaction file which gives scores for 
257
 
        various rarefactions and depths
258
 
        
259
 
        mapping_file - file that has ID's and categories that the ID's
260
 
        fall in
261
 
        
262
 
        category - the category to be compared, is a string
263
 
        
264
 
        depth - the depth of the rarefaction_file to use, is an integer
265
 
    
266
 
    outputs:
267
 
        results - a nested dictionary which specifies the category as
268
 
        the top level key, and as its value, dictionaries which give the
269
 
        results of the t_two_sample test for all unique pairs of values
270
 
        in the specified category
271
 
    
272
 
    """
273
 
     
274
 
    rarefaction_data = parse_rarefaction(rarefaction_lines)
275
 
    mapping_data = parse_mapping_file_to_dict(mapping_lines)[0]
276
 
    value_pairs = make_value_pairs_from_category(mapping_data, category)
277
 
    
278
 
    category_values_Ids = make_category_values_Id_dict(mapping_data, 
279
 
                                                       category)
280
 
    
281
 
    SampleId_pairs = map_category_value_pairs_to_Ids(value_pairs,
282
 
                                                    category_values_Ids)
283
 
    
284
 
    map_from_Id_to_col = make_SampleIds_rarefaction_columns_dict(
285
 
                                                       rarefaction_data)
286
 
    
287
 
    reduced_rarefaction_mtx = extract_rarefaction_scores_at_depth(depth,
288
 
                                                       rarefaction_data)
289
 
    
290
 
    results = {category:{}}
291
 
    
292
 
    for pair in range(len(SampleId_pairs)):
293
 
        i=(convert_SampleIds_to_rarefaction_mtx(SampleId_pairs[pair][0],
294
 
                           reduced_rarefaction_mtx, map_from_Id_to_col))
295
 
        
296
 
        j=(convert_SampleIds_to_rarefaction_mtx(SampleId_pairs[pair][1],
297
 
                           reduced_rarefaction_mtx, map_from_Id_to_col))
298
 
        
299
 
        results[category][(str(value_pairs[pair][0]),
300
 
                           str(value_pairs[pair][1]))] =\
301
 
                          t_two_sample(i,j)
302
 
    
303
 
    return results
304