Mercurial > repos > davidvanzessen > change_o

comparison CreateGermlines.py @ 0:dda9b2e72e2b draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Tue, 03 May 2016 09:52:21 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:dda9b2e72e2b
1 #!/usr/bin/env python3
2 """
3 Reconstructs germline sequences from alignment data
4 """
5 # Info
6 __author__ = 'Namita Gupta, Jason Anthony Vander Heiden'
7 from changeo import __version__, __date__
8
9 # Imports
10 import os
11 import sys
12 from argparse import ArgumentParser
13 from collections import OrderedDict
14 from textwrap import dedent
15 from time import time
16
17 # Presto and change imports
18 from presto.Defaults import default_out_args
19 from presto.IO import getOutputHandle, printLog, printProgress
20 from changeo.Commandline import CommonHelpFormatter, getCommonArgParser, parseCommonArgs
21 from changeo.IO import getDbWriter, readDbFile, countDbFile, getRepo
22 from changeo.Receptor import allele_regex, parseAllele
23
24 # Defaults
25 default_germ_types = 'dmask'
26 default_v_field = 'V_CALL'
27 default_seq_field = 'SEQUENCE_IMGT'
28
29
30 def joinGermline(align, repo_dict, germ_types, v_field, seq_field):
31 """
32 Join gapped germline sequences aligned with sample sequences
33
34 Arguments:
35 align = iterable yielding dictionaries of sample sequence data
36 repo_dict = dictionary of IMGT gapped germline sequences
37 germ_types = types of germline sequences to be output
38 (full germline, D-region masked, only V-region germline)
39 v_field = field in which to look for V call
40 seq_field = field in which to look for sequence
41
42 Returns:
43 dictionary of germline_type: germline_sequence
44 """
45 j_field = 'J_CALL'
46 germlines = {'full': '', 'dmask': '', 'vonly': ''}
47 result_log = OrderedDict()
48 result_log['ID'] = align['SEQUENCE_ID']
49
50 # Find germline V-region gene
51 if v_field == 'V_CALL_GENOTYPED':
52 vgene = parseAllele(align[v_field], allele_regex, 'list')
53 vkey = vgene
54 else:
55 vgene = parseAllele(align[v_field], allele_regex, 'first')
56 vkey = (vgene, )
57
58 # Build V-region germline
59 if vgene is not None:
60 result_log['V_CALL'] = ','.join(vkey)
61 if vkey in repo_dict:
62 vseq = repo_dict[vkey]
63 # Germline start
64 try: vstart = int(align['V_GERM_START_IMGT']) - 1
65 except (TypeError, ValueError): vstart = 0
66 # Germline length
67 try: vlen = int(align['V_GERM_LENGTH_IMGT'])
68 except (TypeError, ValueError): vlen = 0
69 # TODO: not sure what this line is doing here. it no make no sense.
70 vpad = vlen - len(vseq[vstart:])
71 if vpad < 0: vpad = 0
72 germ_vseq = vseq[vstart:(vstart + vlen)] + ('N' * vpad)
73 else:
74 result_log['ERROR'] = 'Germline %s not in repertoire' % ','.join(vkey)
75 return result_log, germlines
76 else:
77 result_log['V_CALL'] = None
78 try: vlen = int(align['V_GERM_LENGTH_IMGT'])
79 except (TypeError, ValueError): vlen = 0
80 germ_vseq = 'N' * vlen
81
82 # Find germline D-region gene
83 dgene = parseAllele(align['D_CALL'], allele_regex, 'first')
84
85 # Build D-region germline
86 if dgene is not None:
87 result_log['D_CALL'] = dgene
88 dkey = (dgene, )
89 if dkey in repo_dict:
90 dseq = repo_dict[dkey]
91 # Germline start
92 try: dstart = int(align['D_GERM_START']) - 1
93 except (TypeError, ValueError): dstart = 0
94 # Germline length
95 try: dlen = int(align['D_GERM_LENGTH'])
96 except (TypeError, ValueError): dlen = 0
97 germ_dseq = repo_dict[dkey][dstart:(dstart + dlen)]
98 else:
99 result_log['ERROR'] = 'Germline %s not in repertoire' % dgene
100 return result_log, germlines
101 else:
102 result_log['D_CALL'] = None
103 germ_dseq = ''
104
105 # Find germline J-region gene
106 jgene = parseAllele(align[j_field], allele_regex, 'first')
107
108 # Build D-region germline
109 if jgene is not None:
110 result_log['J_CALL'] = jgene
111 jkey = (jgene, )
112 if jkey in repo_dict:
113 jseq = repo_dict[jkey]
114 # Germline start
115 try: jstart = int(align['J_GERM_START']) - 1
116 except (TypeError, ValueError): jstart = 0
117 # Germline length
118 try: jlen = int(align['J_GERM_LENGTH'])
119 except (TypeError, ValueError): jlen = 0
120 # TODO: not sure what this line is doing either
121 jpad = jlen - len(jseq[jstart:])
122 if jpad < 0: jpad = 0
123 germ_jseq = jseq[jstart:(jstart + jlen)] + ('N' * jpad)
124 else:
125 result_log['ERROR'] = 'Germline %s not in repertoire' % jgene
126 return result_log, germlines
127 else:
128 result_log['J_CALL'] = None
129 try: jlen = int(align['J_GERM_LENGTH'])
130 except (TypeError, ValueError): jlen = 0
131 germ_jseq = 'N' * jlen
132
133 # Assemble pieces starting with V-region
134 germ_seq = germ_vseq
135 regions = 'V' * len(germ_vseq)
136
137 # Nucleotide additions before D (before J for light chains)
138 try: n1_len = int(align['N1_LENGTH'])
139 except (TypeError, ValueError): n1_len = 0
140 if n1_len < 0:
141 result_log['ERROR'] = 'N1_LENGTH is negative'
142 return result_log, germlines
143
144 germ_seq += 'N' * n1_len
145 regions += 'N' * n1_len
146
147 # Add D-region
148 germ_seq += germ_dseq
149 regions += 'D' * len(germ_dseq)
150 #print 'VD>', germ_seq, '\nVD>', regions
151
152 # Nucleotide additions after D (heavy chains only)
153 try: n2_len = int(align['N2_LENGTH'])
154 except (TypeError, ValueError): n2_len = 0
155 if n2_len < 0:
156 result_log['ERROR'] = 'N2_LENGTH is negative'
157 return result_log, germlines
158
159 germ_seq += 'N' * n2_len
160 regions += 'N' * n2_len
161
162 # Add J-region
163 germ_seq += germ_jseq
164 regions += 'J' * len(germ_jseq)
165
166 # Define return germlines
167 germlines['full'] = germ_seq
168 germlines['regions'] = regions
169 if 'dmask' in germ_types:
170 germlines['dmask'] = germ_seq[:len(germ_vseq)] + \
171 'N' * (len(germ_seq) - len(germ_vseq) - len(germ_jseq)) + \
172 germ_seq[-len(germ_jseq):]
173 if 'vonly' in germ_types:
174 germlines['vonly'] = germ_vseq
175
176 # Check that input and germline sequence match
177 if len(align[seq_field]) == 0:
178 result_log['ERROR'] = 'Sequence is missing from %s column' % seq_field
179 elif len(germlines['full']) != len(align[seq_field]):
180 result_log['ERROR'] = 'Germline sequence is %d nucleotides longer than input sequence' % \
181 (len(germlines['full']) - len(align[seq_field]))
182
183 # Convert to uppercase
184 for k, v in germlines.items(): germlines[k] = v.upper()
185
186 return result_log, germlines
187
188
189 def assembleEachGermline(db_file, repo, germ_types, v_field, seq_field, out_args=default_out_args):
190 """
191 Write germline sequences to tab-delimited database file
192
193 Arguments:
194 db_file = input tab-delimited database file
195 repo = folder with germline repertoire files
196 germ_types = types of germline sequences to be output
197 (full germline, D-region masked, only V-region germline)
198 v_field = field in which to look for V call
199 seq_field = field in which to look for sequence
200 out_args = arguments for output preferences
201
202 Returns:
203 None
204 """
205 # Print parameter info
206 log = OrderedDict()
207 log['START'] = 'CreateGermlines'
208 log['DB_FILE'] = os.path.basename(db_file)
209 log['GERM_TYPES'] = germ_types if isinstance(germ_types, str) else ','.join(germ_types)
210 log['CLONED'] = 'False'
211 log['V_FIELD'] = v_field
212 log['SEQ_FIELD'] = seq_field
213 printLog(log)
214
215 # Get repertoire and open Db reader
216 repo_dict = getRepo(repo)
217 reader = readDbFile(db_file, ig=False)
218
219 # Exit if V call field does not exist in reader
220 if v_field not in reader.fieldnames:
221 sys.exit('Error: V field does not exist in input database file.')
222
223 # Define log handle
224 if out_args['log_file'] is None:
225 log_handle = None
226 else:
227 log_handle = open(out_args['log_file'], 'w')
228
229 add_fields = []
230 seq_type = seq_field.split('_')[-1]
231 if 'full' in germ_types: add_fields += ['GERMLINE_' + seq_type]
232 if 'dmask' in germ_types: add_fields += ['GERMLINE_' + seq_type + '_D_MASK']
233 if 'vonly' in germ_types: add_fields += ['GERMLINE_' + seq_type + '_V_REGION']
234
235 # Create output file handle and Db writer
236 pass_handle = getOutputHandle(db_file, 'germ-pass',
237 out_dir=out_args['out_dir'],
238 out_name=out_args['out_name'],
239 out_type=out_args['out_type'])
240 pass_writer = getDbWriter(pass_handle, db_file, add_fields=add_fields)
241
242 if out_args['failed']:
243 fail_handle = getOutputHandle(db_file, 'germ-fail',
244 out_dir=out_args['out_dir'],
245 out_name=out_args['out_name'],
246 out_type=out_args['out_type'])
247 fail_writer = getDbWriter(fail_handle, db_file, add_fields=add_fields)
248 else:
249 fail_handle = None
250 fail_writer = None
251
252 # Initialize time and total count for progress bar
253 start_time = time()
254 rec_count = countDbFile(db_file)
255 pass_count = fail_count = 0
256 # Iterate over rows
257 for i,row in enumerate(reader):
258 # Print progress
259 printProgress(i, rec_count, 0.05, start_time)
260
261 result_log, germlines = joinGermline(row, repo_dict, germ_types, v_field, seq_field)
262
263 # Add germline field(s) to dictionary
264 if 'full' in germ_types: row['GERMLINE_' + seq_type] = germlines['full']
265 if 'dmask' in germ_types: row['GERMLINE_' + seq_type + '_D_MASK'] = germlines['dmask']
266 if 'vonly' in germ_types: row['GERMLINE_' + seq_type + '_V_REGION'] = germlines['vonly']
267
268 # Write row to pass or fail file
269 if 'ERROR' in result_log:
270 fail_count += 1
271 if fail_writer is not None: fail_writer.writerow(row)
272 else:
273 result_log['SEQUENCE'] = row[seq_field]
274 result_log['GERMLINE'] = germlines['full']
275 result_log['REGIONS'] = germlines['regions']
276
277 pass_count += 1
278 pass_writer.writerow(row)
279 printLog(result_log, handle=log_handle)
280
281 # Print log
282 printProgress(i+1, rec_count, 0.05, start_time)
283 log = OrderedDict()
284 log['OUTPUT'] = os.path.basename(pass_handle.name)
285 log['RECORDS'] = rec_count
286 log['PASS'] = pass_count
287 log['FAIL'] = fail_count
288 log['END'] = 'CreateGermlines'
289 printLog(log)
290
291 # Close file handles
292 pass_handle.close()
293 if fail_handle is not None: fail_handle.close()
294 if log_handle is not None: log_handle.close()
295
296
297 def makeCloneGermline(clone, clone_dict, repo_dict, germ_types, v_field, seq_field, counts, writers, out_args):
298 """
299 Determine consensus clone sequence and create germline for clone
300
301 Arguments:
302 clone = clone ID
303 clone_dict = iterable yielding dictionaries of sequence data from clone
304 repo_dict = dictionary of IMGT gapped germline sequences
305 germ_types = types of germline sequences to be output
306 (full germline, D-region masked, only V-region germline)
307 v_field = field in which to look for V call
308 seq_field = field in which to look for sequence
309 counts = dictionary of pass counter and fail counter
310 writers = dictionary with pass and fail DB writers
311 out_args = arguments for output preferences
312
313 Returns:
314 None
315 """
316 seq_type = seq_field.split('_')[-1]
317 j_field = 'J_CALL'
318
319 # Create dictionaries to count observed V/J calls
320 v_dict = OrderedDict()
321 j_dict = OrderedDict()
322
323 # Find longest sequence in clone
324 max_length = 0
325 for val in clone_dict.values():
326 v = val[v_field]
327 v_dict[v] = v_dict.get(v,0) + 1
328 j = val[j_field]
329 j_dict[j] = j_dict.get(j,0) + 1
330 if len(val[seq_field]) > max_length: max_length = len(val[seq_field])
331
332 # Consensus V and J having most observations
333 v_cons = [k for k in list(v_dict.keys()) if v_dict[k] == max(v_dict.values())]
334 j_cons = [k for k in list(j_dict.keys()) if j_dict[k] == max(j_dict.values())]
335 # Consensus sequence(s) with consensus V/J calls and longest sequence
336 cons = [val for val in list(clone_dict.values()) if val.get(v_field,'') in v_cons and \
337 val.get(j_field,'') in j_cons and \
338 len(val[seq_field])==max_length]
339 # Sequence(s) with consensus V/J are not longest
340 if not cons:
341 # Sequence(s) with consensus V/J (not longest)
342 cons = [val for val in list(clone_dict.values()) if val.get(v_field,'') in v_cons and val.get(j_field,'') in j_cons]
343
344 # No sequence has both consensus V and J call
345 if not cons:
346 result_log = OrderedDict()
347 result_log['ID'] = clone
348 result_log['V_CALL'] = ','.join(v_cons)
349 result_log['J_CALL'] = ','.join(j_cons)
350 result_log['ERROR'] = 'No consensus sequence for clone found'
351 else:
352 # Pad end of consensus sequence with gaps to make it the max length
353 cons = cons[0]
354 cons['J_GERM_LENGTH'] = str(int(cons['J_GERM_LENGTH'] or 0) + max_length - len(cons[seq_field]))
355 cons[seq_field] += '.'*(max_length - len(cons[seq_field]))
356 result_log, germlines = joinGermline(cons, repo_dict, germ_types, v_field, seq_field)
357 result_log['ID'] = clone
358 result_log['CONSENSUS'] = cons['SEQUENCE_ID']
359 else:
360 cons = cons[0]
361 result_log, germlines = joinGermline(cons, repo_dict, germ_types, v_field, seq_field)
362 result_log['ID'] = clone
363 result_log['CONSENSUS'] = cons['SEQUENCE_ID']
364
365 # Write sequences of clone
366 for val in clone_dict.values():
367 if 'ERROR' not in result_log:
368 # Update lengths padded to longest sequence in clone
369 val['J_GERM_LENGTH'] = str(int(val['J_GERM_LENGTH'] or 0) + max_length - len(val[seq_field]))
370 val[seq_field] += '.'*(max_length - len(val[seq_field]))
371
372 # Add column(s) to tab-delimited database file
373 if 'full' in germ_types: val['GERMLINE_' + seq_type] = germlines['full']
374 if 'dmask' in germ_types: val['GERMLINE_' + seq_type + '_D_MASK'] = germlines['dmask']
375 if 'vonly' in germ_types: val['GERMLINE_' + seq_type + '_V_REGION'] = germlines['vonly']
376
377 result_log['SEQUENCE'] = cons[seq_field]
378 result_log['GERMLINE'] = germlines['full']
379 result_log['REGIONS'] = germlines['regions']
380
381 # Write to pass file
382 counts['pass'] += 1
383 writers['pass'].writerow(val)
384 else:
385 # Write to fail file
386 counts['fail'] += 1
387 if writers['fail'] is not None: writers['fail'].writerow(val)
388 # Return log
389 return result_log
390
391
392 def assembleCloneGermline(db_file, repo, germ_types, v_field, seq_field, out_args=default_out_args):
393 """
394 Assemble one germline sequence for each clone in a tab-delimited database file
395
396 Arguments:
397 db_file = input tab-delimited database file
398 repo = folder with germline repertoire files
399 germ_types = types of germline sequences to be output
400 (full germline, D-region masked, only V-region germline)
401 v_field = field in which to look for V call
402 seq_field = field in which to look for sequence
403 out_args = arguments for output preferences
404
405 Returns:
406 None
407 """
408 # Print parameter info
409 log = OrderedDict()
410 log['START'] = 'CreateGermlines'
411 log['DB_FILE'] = os.path.basename(db_file)
412 log['GERM_TYPES'] = germ_types if isinstance(germ_types, str) else ','.join(germ_types)
413 log['CLONED'] = 'True'
414 log['V_FIELD'] = v_field
415 log['SEQ_FIELD'] = seq_field
416 printLog(log)
417
418 # Get repertoire and open Db reader
419 repo_dict = getRepo(repo)
420 reader = readDbFile(db_file, ig=False)
421
422 # Exit if V call field does not exist in reader
423 if v_field not in reader.fieldnames:
424 sys.exit('Error: V field does not exist in input database file.')
425
426 # Define log handle
427 if out_args['log_file'] is None:
428 log_handle = None
429 else:
430 log_handle = open(out_args['log_file'], 'w')
431
432 add_fields = []
433 seq_type = seq_field.split('_')[-1]
434 if 'full' in germ_types: add_fields += ['GERMLINE_' + seq_type]
435 if 'dmask' in germ_types: add_fields += ['GERMLINE_' + seq_type + '_D_MASK']
436 if 'vonly' in germ_types: add_fields += ['GERMLINE_' + seq_type + '_V_REGION']
437
438 # Create output file handle and Db writer
439 writers = {}
440 pass_handle = getOutputHandle(db_file, 'germ-pass', out_dir=out_args['out_dir'],
441 out_name=out_args['out_name'], out_type=out_args['out_type'])
442 writers['pass'] = getDbWriter(pass_handle, db_file, add_fields=add_fields)
443
444 if out_args['failed']:
445 fail_handle = getOutputHandle(db_file, 'germ-fail', out_dir=out_args['out_dir'],
446 out_name=out_args['out_name'], out_type=out_args['out_type'])
447 writers['fail'] = getDbWriter(fail_handle, db_file, add_fields=add_fields)
448 else:
449 fail_handle = None
450 writers['fail'] = None
451
452 # Initialize time and total count for progress bar
453 start_time = time()
454 rec_count = countDbFile(db_file)
455 counts = {}
456 clone_count = counts['pass'] = counts['fail'] = 0
457 # Iterate over rows
458 clone = 'initial'
459 clone_dict = OrderedDict()
460 for i,row in enumerate(reader):
461 # Print progress
462 printProgress(i, rec_count, 0.05, start_time)
463
464 # Clone isn't over yet
465 if row.get('CLONE','') == clone:
466 clone_dict[row["SEQUENCE_ID"]] = row
467 # Clone just finished
468 elif clone_dict:
469 clone_count += 1
470 result_log = makeCloneGermline(clone, clone_dict, repo_dict, germ_types,
471 v_field, seq_field, counts, writers, out_args)
472 printLog(result_log, handle=log_handle)
473 # Now deal with current row (first of next clone)
474 clone = row['CLONE']
475 clone_dict = OrderedDict([(row['SEQUENCE_ID'],row)])
476 # Last case is only for first row of file
477 else:
478 clone = row['CLONE']
479 clone_dict = OrderedDict([(row['SEQUENCE_ID'],row)])
480 clone_count += 1
481 result_log = makeCloneGermline(clone, clone_dict, repo_dict, germ_types, v_field,
482 seq_field, counts, writers, out_args)
483 printLog(result_log, handle=log_handle)
484
485 # Print log
486 printProgress(i+1, rec_count, 0.05, start_time)
487 log = OrderedDict()
488 log['OUTPUT'] = os.path.basename(pass_handle.name)
489 log['CLONES'] = clone_count
490 log['RECORDS'] = rec_count
491 log['PASS'] = counts['pass']
492 log['FAIL'] = counts['fail']
493 log['END'] = 'CreateGermlines'
494 printLog(log)
495
496 # Close file handles
497 pass_handle.close()
498 if fail_handle is not None: fail_handle.close()
499 if log_handle is not None: log_handle.close()
500
501
502 def getArgParser():
503 """
504 Defines the ArgumentParser
505
506 Arguments:
507 None
508
509 Returns:
510 an ArgumentParser object
511 """
512 # Define input and output field help message
513 fields = dedent(
514 '''
515 output files:
516 germ-pass
517 database with assigned germline sequences.
518 germ-fail
519 database with records failing germline assignment.
520
521 required fields:
522 SEQUENCE_ID, SEQUENCE_INPUT, SEQUENCE_VDJ or SEQUENCE_IMGT,
523 V_CALL or V_CALL_GENOTYPED, D_CALL, J_CALL,
524 V_SEQ_START, V_SEQ_LENGTH, V_GERM_START_IMGT, V_GERM_LENGTH_IMGT,
525 D_SEQ_START, D_SEQ_LENGTH, D_GERM_START, D_GERM_LENGTH,
526 J_SEQ_START, J_SEQ_LENGTH, J_GERM_START, J_GERM_LENGTH
527
528 optional fields:
529 CLONE
530
531 output fields:
532 GERMLINE_VDJ, GERMLINE_VDJ_D_MASK, GERMLINE_VDJ_V_REGION,
533 GERMLINE_IMGT, GERMLINE_IMGT_D_MASK, GERMLINE_IMGT_V_REGION
534 ''')
535
536 # Parent parser
537 parser_parent = getCommonArgParser(seq_in=False, seq_out=False, db_in=True,
538 annotation=False)
539 # Define argument parser
540 parser = ArgumentParser(description=__doc__, epilog=fields,
541 parents=[parser_parent],
542 formatter_class=CommonHelpFormatter)
543 parser.add_argument('--version', action='version',
544 version='%(prog)s:' + ' %s-%s' %(__version__, __date__))
545
546 parser.add_argument('-r', nargs='+', action='store', dest='repo', required=True,
547 help='List of folders and/or fasta files with germline sequences.')
548 parser.add_argument('-g', action='store', dest='germ_types', default=default_germ_types,
549 nargs='+', choices=('full', 'dmask', 'vonly'),
550 help='Specify type(s) of germlines to include full germline, \
551 germline with D-region masked, or germline for V region only.')
552 parser.add_argument('--cloned', action='store_true', dest='cloned',
553 help='Specify to create only one germline per clone \
554 (assumes input file is sorted by clone column)')
555 parser.add_argument('--vf', action='store', dest='v_field', default=default_v_field,
556 help='Specify field to use for germline V call')
557 parser.add_argument('--sf', action='store', dest='seq_field', default=default_seq_field,
558 help='Specify field to use for sequence')
559
560 return parser
561
562
563 if __name__ == "__main__":
564 """
565 Parses command line arguments and calls main
566 """
567
568 # Parse command line arguments
569 parser = getArgParser()
570 args = parser.parse_args()
571 args_dict = parseCommonArgs(args)
572 del args_dict['db_files']
573 del args_dict['cloned']
574 args_dict['v_field'] = args_dict['v_field'].upper()
575 args_dict['seq_field'] = args_dict['seq_field'].upper()
576
577 for f in args.__dict__['db_files']:
578 args_dict['db_file'] = f
579 if args.__dict__['cloned']:
580 assembleCloneGermline(**args_dict)
581 else:
582 assembleEachGermline(**args_dict)