annotate change_o/DefineClones.py @ 72:c787e79e8b90 draft

Uploaded
author davidvanzessen
date Tue, 05 Feb 2019 03:26:41 -0500
parents 22dddabe3637
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1 #!/usr/bin/env python3
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
2 """
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
3 Assign Ig sequences into clones
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
4 """
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
5 # Info
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
6 __author__ = 'Namita Gupta, Jason Anthony Vander Heiden, Gur Yaari, Mohamed Uduman'
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
7 from changeo import __version__, __date__
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
8
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
9 # Imports
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
10 import os
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
11 import re
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
12 import sys
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
13 import csv
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
14 import numpy as np
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
15 from argparse import ArgumentParser
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
16 from collections import OrderedDict
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
17 from itertools import chain
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
18 from textwrap import dedent
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
19 from time import time
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
20 from Bio import pairwise2
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
21 from Bio.Seq import translate
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
22
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
23 # Presto and changeo imports
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
24 from presto.Defaults import default_out_args
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
25 from presto.IO import getFileType, getOutputHandle, printLog, printProgress
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
26 from presto.Multiprocessing import manageProcesses
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
27 from presto.Sequence import getDNAScoreDict
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
28 from changeo.Commandline import CommonHelpFormatter, getCommonArgParser, parseCommonArgs
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
29 from changeo.Distance import distance_models, calcDistances, formClusters
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
30 from changeo.IO import getDbWriter, readDbFile, countDbFile
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
31 from changeo.Multiprocessing import DbData, DbResult
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
32
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
33 ## Set maximum field size for csv.reader
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
34 csv.field_size_limit(sys.maxsize)
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
35
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
36 # Defaults
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
37 default_translate = False
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
38 default_distance = 0.0
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
39 default_index_mode = 'gene'
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
40 default_index_action = 'set'
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
41 default_bygroup_model = 'ham'
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
42 default_hclust_model = 'chen2010'
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
43 default_seq_field = 'JUNCTION'
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
44 default_norm = 'len'
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
45 default_sym = 'avg'
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
46 default_linkage = 'single'
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
47 choices_bygroup_model = ('ham', 'aa', 'hh_s1f', 'hh_s5f', 'mk_rs1nf', 'mk_rs5nf', 'hs1f_compat', 'm1n_compat')
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
48
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
49
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
50 def indexByIdentity(index, key, rec, fields=None):
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
51 """
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
52 Updates a preclone index with a simple key
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
53
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
54 Arguments:
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
55 index = preclone index from indexJunctions
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
56 key = index key
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
57 rec = IgRecord to add to the index
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
58 fields = additional annotation fields to use to group preclones;
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
59 if None use only V, J and junction length
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
60
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
61 Returns:
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
62 None. Updates index with new key and records.
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
63 """
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
64 index.setdefault(tuple(key), []).append(rec)
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
65
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
66
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
67 def indexByUnion(index, key, rec, fields=None):
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
68 """
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
69 Updates a preclone index with the union of nested keys
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
70
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
71 Arguments:
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
72 index = preclone index from indexJunctions
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
73 key = index key
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
74 rec = IgRecord to add to the index
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
75 fields = additional annotation fields to use to group preclones;
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
76 if None use only V, J and junction length
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
77
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
78 Returns:
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
79 None. Updates index with new key and records.
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
80 """
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
81 # List of values for this/new key
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
82 val = [rec]
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
83 f_range = list(range(2, 3 + (len(fields) if fields else 0)))
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
84
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
85 # See if field/junction length combination exists in index
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
86 outer_dict = index
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
87 for field in f_range:
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
88 try:
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
89 outer_dict = outer_dict[key[field]]
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
90 except (KeyError):
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
91 outer_dict = None
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
92 break
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
93 # If field combination exists, look through Js
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
94 j_matches = []
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
95 if outer_dict is not None:
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
96 for j in outer_dict.keys():
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
97 if not set(key[1]).isdisjoint(set(j)):
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
98 key[1] = tuple(set(key[1]).union(set(j)))
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
99 j_matches += [j]
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
100 # If J overlap exists, look through Vs for each J
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
101 for j in j_matches:
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
102 v_matches = []
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
103 # Collect V matches for this J
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
104 for v in outer_dict[j].keys():
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
105 if not set(key[0]).isdisjoint(set(v)):
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
106 key[0] = tuple(set(key[0]).union(set(v)))
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
107 v_matches += [v]
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
108 # If there are V overlaps for this J, pop them out
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
109 if v_matches:
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
110 val += list(chain(*(outer_dict[j].pop(v) for v in v_matches)))
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
111 # If the J dict is now empty, remove it
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
112 if not outer_dict[j]:
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
113 outer_dict.pop(j, None)
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
114
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
115 # Add value(s) into index nested dictionary
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
116 # OMG Python pointers are the best!
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
117 # Add field dictionaries into index
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
118 outer_dict = index
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
119 for field in f_range:
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
120 outer_dict.setdefault(key[field], {})
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
121 outer_dict = outer_dict[key[field]]
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
122 # Add J, then V into index
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
123 if key[1] in outer_dict:
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
124 outer_dict[key[1]].update({key[0]: val})
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
125 else:
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
126 outer_dict[key[1]] = {key[0]: val}
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
127
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
128
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
129 def indexJunctions(db_iter, fields=None, mode=default_index_mode,
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
130 action=default_index_action):
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
131 """
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
132 Identifies preclonal groups by V, J and junction length
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
133
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
134 Arguments:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
135 db_iter = an iterator of IgRecords defined by readDbFile
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
136 fields = additional annotation fields to use to group preclones;
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
137 if None use only V, J and junction length
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
138 mode = specificity of alignment call to use for assigning preclones;
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
139 one of ('allele', 'gene')
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
140 action = how to handle multiple value fields when assigning preclones;
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
141 one of ('first', 'set')
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
142
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
143 Returns:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
144 a dictionary of {(V, J, junction length):[IgRecords]}
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
145 """
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
146 # print(fields)
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
147 # Define functions for grouping keys
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
148 if mode == 'allele' and fields is None:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
149 def _get_key(rec, act):
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
150 return [rec.getVAllele(act), rec.getJAllele(act),
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
151 None if rec.junction is None else len(rec.junction)]
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
152 elif mode == 'gene' and fields is None:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
153 def _get_key(rec, act):
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
154 return [rec.getVGene(act), rec.getJGene(act),
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
155 None if rec.junction is None else len(rec.junction)]
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
156 elif mode == 'allele' and fields is not None:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
157 def _get_key(rec, act):
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
158 vdj = [rec.getVAllele(act), rec.getJAllele(act),
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
159 None if rec.junction is None else len(rec.junction)]
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
160 ann = [rec.toDict().get(k, None) for k in fields]
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
161 return list(chain(vdj, ann))
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
162 elif mode == 'gene' and fields is not None:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
163 def _get_key(rec, act):
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
164 vdj = [rec.getVGene(act), rec.getJGene(act),
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
165 None if rec.junction is None else len(rec.junction)]
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
166 ann = [rec.toDict().get(k, None) for k in fields]
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
167 return list(chain(vdj, ann))
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
168
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
169 # Function to flatten nested dictionary
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
170 def _flatten_dict(d, parent_key=''):
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
171 items = []
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
172 for k, v in d.items():
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
173 new_key = parent_key + [k] if parent_key else [k]
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
174 if isinstance(v, dict):
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
175 items.extend(_flatten_dict(v, new_key).items())
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
176 else:
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
177 items.append((new_key, v))
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
178 flat_dict = {None if None in i[0] else tuple(i[0]): i[1] for i in items}
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
179 return flat_dict
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
180
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
181 if action == 'first':
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
182 index_func = indexByIdentity
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
183 elif action == 'set':
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
184 index_func = indexByUnion
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
185 else:
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
186 sys.stderr.write('Unrecognized action: %s.\n' % action)
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
187
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
188 start_time = time()
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
189 clone_index = {}
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
190 rec_count = 0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
191 for rec in db_iter:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
192 key = _get_key(rec, action)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
193
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
194 # Print progress
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
195 if rec_count == 0:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
196 print('PROGRESS> Grouping sequences')
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
197
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
198 printProgress(rec_count, step=1000, start_time=start_time)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
199 rec_count += 1
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
200
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
201 # Assigned passed preclone records to key and failed to index None
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
202 if all([k is not None and k != '' for k in key]):
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
203 # Update index dictionary
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
204 index_func(clone_index, key, rec, fields)
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
205 else:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
206 clone_index.setdefault(None, []).append(rec)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
207
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
208 printProgress(rec_count, step=1000, start_time=start_time, end=True)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
209
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
210 if action == 'set':
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
211 clone_index = _flatten_dict(clone_index)
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
212
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
213 return clone_index
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
214
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
215
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
216 def distanceClones(records, model=default_bygroup_model, distance=default_distance,
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
217 dist_mat=None, norm=default_norm, sym=default_sym,
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
218 linkage=default_linkage, seq_field=default_seq_field):
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
219 """
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
220 Separates a set of IgRecords into clones
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
221
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
222 Arguments:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
223 records = an iterator of IgRecords
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
224 model = substitution model used to calculate distance
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
225 distance = the distance threshold to assign clonal groups
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
226 dist_mat = pandas DataFrame of pairwise nucleotide or amino acid distances
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
227 norm = normalization method
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
228 sym = symmetry method
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
229 linkage = type of linkage
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
230 seq_field = sequence field used to calculate distance between records
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
231
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
232 Returns:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
233 a dictionary of lists defining {clone number: [IgRecords clonal group]}
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
234 """
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
235 # Get distance matrix if not provided
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
236 if dist_mat is None:
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
237 try:
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
238 dist_mat = distance_models[model]
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
239 except KeyError:
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
240 sys.exit('Unrecognized distance model: %s' % args_dict['model'])
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
241
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
242 # TODO: can be cleaned up with abstract model class
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
243 # Determine length of n-mers
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
244 if model in ['hs1f_compat', 'm1n_compat', 'aa', 'ham', 'hh_s1f', 'mk_rs1nf']:
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
245 nmer_len = 1
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
246 elif model in ['hh_s5f', 'mk_rs5nf']:
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
247 nmer_len = 5
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
248 else:
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
249 sys.exit('Unrecognized distance model: %s.\n' % model)
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
250
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
251 # Define unique junction mapping
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
252 seq_map = {}
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
253 for ig in records:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
254 seq = ig.getSeqField(seq_field)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
255 # Check if sequence length is 0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
256 if len(seq) == 0:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
257 return None
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
258
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
259 seq = re.sub('[\.-]', 'N', str(seq))
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
260 if model == 'aa': seq = translate(seq)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
261
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
262 seq_map.setdefault(seq, []).append(ig)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
263
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
264 # Process records
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
265 if len(seq_map) == 1:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
266 return {1:records}
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
267
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
268 # Define sequences
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
269 seqs = list(seq_map.keys())
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
270
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
271 # Calculate pairwise distance matrix
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
272 dists = calcDistances(seqs, nmer_len, dist_mat, sym=sym, norm=norm)
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
273
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
274 # Perform hierarchical clustering
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
275 clusters = formClusters(dists, linkage, distance)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
276
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
277 # Turn clusters into clone dictionary
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
278 clone_dict = {}
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
279 for i, c in enumerate(clusters):
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
280 clone_dict.setdefault(c, []).extend(seq_map[seqs[i]])
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
281
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
282 return clone_dict
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
283
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
284
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
285 def distChen2010(records):
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
286 """
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
287 Calculate pairwise distances as defined in Chen 2010
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
288
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
289 Arguments:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
290 records = list of IgRecords where first is query to be compared to others in list
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
291
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
292 Returns:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
293 list of distances
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
294 """
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
295 # Pull out query sequence and V/J information
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
296 query = records.popitem(last=False)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
297 query_cdr3 = query.junction[3:-3]
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
298 query_v_allele = query.getVAllele()
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
299 query_v_gene = query.getVGene()
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
300 query_v_family = query.getVFamily()
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
301 query_j_allele = query.getJAllele()
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
302 query_j_gene = query.getJGene()
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
303 # Create alignment scoring dictionary
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
304 score_dict = getDNAScoreDict()
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
305
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
306 scores = [0]*len(records)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
307 for i in range(len(records)):
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
308 ld = pairwise2.align.globalds(query_cdr3, records[i].junction[3:-3],
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
309 score_dict, -1, -1, one_alignment_only=True)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
310 # Check V similarity
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
311 if records[i].getVAllele() == query_v_allele: ld += 0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
312 elif records[i].getVGene() == query_v_gene: ld += 1
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
313 elif records[i].getVFamily() == query_v_family: ld += 3
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
314 else: ld += 5
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
315 # Check J similarity
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
316 if records[i].getJAllele() == query_j_allele: ld += 0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
317 elif records[i].getJGene() == query_j_gene: ld += 1
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
318 else: ld += 3
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
319 # Divide by length
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
320 scores[i] = ld/max(len(records[i].junction[3:-3]), query_cdr3)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
321
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
322 return scores
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
323
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
324
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
325 def distAdemokun2011(records):
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
326 """
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
327 Calculate pairwise distances as defined in Ademokun 2011
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
328
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
329 Arguments:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
330 records = list of IgRecords where first is query to be compared to others in list
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
331
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
332 Returns:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
333 list of distances
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
334 """
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
335 # Pull out query sequence and V family information
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
336 query = records.popitem(last=False)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
337 query_cdr3 = query.junction[3:-3]
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
338 query_v_family = query.getVFamily()
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
339 # Create alignment scoring dictionary
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
340 score_dict = getDNAScoreDict()
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
341
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
342 scores = [0]*len(records)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
343 for i in range(len(records)):
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
344
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
345 if abs(len(query_cdr3) - len(records[i].junction[3:-3])) > 10:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
346 scores[i] = 1
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
347 elif query_v_family != records[i].getVFamily():
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
348 scores[i] = 1
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
349 else:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
350 ld = pairwise2.align.globalds(query_cdr3, records[i].junction[3:-3],
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
351 score_dict, -1, -1, one_alignment_only=True)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
352 scores[i] = ld/min(len(records[i].junction[3:-3]), query_cdr3)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
353
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
354 return scores
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
355
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
356
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
357 def hierClust(dist_mat, method='chen2010'):
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
358 """
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
359 Calculate hierarchical clustering
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
360
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
361 Arguments:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
362 dist_mat = square-formed distance matrix of pairwise CDR3 comparisons
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
363
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
364 Returns:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
365 list of cluster ids
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
366 """
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
367 if method == 'chen2010':
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
368 clusters = formClusters(dist_mat, 'average', 0.32)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
369 elif method == 'ademokun2011':
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
370 clusters = formClusters(dist_mat, 'complete', 0.25)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
371 else: clusters = np.ones(dist_mat.shape[0])
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
372
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
373 return clusters
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
374
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
375 # TODO: Merge duplicate feed, process and collect functions.
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
376 def feedQueue(alive, data_queue, db_file, group_func, group_args={}):
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
377 """
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
378 Feeds the data queue with Ig records
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
379
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
380 Arguments:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
381 alive = a multiprocessing.Value boolean controlling whether processing continues
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
382 if False exit process
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
383 data_queue = a multiprocessing.Queue to hold data for processing
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
384 db_file = the Ig record database file
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
385 group_func = the function to use for assigning preclones
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
386 group_args = a dictionary of arguments to pass to group_func
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
387
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
388 Returns:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
389 None
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
390 """
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
391 # Open input file and perform grouping
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
392 try:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
393 # Iterate over Ig records and assign groups
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
394 db_iter = readDbFile(db_file)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
395 clone_dict = group_func(db_iter, **group_args)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
396 except:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
397 #sys.stderr.write('Exception in feeder grouping step\n')
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
398 alive.value = False
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
399 raise
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
400
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
401 # Add groups to data queue
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
402 try:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
403 #print 'START FEED', alive.value
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
404 # Iterate over groups and feed data queue
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
405 clone_iter = iter(clone_dict.items())
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
406 while alive.value:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
407 # Get data from queue
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
408 if data_queue.full(): continue
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
409 else: data = next(clone_iter, None)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
410 # Exit upon reaching end of iterator
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
411 if data is None: break
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
412 #print "FEED", alive.value, k
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
413
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
414 # Feed queue
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
415 data_queue.put(DbData(*data))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
416 else:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
417 sys.stderr.write('PID %s: Error in sibling process detected. Cleaning up.\n' \
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
418 % os.getpid())
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
419 return None
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
420 except:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
421 #sys.stderr.write('Exception in feeder queue feeding step\n')
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
422 alive.value = False
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
423 raise
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
424
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
425 return None
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
426
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
427
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
428 def feedQueueClust(alive, data_queue, db_file, group_func=None, group_args={}):
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
429 """
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
430 Feeds the data queue with Ig records
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
431
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
432 Arguments:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
433 alive = a multiprocessing.Value boolean controlling whether processing continues
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
434 if False exit process
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
435 data_queue = a multiprocessing.Queue to hold data for processing
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
436 db_file = the Ig record database file
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
437
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
438 Returns:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
439 None
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
440 """
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
441 # Open input file and perform grouping
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
442 try:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
443 # Iterate over Ig records and order by junction length
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
444 records = {}
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
445 db_iter = readDbFile(db_file)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
446 for rec in db_iter:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
447 records[rec.id] = rec
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
448 records = OrderedDict(sorted(list(records.items()), key=lambda i: i[1].junction_length))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
449 dist_dict = {}
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
450 for __ in range(len(records)):
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
451 k,v = records.popitem(last=False)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
452 dist_dict[k] = [v].append(list(records.values()))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
453 except:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
454 #sys.stderr.write('Exception in feeder grouping step\n')
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
455 alive.value = False
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
456 raise
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
457
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
458 # Add groups to data queue
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
459 try:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
460 # print 'START FEED', alive.value
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
461 # Iterate over groups and feed data queue
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
462 dist_iter = iter(dist_dict.items())
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
463 while alive.value:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
464 # Get data from queue
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
465 if data_queue.full(): continue
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
466 else: data = next(dist_iter, None)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
467 # Exit upon reaching end of iterator
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
468 if data is None: break
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
469 #print "FEED", alive.value, k
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
470
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
471 # Feed queue
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
472 data_queue.put(DbData(*data))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
473 else:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
474 sys.stderr.write('PID %s: Error in sibling process detected. Cleaning up.\n' \
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
475 % os.getpid())
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
476 return None
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
477 except:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
478 #sys.stderr.write('Exception in feeder queue feeding step\n')
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
479 alive.value = False
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
480 raise
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
481
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
482 return None
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
483
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
484
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
485 def processQueue(alive, data_queue, result_queue, clone_func, clone_args):
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
486 """
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
487 Pulls from data queue, performs calculations, and feeds results queue
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
488
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
489 Arguments:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
490 alive = a multiprocessing.Value boolean controlling whether processing continues
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
491 if False exit process
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
492 data_queue = a multiprocessing.Queue holding data to process
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
493 result_queue = a multiprocessing.Queue to hold processed results
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
494 clone_func = the function to call for clonal assignment
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
495 clone_args = a dictionary of arguments to pass to clone_func
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
496
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
497 Returns:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
498 None
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
499 """
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
500 try:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
501 # Iterator over data queue until sentinel object reached
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
502 while alive.value:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
503 # Get data from queue
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
504 if data_queue.empty(): continue
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
505 else: data = data_queue.get()
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
506 # Exit upon reaching sentinel
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
507 if data is None: break
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
508
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
509 # Define result object for iteration and get data records
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
510 records = data.data
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
511 # print(data.id)
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
512 result = DbResult(data.id, records)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
513
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
514 # Check for invalid data (due to failed indexing) and add failed result
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
515 if not data:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
516 result_queue.put(result)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
517 continue
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
518
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
519 # Add V(D)J to log
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
520 result.log['ID'] = ','.join([str(x) for x in data.id])
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
521 result.log['VALLELE'] = ','.join(set([(r.getVAllele() or '') for r in records]))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
522 result.log['DALLELE'] = ','.join(set([(r.getDAllele() or '') for r in records]))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
523 result.log['JALLELE'] = ','.join(set([(r.getJAllele() or '') for r in records]))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
524 result.log['JUNCLEN'] = ','.join(set([(str(len(r.junction)) or '0') for r in records]))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
525 result.log['SEQUENCES'] = len(records)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
526
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
527 # Checking for preclone failure and assign clones
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
528 clones = clone_func(records, **clone_args) if data else None
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
529
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
530 # import cProfile
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
531 # prof = cProfile.Profile()
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
532 # clones = prof.runcall(clone_func, records, **clone_args)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
533 # prof.dump_stats('worker-%d.prof' % os.getpid())
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
534
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
535 if clones is not None:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
536 result.results = clones
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
537 result.valid = True
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
538 result.log['CLONES'] = len(clones)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
539 else:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
540 result.log['CLONES'] = 0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
541
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
542 # Feed results to result queue
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
543 result_queue.put(result)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
544 else:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
545 sys.stderr.write('PID %s: Error in sibling process detected. Cleaning up.\n' \
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
546 % os.getpid())
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
547 return None
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
548 except:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
549 #sys.stderr.write('Exception in worker\n')
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
550 alive.value = False
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
551 raise
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
552
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
553 return None
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
554
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
555
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
556 def processQueueClust(alive, data_queue, result_queue, clone_func, clone_args):
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
557 """
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
558 Pulls from data queue, performs calculations, and feeds results queue
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
559
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
560 Arguments:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
561 alive = a multiprocessing.Value boolean controlling whether processing continues
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
562 if False exit process
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
563 data_queue = a multiprocessing.Queue holding data to process
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
564 result_queue = a multiprocessing.Queue to hold processed results
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
565 clone_func = the function to call for calculating pairwise distances between sequences
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
566 clone_args = a dictionary of arguments to pass to clone_func
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
567
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
568 Returns:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
569 None
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
570 """
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
571
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
572 try:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
573 # print 'START WORK', alive.value
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
574 # Iterator over data queue until sentinel object reached
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
575 while alive.value:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
576 # Get data from queue
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
577 if data_queue.empty(): continue
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
578 else: data = data_queue.get()
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
579 # Exit upon reaching sentinel
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
580 if data is None: break
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
581 # print "WORK", alive.value, data['id']
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
582
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
583 # Define result object for iteration and get data records
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
584 records = data.data
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
585 result = DbResult(data.id, records)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
586
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
587 # Create row of distance matrix and check for error
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
588 dist_row = clone_func(records, **clone_args) if data else None
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
589 if dist_row is not None:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
590 result.results = dist_row
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
591 result.valid = True
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
592
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
593 # Feed results to result queue
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
594 result_queue.put(result)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
595 else:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
596 sys.stderr.write('PID %s: Error in sibling process detected. Cleaning up.\n' \
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
597 % os.getpid())
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
598 return None
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
599 except:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
600 #sys.stderr.write('Exception in worker\n')
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
601 alive.value = False
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
602 raise
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
603
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
604 return None
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
605
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
606
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
607 def collectQueue(alive, result_queue, collect_queue, db_file, out_args, cluster_func=None, cluster_args={}):
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
608 """
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
609 Assembles results from a queue of individual sequence results and manages log/file I/O
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
610
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
611 Arguments:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
612 alive = a multiprocessing.Value boolean controlling whether processing continues
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
613 if False exit process
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
614 result_queue = a multiprocessing.Queue holding processQueue results
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
615 collect_queue = a multiprocessing.Queue to store collector return values
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
616 db_file = the input database file name
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
617 out_args = common output argument dictionary from parseCommonArgs
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
618 cluster_func = the function to call for carrying out clustering on distance matrix
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
619 cluster_args = a dictionary of arguments to pass to cluster_func
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
620
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
621 Returns:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
622 None
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
623 (adds 'log' and 'out_files' to collect_dict)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
624 """
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
625 # Open output files
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
626 try:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
627 # Count records and define output format
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
628 out_type = getFileType(db_file) if out_args['out_type'] is None \
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
629 else out_args['out_type']
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
630 result_count = countDbFile(db_file)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
631
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
632 # Defined successful output handle
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
633 pass_handle = getOutputHandle(db_file,
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
634 out_label='clone-pass',
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
635 out_dir=out_args['out_dir'],
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
636 out_name=out_args['out_name'],
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
637 out_type=out_type)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
638 pass_writer = getDbWriter(pass_handle, db_file, add_fields='CLONE')
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
639
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
640 # Defined failed alignment output handle
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
641 if out_args['failed']:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
642 fail_handle = getOutputHandle(db_file,
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
643 out_label='clone-fail',
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
644 out_dir=out_args['out_dir'],
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
645 out_name=out_args['out_name'],
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
646 out_type=out_type)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
647 fail_writer = getDbWriter(fail_handle, db_file)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
648 else:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
649 fail_handle = None
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
650 fail_writer = None
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
651
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
652 # Define log handle
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
653 if out_args['log_file'] is None:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
654 log_handle = None
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
655 else:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
656 log_handle = open(out_args['log_file'], 'w')
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
657 except:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
658 #sys.stderr.write('Exception in collector file opening step\n')
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
659 alive.value = False
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
660 raise
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
661
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
662 # Get results from queue and write to files
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
663 try:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
664 #print 'START COLLECT', alive.value
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
665 # Iterator over results queue until sentinel object reached
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
666 start_time = time()
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
667 rec_count = clone_count = pass_count = fail_count = 0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
668 while alive.value:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
669 # Get result from queue
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
670 if result_queue.empty(): continue
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
671 else: result = result_queue.get()
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
672 # Exit upon reaching sentinel
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
673 if result is None: break
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
674 #print "COLLECT", alive.value, result['id']
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
675
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
676 # Print progress for previous iteration and update record count
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
677 if rec_count == 0:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
678 print('PROGRESS> Assigning clones')
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
679 printProgress(rec_count, result_count, 0.05, start_time)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
680 rec_count += len(result.data)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
681
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
682 # Write passed and failed records
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
683 if result:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
684 for clone in result.results.values():
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
685 clone_count += 1
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
686 for i, rec in enumerate(clone):
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
687 rec.annotations['CLONE'] = clone_count
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
688 pass_writer.writerow(rec.toDict())
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
689 pass_count += 1
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
690 result.log['CLONE%i-%i' % (clone_count, i + 1)] = str(rec.junction)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
691
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
692 else:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
693 for i, rec in enumerate(result.data):
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
694 if fail_writer is not None: fail_writer.writerow(rec.toDict())
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
695 fail_count += 1
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
696 result.log['CLONE0-%i' % (i + 1)] = str(rec.junction)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
697
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
698 # Write log
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
699 printLog(result.log, handle=log_handle)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
700 else:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
701 sys.stderr.write('PID %s: Error in sibling process detected. Cleaning up.\n' \
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
702 % os.getpid())
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
703 return None
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
704
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
705 # Print total counts
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
706 printProgress(rec_count, result_count, 0.05, start_time)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
707
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
708 # Close file handles
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
709 pass_handle.close()
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
710 if fail_handle is not None: fail_handle.close()
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
711 if log_handle is not None: log_handle.close()
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
712
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
713 # Update return list
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
714 log = OrderedDict()
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
715 log['OUTPUT'] = os.path.basename(pass_handle.name)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
716 log['CLONES'] = clone_count
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
717 log['RECORDS'] = rec_count
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
718 log['PASS'] = pass_count
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
719 log['FAIL'] = fail_count
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
720 collect_dict = {'log':log, 'out_files': [pass_handle.name]}
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
721 collect_queue.put(collect_dict)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
722 except:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
723 #sys.stderr.write('Exception in collector result processing step\n')
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
724 alive.value = False
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
725 raise
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
726
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
727 return None
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
728
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
729
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
730 def collectQueueClust(alive, result_queue, collect_queue, db_file, out_args, cluster_func, cluster_args):
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
731 """
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
732 Assembles results from a queue of individual sequence results and manages log/file I/O
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
733
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
734 Arguments:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
735 alive = a multiprocessing.Value boolean controlling whether processing continues
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
736 if False exit process
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
737 result_queue = a multiprocessing.Queue holding processQueue results
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
738 collect_queue = a multiprocessing.Queue to store collector return values
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
739 db_file = the input database file name
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
740 out_args = common output argument dictionary from parseCommonArgs
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
741 cluster_func = the function to call for carrying out clustering on distance matrix
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
742 cluster_args = a dictionary of arguments to pass to cluster_func
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
743
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
744 Returns:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
745 None
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
746 (adds 'log' and 'out_files' to collect_dict)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
747 """
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
748 # Open output files
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
749 try:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
750
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
751 # Iterate over Ig records to count and order by junction length
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
752 result_count = 0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
753 records = {}
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
754 # print 'Reading file...'
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
755 db_iter = readDbFile(db_file)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
756 for rec in db_iter:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
757 records[rec.id] = rec
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
758 result_count += 1
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
759 records = OrderedDict(sorted(list(records.items()), key=lambda i: i[1].junction_length))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
760
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
761 # Define empty matrix to store assembled results
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
762 dist_mat = np.zeros((result_count,result_count))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
763
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
764 # Count records and define output format
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
765 out_type = getFileType(db_file) if out_args['out_type'] is None \
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
766 else out_args['out_type']
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
767
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
768 # Defined successful output handle
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
769 pass_handle = getOutputHandle(db_file,
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
770 out_label='clone-pass',
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
771 out_dir=out_args['out_dir'],
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
772 out_name=out_args['out_name'],
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
773 out_type=out_type)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
774 pass_writer = getDbWriter(pass_handle, db_file, add_fields='CLONE')
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
775
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
776 # Defined failed cloning output handle
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
777 if out_args['failed']:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
778 fail_handle = getOutputHandle(db_file,
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
779 out_label='clone-fail',
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
780 out_dir=out_args['out_dir'],
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
781 out_name=out_args['out_name'],
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
782 out_type=out_type)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
783 fail_writer = getDbWriter(fail_handle, db_file)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
784 else:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
785 fail_handle = None
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
786 fail_writer = None
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
787
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
788 # Open log file
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
789 if out_args['log_file'] is None:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
790 log_handle = None
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
791 else:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
792 log_handle = open(out_args['log_file'], 'w')
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
793 except:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
794 alive.value = False
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
795 raise
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
796
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
797 try:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
798 # Iterator over results queue until sentinel object reached
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
799 start_time = time()
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
800 row_count = rec_count = 0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
801 while alive.value:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
802 # Get result from queue
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
803 if result_queue.empty(): continue
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
804 else: result = result_queue.get()
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
805 # Exit upon reaching sentinel
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
806 if result is None: break
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
807
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
808 # Print progress for previous iteration
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
809 if row_count == 0:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
810 print('PROGRESS> Assigning clones')
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
811 printProgress(row_count, result_count, 0.05, start_time)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
812
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
813 # Update counts for iteration
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
814 row_count += 1
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
815 rec_count += len(result)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
816
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
817 # Add result row to distance matrix
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
818 if result:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
819 dist_mat[list(range(result_count-len(result),result_count)),result_count-len(result)] = result.results
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
820
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
821 else:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
822 sys.stderr.write('PID %s: Error in sibling process detected. Cleaning up.\n' \
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
823 % os.getpid())
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
824 return None
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
825
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
826 # Calculate linkage and carry out clustering
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
827 # print dist_mat
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
828 clusters = cluster_func(dist_mat, **cluster_args) if dist_mat is not None else None
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
829 clones = {}
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
830 # print clusters
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
831 for i, c in enumerate(clusters):
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
832 clones.setdefault(c, []).append(records[list(records.keys())[i]])
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
833
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
834 # Write passed and failed records
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
835 clone_count = pass_count = fail_count = 0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
836 if clones:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
837 for clone in clones.values():
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
838 clone_count += 1
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
839 for i, rec in enumerate(clone):
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
840 rec.annotations['CLONE'] = clone_count
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
841 pass_writer.writerow(rec.toDict())
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
842 pass_count += 1
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
843 #result.log['CLONE%i-%i' % (clone_count, i + 1)] = str(rec.junction)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
844
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
845 else:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
846 for i, rec in enumerate(result.data):
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
847 fail_writer.writerow(rec.toDict())
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
848 fail_count += 1
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
849 #result.log['CLONE0-%i' % (i + 1)] = str(rec.junction)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
850
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
851 # Print final progress
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
852 printProgress(row_count, result_count, 0.05, start_time)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
853
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
854 # Close file handles
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
855 pass_handle.close()
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
856 if fail_handle is not None: fail_handle.close()
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
857 if log_handle is not None: log_handle.close()
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
858
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
859 # Update return list
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
860 log = OrderedDict()
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
861 log['OUTPUT'] = os.path.basename(pass_handle.name)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
862 log['CLONES'] = clone_count
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
863 log['RECORDS'] = rec_count
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
864 log['PASS'] = pass_count
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
865 log['FAIL'] = fail_count
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
866 collect_dict = {'log':log, 'out_files': [pass_handle.name]}
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
867 collect_queue.put(collect_dict)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
868 except:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
869 alive.value = False
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
870 raise
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
871
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
872 return None
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
873
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
874
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
875 def defineClones(db_file, feed_func, work_func, collect_func, clone_func, cluster_func=None,
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
876 group_func=None, group_args={}, clone_args={}, cluster_args={},
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
877 out_args=default_out_args, nproc=None, queue_size=None):
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
878 """
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
879 Define clonally related sequences
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
880
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
881 Arguments:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
882 db_file = filename of input database
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
883 feed_func = the function that feeds the queue
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
884 work_func = the worker function that will run on each CPU
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
885 collect_func = the function that collects results from the workers
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
886 group_func = the function to use for assigning preclones
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
887 clone_func = the function to use for determining clones within preclonal groups
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
888 group_args = a dictionary of arguments to pass to group_func
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
889 clone_args = a dictionary of arguments to pass to clone_func
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
890 out_args = common output argument dictionary from parseCommonArgs
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
891 nproc = the number of processQueue processes;
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
892 if None defaults to the number of CPUs
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
893 queue_size = maximum size of the argument queue;
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
894 if None defaults to 2*nproc
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
895
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
896 Returns:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
897 a list of successful output file names
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
898 """
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
899 # Print parameter info
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
900 log = OrderedDict()
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
901 log['START'] = 'DefineClones'
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
902 log['DB_FILE'] = os.path.basename(db_file)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
903 if group_func is not None:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
904 log['GROUP_FUNC'] = group_func.__name__
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
905 log['GROUP_ARGS'] = group_args
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
906 log['CLONE_FUNC'] = clone_func.__name__
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
907
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
908 # TODO: this is yucky, but can be fixed by using a model class
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
909 clone_log = clone_args.copy()
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
910 if 'dist_mat' in clone_log: del clone_log['dist_mat']
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
911 log['CLONE_ARGS'] = clone_log
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
912
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
913 if cluster_func is not None:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
914 log['CLUSTER_FUNC'] = cluster_func.__name__
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
915 log['CLUSTER_ARGS'] = cluster_args
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
916 log['NPROC'] = nproc
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
917 printLog(log)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
918
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
919 # Define feeder function and arguments
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
920 feed_args = {'db_file': db_file,
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
921 'group_func': group_func,
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
922 'group_args': group_args}
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
923 # Define worker function and arguments
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
924 work_args = {'clone_func': clone_func,
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
925 'clone_args': clone_args}
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
926 # Define collector function and arguments
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
927 collect_args = {'db_file': db_file,
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
928 'out_args': out_args,
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
929 'cluster_func': cluster_func,
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
930 'cluster_args': cluster_args}
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
931
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
932 # Call process manager
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
933 result = manageProcesses(feed_func, work_func, collect_func,
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
934 feed_args, work_args, collect_args,
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
935 nproc, queue_size)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
936
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
937 # Print log
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
938 result['log']['END'] = 'DefineClones'
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
939 printLog(result['log'])
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
940
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
941 return result['out_files']
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
942
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
943
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
944 def getArgParser():
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
945 """
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
946 Defines the ArgumentParser
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
947
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
948 Arguments:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
949 None
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
950
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
951 Returns:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
952 an ArgumentParser object
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
953 """
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
954 # Define input and output fields
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
955 fields = dedent(
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
956 '''
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
957 output files:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
958 clone-pass
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
959 database with assigned clonal group numbers.
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
960 clone-fail
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
961 database with records failing clonal grouping.
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
962
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
963 required fields:
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
964 SEQUENCE_ID, V_CALL or V_CALL_GENOTYPED, D_CALL, J_CALL, JUNCTION
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
965
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
966 <field>
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
967 sequence field specified by the --sf parameter
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
968
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
969 output fields:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
970 CLONE
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
971 ''')
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
972
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
973 # Define ArgumentParser
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
974 parser = ArgumentParser(description=__doc__, epilog=fields,
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
975 formatter_class=CommonHelpFormatter)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
976 parser.add_argument('--version', action='version',
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
977 version='%(prog)s:' + ' %s-%s' %(__version__, __date__))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
978 subparsers = parser.add_subparsers(title='subcommands', dest='command', metavar='',
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
979 help='Cloning method')
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
980 # TODO: This is a temporary fix for Python issue 9253
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
981 subparsers.required = True
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
982
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
983 # Parent parser
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
984 parser_parent = getCommonArgParser(seq_in=False, seq_out=False, db_in=True,
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
985 multiproc=True)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
986
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
987 # Distance cloning method
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
988 parser_bygroup = subparsers.add_parser('bygroup', parents=[parser_parent],
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
989 formatter_class=CommonHelpFormatter,
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
990 help='''Defines clones as having same V assignment,
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
991 J assignment, and junction length with
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
992 specified substitution distance model.''',
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
993 description='''Defines clones as having same V assignment,
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
994 J assignment, and junction length with
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
995 specified substitution distance model.''')
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
996 parser_bygroup.add_argument('-f', nargs='+', action='store', dest='fields', default=None,
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
997 help='Additional fields to use for grouping clones (non VDJ)')
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
998 parser_bygroup.add_argument('--mode', action='store', dest='mode',
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
999 choices=('allele', 'gene'), default=default_index_mode,
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1000 help='''Specifies whether to use the V(D)J allele or gene for
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1001 initial grouping.''')
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1002 parser_bygroup.add_argument('--act', action='store', dest='action',
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1003 choices=('first', 'set'), default=default_index_action,
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1004 help='''Specifies how to handle multiple V(D)J assignments
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1005 for initial grouping.''')
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1006 parser_bygroup.add_argument('--model', action='store', dest='model',
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1007 choices=choices_bygroup_model,
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1008 default=default_bygroup_model,
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1009 help='''Specifies which substitution model to use for calculating distance
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1010 between sequences. The "ham" model is nucleotide Hamming distance and
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1011 "aa" is amino acid Hamming distance. The "hh_s1f" and "hh_s5f" models are
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1012 human specific single nucleotide and 5-mer content models, respectively,
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1013 from Yaari et al, 2013. The "mk_rs1nf" and "mk_rs5nf" models are
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1014 mouse specific single nucleotide and 5-mer content models, respectively,
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1015 from Cui et al, 2016. The "m1n_compat" and "hs1f_compat" models are
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1016 deprecated models provided backwards compatibility with the "m1n" and
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1017 "hs1f" models in Change-O v0.3.3 and SHazaM v0.1.4. Both
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1018 5-mer models should be considered experimental.''')
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1019 parser_bygroup.add_argument('--dist', action='store', dest='distance', type=float,
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1020 default=default_distance,
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1021 help='The distance threshold for clonal grouping')
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1022 parser_bygroup.add_argument('--norm', action='store', dest='norm',
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1023 choices=('len', 'mut', 'none'), default=default_norm,
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1024 help='''Specifies how to normalize distances. One of none
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1025 (do not normalize), len (normalize by length),
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1026 or mut (normalize by number of mutations between sequences).''')
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1027 parser_bygroup.add_argument('--sym', action='store', dest='sym',
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1028 choices=('avg', 'min'), default=default_sym,
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1029 help='''Specifies how to combine asymmetric distances. One of avg
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1030 (average of A->B and B->A) or min (minimum of A->B and B->A).''')
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1031 parser_bygroup.add_argument('--link', action='store', dest='linkage',
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1032 choices=('single', 'average', 'complete'), default=default_linkage,
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1033 help='''Type of linkage to use for hierarchical clustering.''')
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1034 parser_bygroup.add_argument('--sf', action='store', dest='seq_field',
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1035 default=default_seq_field,
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1036 help='''The name of the field to be used to calculate
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1037 distance between records''')
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1038 parser_bygroup.set_defaults(feed_func=feedQueue)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1039 parser_bygroup.set_defaults(work_func=processQueue)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1040 parser_bygroup.set_defaults(collect_func=collectQueue)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1041 parser_bygroup.set_defaults(group_func=indexJunctions)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1042 parser_bygroup.set_defaults(clone_func=distanceClones)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1043
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1044 # Chen2010
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1045 parser_chen = subparsers.add_parser('chen2010', parents=[parser_parent],
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1046 formatter_class=CommonHelpFormatter,
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1047 help='''Defines clones by method specified in Chen, 2010.''',
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1048 description='''Defines clones by method specified in Chen, 2010.''')
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1049 parser_chen.set_defaults(feed_func=feedQueueClust)
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1050 parser_chen.set_defaults(work_func=processQueueClust)
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1051 parser_chen.set_defaults(collect_func=collectQueueClust)
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1052 parser_chen.set_defaults(cluster_func=hierClust)
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1053
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1054 # Ademokun2011
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1055 parser_ade = subparsers.add_parser('ademokun2011', parents=[parser_parent],
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1056 formatter_class=CommonHelpFormatter,
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1057 help='''Defines clones by method specified in Ademokun, 2011.''',
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1058 description='''Defines clones by method specified in Ademokun, 2011.''')
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1059 parser_ade.set_defaults(feed_func=feedQueueClust)
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1060 parser_ade.set_defaults(work_func=processQueueClust)
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1061 parser_ade.set_defaults(collect_func=collectQueueClust)
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1062 parser_ade.set_defaults(cluster_func=hierClust)
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1063
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1064 return parser
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1065
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1066
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1067 if __name__ == '__main__':
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1068 """
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1069 Parses command line arguments and calls main function
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1070 """
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1071 # Parse arguments
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1072 parser = getArgParser()
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1073 args = parser.parse_args()
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1074 args_dict = parseCommonArgs(args)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1075 # Convert case of fields
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1076 if 'seq_field' in args_dict:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1077 args_dict['seq_field'] = args_dict['seq_field'].upper()
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1078 if 'fields' in args_dict and args_dict['fields'] is not None:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1079 args_dict['fields'] = [f.upper() for f in args_dict['fields']]
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1080
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1081 # Define clone_args
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1082 if args.command == 'bygroup':
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1083 args_dict['group_args'] = {'fields': args_dict['fields'],
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1084 'action': args_dict['action'],
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1085 'mode':args_dict['mode']}
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1086 args_dict['clone_args'] = {'model': args_dict['model'],
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1087 'distance': args_dict['distance'],
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1088 'norm': args_dict['norm'],
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1089 'sym': args_dict['sym'],
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1090 'linkage': args_dict['linkage'],
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1091 'seq_field': args_dict['seq_field']}
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1092
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1093 # Get distance matrix
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1094 try:
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1095 args_dict['clone_args']['dist_mat'] = distance_models[args_dict['model']]
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1096 except KeyError:
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1097 sys.exit('Unrecognized distance model: %s' % args_dict['model'])
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1098
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1099 del args_dict['fields']
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1100 del args_dict['action']
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1101 del args_dict['mode']
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1102 del args_dict['model']
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1103 del args_dict['distance']
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1104 del args_dict['norm']
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1105 del args_dict['sym']
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1106 del args_dict['linkage']
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1107 del args_dict['seq_field']
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1108
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1109 # Define clone_args
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1110 if args.command == 'chen2010':
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1111 args_dict['clone_func'] = distChen2010
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1112 args_dict['cluster_args'] = {'method': args.command }
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1113
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1114 if args.command == 'ademokun2011':
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1115 args_dict['clone_func'] = distAdemokun2011
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1116 args_dict['cluster_args'] = {'method': args.command }
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1117
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1118 # Call defineClones
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1119 del args_dict['command']
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1120 del args_dict['db_files']
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1121 for f in args.__dict__['db_files']:
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1122 args_dict['db_file'] = f
52
22dddabe3637 Uploaded
davidvanzessen
parents: 0
diff changeset
1123 defineClones(**args_dict)