Mercurial > repos > davidvanzessen > shm_csr
changeset 52:22dddabe3637 draft
Uploaded
| author | davidvanzessen |
|---|---|
| date | Tue, 23 May 2017 08:32:58 -0400 |
| parents | 8fa8836bd605 |
| children | 3be28ac82909 |
| files | change_o/DefineClones.py change_o/MakeDb.py |
| diffstat | 2 files changed, 547 insertions(+), 946 deletions(-) [+] |
line wrap: on
line diff
--- a/change_o/DefineClones.py Tue May 16 08:52:48 2017 -0400 +++ b/change_o/DefineClones.py Tue May 23 08:32:58 2017 -0400 @@ -10,6 +10,7 @@ import os import re import sys +import csv import numpy as np from argparse import ArgumentParser from collections import OrderedDict @@ -25,56 +26,108 @@ from presto.Multiprocessing import manageProcesses from presto.Sequence import getDNAScoreDict from changeo.Commandline import CommonHelpFormatter, getCommonArgParser, parseCommonArgs -from changeo.Distance import getDNADistMatrix, getAADistMatrix, \ - hs1f_model, m1n_model, hs5f_model, \ - calcDistances, formClusters +from changeo.Distance import distance_models, calcDistances, formClusters from changeo.IO import getDbWriter, readDbFile, countDbFile from changeo.Multiprocessing import DbData, DbResult +## Set maximum field size for csv.reader +csv.field_size_limit(sys.maxsize) + # Defaults default_translate = False default_distance = 0.0 -default_bygroup_model = 'hs1f' +default_index_mode = 'gene' +default_index_action = 'set' +default_bygroup_model = 'ham' default_hclust_model = 'chen2010' default_seq_field = 'JUNCTION' default_norm = 'len' default_sym = 'avg' default_linkage = 'single' - -# TODO: should be in Distance, but need to be after function definitions -# Amino acid Hamming distance -aa_model = getAADistMatrix(mask_dist=1, gap_dist=0) - -# DNA Hamming distance -ham_model = getDNADistMatrix(mask_dist=0, gap_dist=0) +choices_bygroup_model = ('ham', 'aa', 'hh_s1f', 'hh_s5f', 'mk_rs1nf', 'mk_rs5nf', 'hs1f_compat', 'm1n_compat') -# TODO: this function is an abstraction to facilitate later cleanup -def getModelMatrix(model): +def indexByIdentity(index, key, rec, fields=None): """ - Simple wrapper to get distance matrix from model name + Updates a preclone index with a simple key + + Arguments: + index = preclone index from indexJunctions + key = index key + rec = IgRecord to add to the index + fields = additional annotation fields to use to group preclones; + if None use only V, J and junction length + + Returns: + None. Updates index with new key and records. + """ + index.setdefault(tuple(key), []).append(rec) + + +def indexByUnion(index, key, rec, fields=None): + """ + Updates a preclone index with the union of nested keys Arguments: - model = model name + index = preclone index from indexJunctions + key = index key + rec = IgRecord to add to the index + fields = additional annotation fields to use to group preclones; + if None use only V, J and junction length - Return: - a pandas.DataFrame containing the character distance matrix + Returns: + None. Updates index with new key and records. """ - if model == 'aa': - return(aa_model) - elif model == 'ham': - return(ham_model) - elif model == 'm1n': - return(m1n_model) - elif model == 'hs1f': - return(hs1f_model) - elif model == 'hs5f': - return(hs5f_model) + # List of values for this/new key + val = [rec] + f_range = list(range(2, 3 + (len(fields) if fields else 0))) + + # See if field/junction length combination exists in index + outer_dict = index + for field in f_range: + try: + outer_dict = outer_dict[key[field]] + except (KeyError): + outer_dict = None + break + # If field combination exists, look through Js + j_matches = [] + if outer_dict is not None: + for j in outer_dict.keys(): + if not set(key[1]).isdisjoint(set(j)): + key[1] = tuple(set(key[1]).union(set(j))) + j_matches += [j] + # If J overlap exists, look through Vs for each J + for j in j_matches: + v_matches = [] + # Collect V matches for this J + for v in outer_dict[j].keys(): + if not set(key[0]).isdisjoint(set(v)): + key[0] = tuple(set(key[0]).union(set(v))) + v_matches += [v] + # If there are V overlaps for this J, pop them out + if v_matches: + val += list(chain(*(outer_dict[j].pop(v) for v in v_matches))) + # If the J dict is now empty, remove it + if not outer_dict[j]: + outer_dict.pop(j, None) + + # Add value(s) into index nested dictionary + # OMG Python pointers are the best! + # Add field dictionaries into index + outer_dict = index + for field in f_range: + outer_dict.setdefault(key[field], {}) + outer_dict = outer_dict[key[field]] + # Add J, then V into index + if key[1] in outer_dict: + outer_dict[key[1]].update({key[0]: val}) else: - sys.stderr.write('Unrecognized distance model: %s.\n' % model) + outer_dict[key[1]] = {key[0]: val} -def indexJunctions(db_iter, fields=None, mode='gene', action='first'): +def indexJunctions(db_iter, fields=None, mode=default_index_mode, + action=default_index_action): """ Identifies preclonal groups by V, J and junction length @@ -90,27 +143,47 @@ Returns: a dictionary of {(V, J, junction length):[IgRecords]} """ + # print(fields) # Define functions for grouping keys if mode == 'allele' and fields is None: def _get_key(rec, act): - return (rec.getVAllele(act), rec.getJAllele(act), - None if rec.junction is None else len(rec.junction)) + return [rec.getVAllele(act), rec.getJAllele(act), + None if rec.junction is None else len(rec.junction)] elif mode == 'gene' and fields is None: def _get_key(rec, act): - return (rec.getVGene(act), rec.getJGene(act), - None if rec.junction is None else len(rec.junction)) + return [rec.getVGene(act), rec.getJGene(act), + None if rec.junction is None else len(rec.junction)] elif mode == 'allele' and fields is not None: def _get_key(rec, act): vdj = [rec.getVAllele(act), rec.getJAllele(act), None if rec.junction is None else len(rec.junction)] ann = [rec.toDict().get(k, None) for k in fields] - return tuple(chain(vdj, ann)) + return list(chain(vdj, ann)) elif mode == 'gene' and fields is not None: def _get_key(rec, act): vdj = [rec.getVGene(act), rec.getJGene(act), None if rec.junction is None else len(rec.junction)] ann = [rec.toDict().get(k, None) for k in fields] - return tuple(chain(vdj, ann)) + return list(chain(vdj, ann)) + + # Function to flatten nested dictionary + def _flatten_dict(d, parent_key=''): + items = [] + for k, v in d.items(): + new_key = parent_key + [k] if parent_key else [k] + if isinstance(v, dict): + items.extend(_flatten_dict(v, new_key).items()) + else: + items.append((new_key, v)) + flat_dict = {None if None in i[0] else tuple(i[0]): i[1] for i in items} + return flat_dict + + if action == 'first': + index_func = indexByIdentity + elif action == 'set': + index_func = indexByUnion + else: + sys.stderr.write('Unrecognized action: %s.\n' % action) start_time = time() clone_index = {} @@ -127,35 +200,16 @@ # Assigned passed preclone records to key and failed to index None if all([k is not None and k != '' for k in key]): - #print key - # TODO: Has much slow. Should have less slow. - if action == 'set': - - f_range = list(range(2, 3 + (len(fields) if fields else 0))) - vdj_range = list(range(2)) - - # Check for any keys that have matching columns and junction length and overlapping genes/alleles - to_remove = [] - if len(clone_index) > (1 if None in clone_index else 0) and key not in clone_index: - key = list(key) - for k in clone_index: - if k is not None and all([key[i] == k[i] for i in f_range]): - if all([not set(key[i]).isdisjoint(set(k[i])) for i in vdj_range]): - for i in vdj_range: key[i] = tuple(set(key[i]).union(set(k[i]))) - to_remove.append(k) - - # Remove original keys, replace with union of all genes/alleles and append values to new key - val = [rec] - val += list(chain(*(clone_index.pop(k) for k in to_remove))) - clone_index[tuple(key)] = clone_index.get(tuple(key),[]) + val - - elif action == 'first': - clone_index.setdefault(key, []).append(rec) + # Update index dictionary + index_func(clone_index, key, rec, fields) else: clone_index.setdefault(None, []).append(rec) printProgress(rec_count, step=1000, start_time=start_time, end=True) + if action == 'set': + clone_index = _flatten_dict(clone_index) + return clone_index @@ -179,15 +233,20 @@ a dictionary of lists defining {clone number: [IgRecords clonal group]} """ # Get distance matrix if not provided - if dist_mat is None: dist_mat = getModelMatrix(model) + if dist_mat is None: + try: + dist_mat = distance_models[model] + except KeyError: + sys.exit('Unrecognized distance model: %s' % args_dict['model']) + # TODO: can be cleaned up with abstract model class # Determine length of n-mers - if model in ['hs1f', 'm1n', 'aa', 'ham']: + if model in ['hs1f_compat', 'm1n_compat', 'aa', 'ham', 'hh_s1f', 'mk_rs1nf']: nmer_len = 1 - elif model in ['hs5f']: + elif model in ['hh_s5f', 'mk_rs5nf']: nmer_len = 5 else: - sys.stderr.write('Unrecognized distance model: %s.\n' % model) + sys.exit('Unrecognized distance model: %s.\n' % model) # Define unique junction mapping seq_map = {} @@ -197,7 +256,7 @@ if len(seq) == 0: return None - seq = re.sub('[\.-]','N', str(seq)) + seq = re.sub('[\.-]', 'N', str(seq)) if model == 'aa': seq = translate(seq) seq_map.setdefault(seq, []).append(ig) @@ -210,7 +269,7 @@ seqs = list(seq_map.keys()) # Calculate pairwise distance matrix - dists = calcDistances(seqs, nmer_len, dist_mat, norm, sym) + dists = calcDistances(seqs, nmer_len, dist_mat, sym=sym, norm=norm) # Perform hierarchical clustering clusters = formClusters(dists, linkage, distance) @@ -449,6 +508,7 @@ # Define result object for iteration and get data records records = data.data + # print(data.id) result = DbResult(data.id, records) # Check for invalid data (due to failed indexing) and add failed result @@ -901,7 +961,7 @@ database with records failing clonal grouping. required fields: - SEQUENCE_ID, V_CALL or V_CALL_GENOTYPED, D_CALL, J_CALL, JUNCTION_LENGTH + SEQUENCE_ID, V_CALL or V_CALL_GENOTYPED, D_CALL, J_CALL, JUNCTION <field> sequence field specified by the --sf parameter @@ -926,32 +986,36 @@ # Distance cloning method parser_bygroup = subparsers.add_parser('bygroup', parents=[parser_parent], - formatter_class=CommonHelpFormatter, - help='''Defines clones as having same V assignment, - J assignment, and junction length with - specified substitution distance model.''') + formatter_class=CommonHelpFormatter, + help='''Defines clones as having same V assignment, + J assignment, and junction length with + specified substitution distance model.''', + description='''Defines clones as having same V assignment, + J assignment, and junction length with + specified substitution distance model.''') parser_bygroup.add_argument('-f', nargs='+', action='store', dest='fields', default=None, help='Additional fields to use for grouping clones (non VDJ)') parser_bygroup.add_argument('--mode', action='store', dest='mode', - choices=('allele', 'gene'), default='gene', + choices=('allele', 'gene'), default=default_index_mode, help='''Specifies whether to use the V(D)J allele or gene for initial grouping.''') - parser_bygroup.add_argument('--act', action='store', dest='action', default='set', - choices=('first', 'set'), + parser_bygroup.add_argument('--act', action='store', dest='action', + choices=('first', 'set'), default=default_index_action, help='''Specifies how to handle multiple V(D)J assignments for initial grouping.''') parser_bygroup.add_argument('--model', action='store', dest='model', - choices=('aa', 'ham', 'm1n', 'hs1f', 'hs5f'), + choices=choices_bygroup_model, default=default_bygroup_model, - help='''Specifies which substitution model to use for - calculating distance between sequences. Where m1n is the - mouse single nucleotide transition/trasversion model - of Smith et al, 1996; hs1f is the human single - nucleotide model derived from Yaari et al, 2013; hs5f - is the human S5F model of Yaari et al, 2013; ham is - nucleotide Hamming distance; and aa is amino acid - Hamming distance. The hs5f data should be - considered experimental.''') + help='''Specifies which substitution model to use for calculating distance + between sequences. The "ham" model is nucleotide Hamming distance and + "aa" is amino acid Hamming distance. The "hh_s1f" and "hh_s5f" models are + human specific single nucleotide and 5-mer content models, respectively, + from Yaari et al, 2013. The "mk_rs1nf" and "mk_rs5nf" models are + mouse specific single nucleotide and 5-mer content models, respectively, + from Cui et al, 2016. The "m1n_compat" and "hs1f_compat" models are + deprecated models provided backwards compatibility with the "m1n" and + "hs1f" models in Change-O v0.3.3 and SHazaM v0.1.4. Both + 5-mer models should be considered experimental.''') parser_bygroup.add_argument('--dist', action='store', dest='distance', type=float, default=default_distance, help='The distance threshold for clonal grouping') @@ -977,22 +1041,25 @@ parser_bygroup.set_defaults(group_func=indexJunctions) parser_bygroup.set_defaults(clone_func=distanceClones) - - # Hierarchical clustering cloning method - parser_hclust = subparsers.add_parser('hclust', parents=[parser_parent], + # Chen2010 + parser_chen = subparsers.add_parser('chen2010', parents=[parser_parent], formatter_class=CommonHelpFormatter, - help='Defines clones by specified distance metric on CDR3s and \ - cutting of hierarchical clustering tree') -# parser_hclust.add_argument('-f', nargs='+', action='store', dest='fields', default=None, -# help='Fields to use for grouping clones (non VDJ)') - parser_hclust.add_argument('--method', action='store', dest='method', - choices=('chen2010', 'ademokun2011'), default=default_hclust_model, - help='Specifies which cloning method to use for calculating distance \ - between CDR3s, computing linkage, and cutting clusters') - parser_hclust.set_defaults(feed_func=feedQueueClust) - parser_hclust.set_defaults(work_func=processQueueClust) - parser_hclust.set_defaults(collect_func=collectQueueClust) - parser_hclust.set_defaults(cluster_func=hierClust) + help='''Defines clones by method specified in Chen, 2010.''', + description='''Defines clones by method specified in Chen, 2010.''') + parser_chen.set_defaults(feed_func=feedQueueClust) + parser_chen.set_defaults(work_func=processQueueClust) + parser_chen.set_defaults(collect_func=collectQueueClust) + parser_chen.set_defaults(cluster_func=hierClust) + + # Ademokun2011 + parser_ade = subparsers.add_parser('ademokun2011', parents=[parser_parent], + formatter_class=CommonHelpFormatter, + help='''Defines clones by method specified in Ademokun, 2011.''', + description='''Defines clones by method specified in Ademokun, 2011.''') + parser_ade.set_defaults(feed_func=feedQueueClust) + parser_ade.set_defaults(work_func=processQueueClust) + parser_ade.set_defaults(collect_func=collectQueueClust) + parser_ade.set_defaults(cluster_func=hierClust) return parser @@ -1023,8 +1090,11 @@ 'linkage': args_dict['linkage'], 'seq_field': args_dict['seq_field']} - # TODO: can be cleaned up with abstract model class - args_dict['clone_args']['dist_mat'] = getModelMatrix(args_dict['model']) + # Get distance matrix + try: + args_dict['clone_args']['dist_mat'] = distance_models[args_dict['model']] + except KeyError: + sys.exit('Unrecognized distance model: %s' % args_dict['model']) del args_dict['fields'] del args_dict['action'] @@ -1037,16 +1107,17 @@ del args_dict['seq_field'] # Define clone_args - if args.command == 'hclust': - dist_funcs = {'chen2010':distChen2010, 'ademokun2011':distAdemokun2011} - args_dict['clone_func'] = dist_funcs[args_dict['method']] - args_dict['cluster_args'] = {'method': args_dict['method']} - #del args_dict['fields'] - del args_dict['method'] + if args.command == 'chen2010': + args_dict['clone_func'] = distChen2010 + args_dict['cluster_args'] = {'method': args.command } + + if args.command == 'ademokun2011': + args_dict['clone_func'] = distAdemokun2011 + args_dict['cluster_args'] = {'method': args.command } # Call defineClones del args_dict['command'] del args_dict['db_files'] for f in args.__dict__['db_files']: args_dict['db_file'] = f - defineClones(**args_dict) \ No newline at end of file + defineClones(**args_dict)
--- a/change_o/MakeDb.py Tue May 16 08:52:48 2017 -0400 +++ b/change_o/MakeDb.py Tue May 23 08:32:58 2017 -0400 @@ -7,765 +7,164 @@ from changeo import __version__, __date__ # Imports -import csv import os -import re import sys -import pandas as pd -import tarfile -import zipfile from argparse import ArgumentParser from collections import OrderedDict -from itertools import groupby -from shutil import rmtree -from tempfile import mkdtemp from textwrap import dedent from time import time from Bio import SeqIO -from Bio.Seq import Seq -from Bio.Alphabet import IUPAC # Presto and changeo imports from presto.Defaults import default_out_args from presto.Annotation import parseAnnotation -from presto.IO import countSeqFile, printLog, printProgress +from presto.IO import countSeqFile, printLog, printMessage, printProgress, readSeqFile from changeo.Commandline import CommonHelpFormatter, getCommonArgParser, parseCommonArgs -from changeo.IO import getDbWriter, countDbFile, getRepo -from changeo.Receptor import IgRecord, parseAllele, v_allele_regex, d_allele_regex, \ - j_allele_regex - -# Default parameters -default_delimiter = ('\t', ',', '-') - - -def gapV(ig_dict, repo_dict): - """ - Insert gaps into V region and update alignment information - - Arguments: - ig_dict : Dictionary of parsed IgBlast output - repo_dict : Dictionary of IMGT gapped germline sequences - - Returns: - dict : Updated with SEQUENCE_IMGT, V_GERM_START_IMGT, and V_GERM_LENGTH_IMGT fields - """ - - seq_imgt = '.' * (int(ig_dict['V_GERM_START_VDJ'])-1) + ig_dict['SEQUENCE_VDJ'] - - # Find gapped germline V segment - vgene = parseAllele(ig_dict['V_CALL'], v_allele_regex, 'first') - vkey = (vgene, ) - #TODO: Figure out else case - if vkey in repo_dict: - vgap = repo_dict[vkey] - # Iterate over gaps in the germline segment - gaps = re.finditer(r'\.', vgap) - gapcount = int(ig_dict['V_GERM_START_VDJ'])-1 - for gap in gaps: - i = gap.start() - # Break if gap begins after V region - if i >= ig_dict['V_GERM_LENGTH_VDJ'] + gapcount: - break - # Insert gap into IMGT sequence - seq_imgt = seq_imgt[:i] + '.' + seq_imgt[i:] - # Update gap counter - gapcount += 1 - ig_dict['SEQUENCE_IMGT'] = seq_imgt - # Update IMGT positioning information for V - ig_dict['V_GERM_START_IMGT'] = 1 - ig_dict['V_GERM_LENGTH_IMGT'] = ig_dict['V_GERM_LENGTH_VDJ'] + gapcount - - return ig_dict +from changeo.IO import countDbFile, extractIMGT, getDbWriter, readRepo +from changeo.Parsers import IgBLASTReader, IMGTReader, IHMMuneReader, getIDforIMGT -def getIMGTJunc(ig_dict, repo_dict): - """ - Identify junction region by IMGT definition - - Arguments: - ig_dict : Dictionary of parsed IgBlast output - repo_dict : Dictionary of IMGT gapped germline sequences - - Returns: - dict : Updated with JUNCTION_LENGTH_IMGT and JUNCTION_IMGT fields - """ - # Find germline J segment - jgene = parseAllele(ig_dict['J_CALL'], j_allele_regex, 'first') - jkey = (jgene, ) - #TODO: Figure out else case - if jkey in repo_dict: - # Get germline J sequence - jgerm = repo_dict[jkey] - jgerm = jgerm[:ig_dict['J_GERM_START']+ig_dict['J_GERM_LENGTH']-1] - # Look for (F|W)GXG aa motif in nt sequence - motif = re.search(r'T(TT|TC|GG)GG[ACGT]{4}GG[AGCT]', jgerm) - aa_end = len(ig_dict['SEQUENCE_IMGT']) - #TODO: Figure out else case - if motif: - # print('\n', motif.group()) - aa_end = motif.start() - len(jgerm) + 3 - # Add fields to dict - ig_dict['JUNCTION'] = ig_dict['SEQUENCE_IMGT'][309:aa_end] - ig_dict['JUNCTION_LENGTH'] = len(ig_dict['JUNCTION']) - - return ig_dict - - -def getRegions(ig_dict): +def getFilePrefix(aligner_output, out_args): """ - Identify FWR and CDR regions by IMGT definition - - Arguments: - ig_dict : Dictionary of parsed alignment output - - Returns: - dict : Updated with FWR1_IMGT, FWR2_IMGT, FWR3_IMGT, FWR4_IMGT, - CDR1_IMGT, CDR2_IMGT, and CDR3_IMGT fields - """ - try: - seq_len = len(ig_dict['SEQUENCE_IMGT']) - ig_dict['FWR1_IMGT'] = ig_dict['SEQUENCE_IMGT'][0:min(78,seq_len)] - except (KeyError, IndexError): - return ig_dict - - try: ig_dict['CDR1_IMGT'] = ig_dict['SEQUENCE_IMGT'][78:min(114, seq_len)] - except (IndexError): return ig_dict - - try: ig_dict['FWR2_IMGT'] = ig_dict['SEQUENCE_IMGT'][114:min(165, seq_len)] - except (IndexError): return ig_dict - - try: ig_dict['CDR2_IMGT'] = ig_dict['SEQUENCE_IMGT'][165:min(195, seq_len)] - except (IndexError): return ig_dict - - try: ig_dict['FWR3_IMGT'] = ig_dict['SEQUENCE_IMGT'][195:min(312, seq_len)] - except (IndexError): return ig_dict - - try: - cdr3_end = 306 + ig_dict['JUNCTION_LENGTH'] - ig_dict['CDR3_IMGT'] = ig_dict['SEQUENCE_IMGT'][312:cdr3_end] - ig_dict['FWR4_IMGT'] = ig_dict['SEQUENCE_IMGT'][cdr3_end:] - except (KeyError, IndexError): - return ig_dict - - return ig_dict - - -def getSeqforIgBlast(seq_file): - """ - Fetch input sequences for IgBlast queries + Get file name prefix and create output directory Arguments: - seq_file = a fasta file of sequences input to IgBlast + aligner_output : aligner output file or directory. + out_args : dictionary of output arguments. Returns: - a dictionary of {ID:Seq} + str : file name prefix. """ - - seq_dict = SeqIO.index(seq_file, "fasta", IUPAC.ambiguous_dna) - - # Create a seq_dict ID translation using IDs truncate up to space or 50 chars - seqs = {} - for seq in seq_dict.values(): - seqs.update({seq.description:str(seq.seq)}) - - return seqs - + # Determine output directory + if not out_args['out_dir']: + out_dir = os.path.dirname(os.path.abspath(aligner_output)) + else: + out_dir = os.path.abspath(out_args['out_dir']) + if not os.path.exists(out_dir): + os.mkdir(out_dir) -def findLine(handle, query): - """ - Finds line with query string in file + # Determine file prefix + if out_args['out_name']: + file_prefix = out_args['out_name'] + else: + file_prefix = os.path.splitext(os.path.split(os.path.abspath(aligner_output))[1])[0] - Arguments: - handle = file handle in which to search for line - query = query string for which to search in file - - Returns: - line from handle in which query string was found - """ - for line in handle: - if(re.match(query, line)): - return line + return os.path.join(out_dir, file_prefix) -def extractIMGT(imgt_output): - """ - Extract necessary files from IMGT results, zipped or unzipped - - Arguments: - imgt_output = zipped file or unzipped folder output by IMGT - - Returns: - sorted list of filenames from which information will be read +def getSeqDict(seq_file): """ - #file_ext = os.path.splitext(imgt_output)[1].lower() - imgt_flags = ('1_Summary', '2_IMGT-gapped', '3_Nt-sequences', '6_Junction') - temp_dir = mkdtemp() - if zipfile.is_zipfile(imgt_output): - # Open zip file - imgt_zip = zipfile.ZipFile(imgt_output, 'r') - # Extract required files - imgt_files = sorted([n for n in imgt_zip.namelist() \ - if os.path.basename(n).startswith(imgt_flags)]) - imgt_zip.extractall(temp_dir, imgt_files) - # Define file list - imgt_files = [os.path.join(temp_dir, f) for f in imgt_files] - elif os.path.isdir(imgt_output): - # Find required files in folder - folder_files = [] - for root, dirs, files in os.walk(imgt_output): - folder_files.extend([os.path.join(os.path.abspath(root), f) for f in files]) - # Define file list - imgt_files = sorted([n for n in folder_files \ - if os.path.basename(n).startswith(imgt_flags)]) - elif tarfile.is_tarfile(imgt_output): - # Open zip file - imgt_tar = tarfile.open(imgt_output, 'r') - # Extract required files - imgt_files = sorted([n for n in imgt_tar.getnames() \ - if os.path.basename(n).startswith(imgt_flags)]) - imgt_tar.extractall(temp_dir, [imgt_tar.getmember(n) for n in imgt_files]) - # Define file list - imgt_files = [os.path.join(temp_dir, f) for f in imgt_files] - else: - sys.exit('ERROR: Unsupported IGMT output file. Must be either a zipped file (.zip), LZMA compressed tarfile (.txz) or a folder.') - - if len(imgt_files) > len(imgt_flags): # e.g. multiple 1_Summary files - sys.exit('ERROR: Wrong files in IMGT output %s.' % imgt_output) - elif len(imgt_files) < len(imgt_flags): - sys.exit('ERROR: Missing necessary file IMGT output %s.' % imgt_output) - - return temp_dir, imgt_files - - -# TODO: return a dictionary with keys determined by the comment strings in the blocks, thus avoiding problems with missing blocks -def readOneIgBlastResult(block): - """ - Parse a single IgBLAST query result + Create a dictionary from a sequence file. Arguments: - block = itertools groupby object of single result - - Returns: - None if no results, otherwise list of DataFrames for each result block - """ - results = list() - i = 0 - for match, subblock in groupby(block, lambda l: l=='\n'): - if not match: - # Strip whitespace and comments - sub = [s.strip() for s in subblock if not s.startswith('#')] - - # Continue on empty block - if not sub: continue - else: i += 1 - - # Split by tabs - sub = [s.split('\t') for s in sub] - - # Append list for "V-(D)-J rearrangement summary" (i == 1) - # And "V-(D)-J junction details" (i == 2) - # Otherwise append DataFrame of subblock - if i == 1 or i == 2: - results.append(sub[0]) - else: - df = pd.DataFrame(sub) - if not df.empty: results.append(df) - - return results if results else None - - -# TODO: needs more speeds. pandas is probably to blame. -def readIgBlast(igblast_output, seq_dict, repo_dict, - score_fields=False, region_fields=False): - """ - Reads IgBlast output - - Arguments: - igblast_output = IgBlast output file (format 7) - seq_dict = a dictionary of {ID:Seq} from input fasta file - repo_dict = dictionary of IMGT gapped germline sequences - score_fields = if True parse alignment scores - region_fields = if True add FWR and CDR region fields + seq_file : sequence file. Returns: - a generator of dictionaries containing alignment data + dict : sequence description as keys with Bio.SeqRecords as values. """ - - # Open IgBlast output file - with open(igblast_output) as f: - # Iterate over individual results (separated by # IGBLASTN) - for k1, block in groupby(f, lambda x: re.match('# IGBLASTN', x)): - block = list(block) - if not k1: - # TODO: move query name extraction into block parser readOneIgBlastResult(). - # Extract sequence ID - query_name = ' '.join(block[0].strip().split(' ')[2:]) - # Initialize db_gen to have ID and input sequence - db_gen = {'SEQUENCE_ID': query_name, - 'SEQUENCE_INPUT': seq_dict[query_name]} - - # Parse further sub-blocks - block_list = readOneIgBlastResult(block) - - # TODO: this is indented pretty far. should be a separate function. or several functions. - # If results exist, parse further to obtain full db_gen - if block_list is not None: - # Parse quality information - db_gen['STOP'] = 'T' if block_list[0][-4] == 'Yes' else 'F' - db_gen['IN_FRAME'] = 'T' if block_list[0][-3] == 'In-frame' else 'F' - db_gen['FUNCTIONAL'] = 'T' if block_list[0][-2] == 'Yes' else 'F' - if block_list[0][-1] == '-': - db_gen['SEQUENCE_INPUT'] = str(Seq(db_gen['SEQUENCE_INPUT'], - IUPAC.ambiguous_dna).reverse_complement()) - - # Parse V, D, and J calls - call_str = ' '.join(block_list[0]) - v_call = parseAllele(call_str, v_allele_regex, action='list') - d_call = parseAllele(call_str, d_allele_regex, action='list') - j_call = parseAllele(call_str, j_allele_regex, action='list') - db_gen['V_CALL'] = ','.join(v_call) if v_call is not None else 'None' - db_gen['D_CALL'] = ','.join(d_call) if d_call is not None else 'None' - db_gen['J_CALL'] = ','.join(j_call) if j_call is not None else 'None' - - # Parse junction sequence - # db_gen['JUNCTION_VDJ'] = re.sub('(N/A)|\[|\(|\)|\]', '', ''.join(block_list[1])) - # db_gen['JUNCTION_LENGTH_VDJ'] = len(db_gen['JUNCTION_VDJ']) - - # TODO: IgBLAST does a stupid and doesn't output block #3 sometimes. why? - # TODO: maybe we should fail these. they look craptastic. - #pd.set_option('display.width', 500) - #print query_name, len(block_list), hit_idx - #for i, x in enumerate(block_list): - # print '[%i]' % i - # print x - - # Parse segment start and stop positions - hit_df = block_list[-1] - - # Alignment info block - # 0: segment - # 1: query id - # 2: subject id - # 3: % identity - # 4: alignment length - # 5: mismatches - # 6: gap opens - # 7: gaps - # 8: q. start - # 9: q. end - # 10: s. start - # 11: s. end - # 12: evalue - # 13: bit score - # 14: query seq - # 15: subject seq - # 16: btop - - # If V call exists, parse V alignment information - seq_vdj = '' - if v_call is not None: - v_align = hit_df[hit_df[0] == 'V'].iloc[0] - # Germline positions - db_gen['V_GERM_START_VDJ'] = int(v_align[10]) - db_gen['V_GERM_LENGTH_VDJ'] = int(v_align[11]) - db_gen['V_GERM_START_VDJ'] + 1 - # Query sequence positions - db_gen['V_SEQ_START'] = int(v_align[8]) - db_gen['V_SEQ_LENGTH'] = int(v_align[9]) - db_gen['V_SEQ_START'] + 1 - - if int(v_align[6]) == 0: - db_gen['INDELS'] = 'F' - else: - db_gen['INDELS'] = 'T' - # Set functional to none so record gets tossed (junction will be wrong) - # db_gen['FUNCTIONAL'] = None - - # V alignment scores - if score_fields: - try: db_gen['V_SCORE'] = float(v_align[13]) - except (TypeError, ValueError): db_gen['V_SCORE'] = 'None' - - try: db_gen['V_IDENTITY'] = float(v_align[3]) / 100.0 - except (TypeError, ValueError): db_gen['V_IDENTITY'] = 'None' - - try: db_gen['V_EVALUE'] = float(v_align[12]) - except (TypeError, ValueError): db_gen['V_EVALUE'] = 'None' - - try: db_gen['V_BTOP'] = v_align[16] - except (TypeError, ValueError): db_gen['V_BTOP'] = 'None' + seq_dict = SeqIO.to_dict(readSeqFile(seq_file), + key_function=lambda x: x.description) - # Update VDJ sequence, removing insertions - start = 0 - for m in re.finditer(r'-', v_align[15]): - ins = m.start() - seq_vdj += v_align[14][start:ins] - start = ins + 1 - seq_vdj += v_align[14][start:] - - # TODO: needs to check that the V results are present before trying to determine N1_LENGTH from them. - # If D call exists, parse D alignment information - if d_call is not None: - d_align = hit_df[hit_df[0] == 'D'].iloc[0] - - # TODO: this is kinda gross. not sure how else to fix the alignment overlap problem though. - # Determine N-region length and amount of J overlap with V or D alignment - overlap = 0 - if v_call is not None: - n1_len = int(d_align[8]) - (db_gen['V_SEQ_START'] + db_gen['V_SEQ_LENGTH']) - if n1_len < 0: - db_gen['N1_LENGTH'] = 0 - overlap = abs(n1_len) - else: - db_gen['N1_LENGTH'] = n1_len - n1_start = (db_gen['V_SEQ_START'] + db_gen['V_SEQ_LENGTH']-1) - n1_end = int(d_align[8])-1 - seq_vdj += db_gen['SEQUENCE_INPUT'][n1_start:n1_end] - - # Query sequence positions - db_gen['D_SEQ_START'] = int(d_align[8]) + overlap - db_gen['D_SEQ_LENGTH'] = max(int(d_align[9]) - db_gen['D_SEQ_START'] + 1, 0) - - # Germline positions - db_gen['D_GERM_START'] = int(d_align[10]) + overlap - db_gen['D_GERM_LENGTH'] = max(int(d_align[11]) - db_gen['D_GERM_START'] + 1, 0) - - # Update VDJ sequence, removing insertions - start = overlap - for m in re.finditer(r'-', d_align[15]): - ins = m.start() - seq_vdj += d_align[14][start:ins] - start = ins + 1 - seq_vdj += d_align[14][start:] - - # TODO: needs to check that the V results are present before trying to determine N1_LENGTH from them. - # If J call exists, parse J alignment information - if j_call is not None: - j_align = hit_df[hit_df[0] == 'J'].iloc[0] - - # TODO: this is kinda gross. not sure how else to fix the alignment overlap problem though. - # Determine N-region length and amount of J overlap with V or D alignment - overlap = 0 - if d_call is not None: - n2_len = int(j_align[8]) - (db_gen['D_SEQ_START'] + db_gen['D_SEQ_LENGTH']) - if n2_len < 0: - db_gen['N2_LENGTH'] = 0 - overlap = abs(n2_len) - else: - db_gen['N2_LENGTH'] = n2_len - n2_start = (db_gen['D_SEQ_START']+db_gen['D_SEQ_LENGTH']-1) - n2_end = int(j_align[8])-1 - seq_vdj += db_gen['SEQUENCE_INPUT'][n2_start:n2_end] - elif v_call is not None: - n1_len = int(j_align[8]) - (db_gen['V_SEQ_START'] + db_gen['V_SEQ_LENGTH']) - if n1_len < 0: - db_gen['N1_LENGTH'] = 0 - overlap = abs(n1_len) - else: - db_gen['N1_LENGTH'] = n1_len - n1_start = (db_gen['V_SEQ_START']+db_gen['V_SEQ_LENGTH']-1) - n1_end = int(j_align[8])-1 - seq_vdj += db_gen['SEQUENCE_INPUT'][n1_start:n1_end] - else: - db_gen['N1_LENGTH'] = 0 - - # Query positions - db_gen['J_SEQ_START'] = int(j_align[8]) + overlap - db_gen['J_SEQ_LENGTH'] = max(int(j_align[9]) - db_gen['J_SEQ_START'] + 1, 0) - - # Germline positions - db_gen['J_GERM_START'] = int(j_align[10]) + overlap - db_gen['J_GERM_LENGTH'] = max(int(j_align[11]) - db_gen['J_GERM_START'] + 1, 0) - - # J alignment scores - if score_fields: - try: db_gen['J_SCORE'] = float(j_align[13]) - except (TypeError, ValueError): db_gen['J_SCORE'] = 'None' - - try: db_gen['J_IDENTITY'] = float(j_align[3]) / 100.0 - except (TypeError, ValueError): db_gen['J_IDENTITY'] = 'None' - - try: db_gen['J_EVALUE'] = float(j_align[12]) - except (TypeError, ValueError): db_gen['J_EVALUE'] = 'None' - - try: db_gen['J_BTOP'] = j_align[16] - except (TypeError, ValueError): db_gen['J_BTOP'] = 'None' - - # Update VDJ sequence, removing insertions - start = overlap - for m in re.finditer(r'-', j_align[15]): - ins = m.start() - seq_vdj += j_align[14][start:ins] - start = ins + 1 - seq_vdj += j_align[14][start:] - - db_gen['SEQUENCE_VDJ'] = seq_vdj - - # Create IMGT-gapped sequence and infer IMGT junction - if v_call is not None: - db_gen = gapV(db_gen, repo_dict) - if j_call is not None: - db_gen = getIMGTJunc(db_gen, repo_dict) - - # FWR and CDR regions - if region_fields: getRegions(db_gen) - - yield IgRecord(db_gen) + return seq_dict -# TODO: should be more readable -def readIMGT(imgt_files, score_fields=False, region_fields=False): - """ - Reads IMGT/HighV-Quest output - - Arguments: - imgt_files = IMGT/HighV-Quest output files 1, 2, 3, and 6 - score_fields = if True parse alignment scores - region_fields = if True add FWR and CDR region fields - - Returns: - a generator of dictionaries containing alignment data +def writeDb(db, fields, file_prefix, total_count, id_dict=None, no_parse=True, partial=False, + out_args=default_out_args): """ - imgt_iters = [csv.DictReader(open(f, 'rU'), delimiter='\t') for f in imgt_files] - # Create a dictionary for each sequence alignment and yield its generator - for sm, gp, nt, jn in zip(*imgt_iters): - if len(set([sm['Sequence ID'], - gp['Sequence ID'], - nt['Sequence ID'], - jn['Sequence ID']])) != 1: - sys.exit('Error: IMGT files are corrupt starting with Summary file record %s' \ - % sm['Sequence ID']) - - db_gen = {'SEQUENCE_ID': sm['Sequence ID'], - 'SEQUENCE_INPUT': sm['Sequence']} - - if 'No results' not in sm['Functionality']: - db_gen['FUNCTIONAL'] = ['?','T','F'][('productive' in sm['Functionality']) + - ('unprod' in sm['Functionality'])] - db_gen['IN_FRAME'] = ['?','T','F'][('in-frame' in sm['JUNCTION frame']) + - ('out-of-frame' in sm['JUNCTION frame'])], - db_gen['STOP'] = ['F','?','T'][('stop codon' in sm['Functionality comment']) + - ('unprod' in sm['Functionality'])] - db_gen['MUTATED_INVARIANT'] = ['F','?','T'][(any(('missing' in sm['Functionality comment'], - 'missing' in sm['V-REGION potential ins/del']))) + - ('unprod' in sm['Functionality'])] - db_gen['INDELS'] = ['F','T'][any((sm['V-REGION potential ins/del'], - sm['V-REGION insertions'], - sm['V-REGION deletions']))] - - db_gen['SEQUENCE_VDJ'] = nt['V-D-J-REGION'] if nt['V-D-J-REGION'] else nt['V-J-REGION'] - db_gen['SEQUENCE_IMGT'] = gp['V-D-J-REGION'] if gp['V-D-J-REGION'] else gp['V-J-REGION'] - - db_gen['V_CALL'] = re.sub('\sor\s', ',', re.sub(',', '', gp['V-GENE and allele'])) - db_gen['D_CALL'] = re.sub('\sor\s', ',', re.sub(',', '', gp['D-GENE and allele'])) - db_gen['J_CALL'] = re.sub('\sor\s', ',', re.sub(',', '', gp['J-GENE and allele'])) - - v_seq_length = len(nt['V-REGION']) if nt['V-REGION'] else 0 - db_gen['V_SEQ_START'] = nt['V-REGION start'] - db_gen['V_SEQ_LENGTH'] = v_seq_length - db_gen['V_GERM_START_IMGT'] = 1 - db_gen['V_GERM_LENGTH_IMGT'] = len(gp['V-REGION']) if gp['V-REGION'] else 0 - - db_gen['N1_LENGTH'] = sum(int(i) for i in [jn["P3'V-nt nb"], - jn['N-REGION-nt nb'], - jn['N1-REGION-nt nb'], - jn["P5'D-nt nb"]] if i) - db_gen['D_SEQ_START'] = sum(int(i) for i in [1, v_seq_length, - jn["P3'V-nt nb"], - jn['N-REGION-nt nb'], - jn['N1-REGION-nt nb'], - jn["P5'D-nt nb"]] if i) - db_gen['D_SEQ_LENGTH'] = int(jn["D-REGION-nt nb"] or 0) - db_gen['D_GERM_START'] = int(jn["5'D-REGION trimmed-nt nb"] or 0) + 1 - db_gen['D_GERM_LENGTH'] = int(jn["D-REGION-nt nb"] or 0) - db_gen['N2_LENGTH'] = sum(int(i) for i in [jn["P3'D-nt nb"], - jn['N2-REGION-nt nb'], - jn["P5'J-nt nb"]] if i) - - db_gen['J_SEQ_START_IMGT'] = sum(int(i) for i in [1, v_seq_length, - jn["P3'V-nt nb"], - jn['N-REGION-nt nb'], - jn['N1-REGION-nt nb'], - jn["P5'D-nt nb"], - jn["D-REGION-nt nb"], - jn["P3'D-nt nb"], - jn['N2-REGION-nt nb'], - jn["P5'J-nt nb"]] if i) - db_gen['J_SEQ_LENGTH'] = len(nt['J-REGION']) if nt['J-REGION'] else 0 - db_gen['J_GERM_START'] = int(jn["5'J-REGION trimmed-nt nb"] or 0) + 1 - db_gen['J_GERM_LENGTH'] = len(gp['J-REGION']) if gp['J-REGION'] else 0 - - db_gen['JUNCTION_LENGTH'] = len(jn['JUNCTION']) if jn['JUNCTION'] else 0 - db_gen['JUNCTION'] = jn['JUNCTION'] - - # Alignment scores - if score_fields: - try: db_gen['V_SCORE'] = float(sm['V-REGION score']) - except (TypeError, ValueError): db_gen['V_SCORE'] = 'None' - - try: db_gen['V_IDENTITY'] = float(sm['V-REGION identity %']) / 100.0 - except (TypeError, ValueError): db_gen['V_IDENTITY'] = 'None' - - try: db_gen['J_SCORE'] = float(sm['J-REGION score']) - except (TypeError, ValueError): db_gen['J_SCORE'] = 'None' - - try: db_gen['J_IDENTITY'] = float(sm['J-REGION identity %']) / 100.0 - except (TypeError, ValueError): db_gen['J_IDENTITY'] = 'None' - - # FWR and CDR regions - if region_fields: getRegions(db_gen) - else: - db_gen['V_CALL'] = 'None' - db_gen['D_CALL'] = 'None' - db_gen['J_CALL'] = 'None' - - yield IgRecord(db_gen) - - -def getIDforIMGT(seq_file): - """ - Create a sequence ID translation using IMGT truncation - - Arguments: - seq_file = a fasta file of sequences input to IMGT - - Returns: - a dictionary of {truncated ID: full seq description} - """ - - # Create a seq_dict ID translation using IDs truncate up to space or 50 chars - ids = {} - for i, rec in enumerate(SeqIO.parse(seq_file, 'fasta', IUPAC.ambiguous_dna)): - if len(rec.description) <= 50: - id_key = rec.description - else: - id_key = re.sub('\||\s|!|&|\*|<|>|\?','_',rec.description[:50]) - ids.update({id_key:rec.description}) - - return ids - - -def writeDb(db_gen, file_prefix, total_count, id_dict={}, no_parse=True, - score_fields=False, region_fields=False, out_args=default_out_args): - """ - Writes tab-delimited database file in output directory + Writes tab-delimited database file in output directory. Arguments: - db_gen = a generator of IgRecord objects containing alignment data - file_prefix = directory and prefix for CLIP tab-delim file - total_count = number of records (for progress bar) - id_dict = a dictionary of {IMGT ID: full seq description} - no_parse = if ID is to be parsed for pRESTO output with default delimiters - score_fields = if True add alignment score fields to output file - region_fields = if True add FWR and CDR region fields to output file - out_args = common output argument dictionary from parseCommonArgs + db : a iterator of IgRecord objects containing alignment data. + fields : a list of ordered field names to write. + file_prefix : directory and prefix for CLIP tab-delim file. + total_count : number of records (for progress bar). + id_dict : a dictionary of the truncated sequence ID mapped to the full sequence ID. + no_parse : if ID is to be parsed for pRESTO output with default delimiters. + partial : if True put incomplete alignments in the pass file. + out_args : common output argument dictionary from parseCommonArgs. Returns: - None + None """ - pass_file = "%s_db-pass.tab" % file_prefix - fail_file = "%s_db-fail.tab" % file_prefix - ordered_fields = ['SEQUENCE_ID', - 'SEQUENCE_INPUT', - 'FUNCTIONAL', - 'IN_FRAME', - 'STOP', - 'MUTATED_INVARIANT', - 'INDELS', - 'V_CALL', - 'D_CALL', - 'J_CALL', - 'SEQUENCE_VDJ', - 'SEQUENCE_IMGT', - 'V_SEQ_START', - 'V_SEQ_LENGTH', - 'V_GERM_START_VDJ', - 'V_GERM_LENGTH_VDJ', - 'V_GERM_START_IMGT', - 'V_GERM_LENGTH_IMGT', - 'N1_LENGTH', - 'D_SEQ_START', - 'D_SEQ_LENGTH', - 'D_GERM_START', - 'D_GERM_LENGTH', - 'N2_LENGTH', - 'J_SEQ_START', - 'J_SEQ_LENGTH', - 'J_GERM_START', - 'J_GERM_LENGTH', - 'JUNCTION_LENGTH', - 'JUNCTION'] + # Function to check for valid records strictly + def _pass_strict(rec): + valid = [rec.v_call and rec.v_call != 'None', + rec.j_call and rec.j_call != 'None', + rec.functional is not None, + rec.seq_vdj, + rec.junction] + return all(valid) - if score_fields: - ordered_fields.extend(['V_SCORE', - 'V_IDENTITY', - 'V_EVALUE', - 'V_BTOP', - 'J_SCORE', - 'J_IDENTITY', - 'J_EVALUE', - 'J_BTOP']) + # Function to check for valid records loosely + def _pass_gentle(rec): + valid = [rec.v_call and rec.v_call != 'None', + rec.d_call and rec.d_call != 'None', + rec.j_call and rec.j_call != 'None'] + return any(valid) - if region_fields: - ordered_fields.extend(['FWR1_IMGT', 'FWR2_IMGT', 'FWR3_IMGT', 'FWR4_IMGT', - 'CDR1_IMGT', 'CDR2_IMGT', 'CDR3_IMGT']) - + # Set pass criteria + _pass = _pass_gentle if partial else _pass_strict - # TODO: This is not the best approach. should pass in output fields. - # Initiate passed handle - pass_handle = None + # Define output file names + pass_file = '%s_db-pass.tab' % file_prefix + fail_file = '%s_db-fail.tab' % file_prefix - # Open failed file - if out_args['failed']: - fail_handle = open(fail_file, 'wt') - fail_writer = getDbWriter(fail_handle, add_fields=['SEQUENCE_ID', 'SEQUENCE_INPUT']) - else: - fail_handle = None - fail_writer = None - - # Initialize counters and file + # Initiate handles, writers and counters + pass_handle = None + fail_handle = None pass_writer = None + fail_writer = None start_time = time() rec_count = pass_count = fail_count = 0 - for record in db_gen: - #printProgress(i + (total_count/2 if id_dict else 0), total_count, 0.05, start_time) - printProgress(rec_count, total_count, 0.05, start_time) - rec_count += 1 - # Count pass or fail - if (record.v_call == 'None' and record.j_call == 'None') or \ - record.functional is None or \ - not record.seq_vdj or \ - not record.junction: - # print(record.v_call, record.j_call, record.functional, record.junction) - fail_count += 1 - if fail_writer is not None: fail_writer.writerow(record.toDict()) - continue - else: - pass_count += 1 - - # Build sample sequence description - if record.id in id_dict: + # Validate and write output + printProgress(0, total_count, 0.05, start_time) + for i, record in enumerate(db, start=1): + + # Replace sequence description with full string, if required + if id_dict is not None and record.id in id_dict: record.id = id_dict[record.id] # Parse sequence description into new columns if not no_parse: - record.annotations = parseAnnotation(record.id, delimiter=out_args['delimiter']) - record.id = record.annotations['ID'] - del record.annotations['ID'] + try: + record.annotations = parseAnnotation(record.id, delimiter=out_args['delimiter']) + record.id = record.annotations['ID'] + del record.annotations['ID'] + + # TODO: This is not the best approach. should pass in output fields. + # If first record, use parsed description to define extra columns + if i == 1: fields.extend(list(record.annotations.keys())) + except IndexError: + # Could not parse pRESTO-style annotations so fall back to no parse + no_parse = True + sys.stderr.write('\nWARNING: Sequence annotation format not recognized. Sequence headers will not be parsed.\n') - # TODO: This is not the best approach. should pass in output fields. - # If first sequence, use parsed description to create new columns and initialize writer - if pass_writer is None: - if not no_parse: ordered_fields.extend(list(record.annotations.keys())) - pass_handle = open(pass_file, 'wt') - pass_writer = getDbWriter(pass_handle, add_fields=ordered_fields) + # Count pass or fail and write to appropriate file + if _pass(record): + # Open pass file + if pass_writer is None: + pass_handle = open(pass_file, 'wt') + pass_writer = getDbWriter(pass_handle, add_fields=fields) - # Write row to tab-delim CLIP file - pass_writer.writerow(record.toDict()) - - # Print log - #printProgress(i+1 + (total_count/2 if id_dict else 0), total_count, 0.05, start_time) - printProgress(rec_count, total_count, 0.05, start_time) + # Write row to pass file + pass_count += 1 + pass_writer.writerow(record.toDict()) + else: + # Open failed file + if out_args['failed'] and fail_writer is None: + fail_handle = open(fail_file, 'wt') + fail_writer = getDbWriter(fail_handle, add_fields=fields) + # Write row to fail file if specified + fail_count += 1 + if fail_writer is not None: + fail_writer.writerow(record.toDict()) + + # Print progress + printProgress(i, total_count, 0.05, start_time) + + # Print consol log log = OrderedDict() log['OUTPUT'] = pass_file log['PASS'] = pass_count @@ -777,146 +176,215 @@ if fail_handle is not None: fail_handle.close() -# TODO: may be able to merge with parseIMGT -def parseIgBlast(igblast_output, seq_file, repo, no_parse=True, score_fields=False, - region_fields=False, out_args=default_out_args): +# TODO: may be able to merge with other mains +def parseIMGT(aligner_output, seq_file=None, no_parse=True, partial=False, + parse_scores=False, parse_regions=False, parse_junction=False, + out_args=default_out_args): """ - Main for IgBlast aligned sample sequences + Main for IMGT aligned sample sequences. Arguments: - igblast_output = IgBlast output file to process - seq_file = fasta file input to IgBlast (from which to get sequence) - repo = folder with germline repertoire files - no_parse = if ID is to be parsed for pRESTO output with default delimiters - score_fields = if True add alignment score fields to output file - region_fields = if True add FWR and CDR region fields to output file - out_args = common output argument dictionary from parseCommonArgs + aligner_output : zipped file or unzipped folder output by IMGT. + seq_file : FASTA file input to IMGT (from which to get seqID). + no_parse : if ID is to be parsed for pRESTO output with default delimiters. + partial : If True put incomplete alignments in the pass file. + parse_scores : if True add alignment score fields to output file. + parse_regions : if True add FWR and CDR region fields to output file. + out_args : common output argument dictionary from parseCommonArgs. Returns: - None + None + """ + # Print parameter info + log = OrderedDict() + log['START'] = 'MakeDb' + log['ALIGNER'] = 'IMGT' + log['ALIGNER_OUTPUT'] = aligner_output + log['SEQ_FILE'] = os.path.basename(seq_file) if seq_file else '' + log['NO_PARSE'] = no_parse + log['PARTIAL'] = partial + log['SCORES'] = parse_scores + log['REGIONS'] = parse_regions + log['JUNCTION'] = parse_junction + printLog(log) + + start_time = time() + printMessage('Loading sequence files', start_time=start_time, width=25) + # Extract IMGT files + temp_dir, imgt_files = extractIMGT(aligner_output) + # Count records in IMGT files + total_count = countDbFile(imgt_files['summary']) + # Get (parsed) IDs from fasta file submitted to IMGT + id_dict = getIDforIMGT(seq_file) if seq_file else {} + printMessage('Done', start_time=start_time, end=True, width=25) + + # Parse IMGT output and write db + with open(imgt_files['summary'], 'r') as summary_handle, \ + open(imgt_files['gapped'], 'r') as gapped_handle, \ + open(imgt_files['ntseq'], 'r') as ntseq_handle, \ + open(imgt_files['junction'], 'r') as junction_handle: + parse_iter = IMGTReader(summary_handle, gapped_handle, ntseq_handle, junction_handle, + parse_scores=parse_scores, parse_regions=parse_regions, + parse_junction=parse_junction) + file_prefix = getFilePrefix(aligner_output, out_args) + writeDb(parse_iter, parse_iter.fields, file_prefix, total_count, id_dict=id_dict, + no_parse=no_parse, partial=partial, out_args=out_args) + + # Cleanup temp directory + temp_dir.cleanup() + + return None + + +# TODO: may be able to merge with other mains +def parseIgBLAST(aligner_output, seq_file, repo, no_parse=True, partial=False, + parse_regions=False, parse_scores=False, parse_igblast_cdr3=False, + out_args=default_out_args): + """ + Main for IgBLAST aligned sample sequences. + + Arguments: + aligner_output : IgBLAST output file to process. + seq_file : fasta file input to IgBlast (from which to get sequence). + repo : folder with germline repertoire files. + no_parse : if ID is to be parsed for pRESTO output with default delimiters. + partial : If True put incomplete alignments in the pass file. + parse_regions : if True add FWR and CDR fields to output file. + parse_scores : if True add alignment score fields to output file. + parse_igblast_cdr3 : if True parse CDR3 sequences generated by IgBLAST + out_args : common output argument dictionary from parseCommonArgs. + + Returns: + None """ # Print parameter info log = OrderedDict() log['START'] = 'MakeDB' log['ALIGNER'] = 'IgBlast' - log['ALIGN_RESULTS'] = os.path.basename(igblast_output) + log['ALIGNER_OUTPUT'] = os.path.basename(aligner_output) log['SEQ_FILE'] = os.path.basename(seq_file) log['NO_PARSE'] = no_parse - log['SCORE_FIELDS'] = score_fields - log['REGION_FIELDS'] = region_fields + log['PARTIAL'] = partial + log['SCORES'] = parse_scores + log['REGIONS'] = parse_regions printLog(log) + start_time = time() + printMessage('Loading sequence files', start_time=start_time, width=25) + # Count records in sequence file + total_count = countSeqFile(seq_file) # Get input sequence dictionary - seq_dict = getSeqforIgBlast(seq_file) + seq_dict = getSeqDict(seq_file) + # Create germline repo dictionary + repo_dict = readRepo(repo) + printMessage('Done', start_time=start_time, end=True, width=25) - # Formalize out_dir and file-prefix - if not out_args['out_dir']: - out_dir = os.path.split(igblast_output)[0] - else: - out_dir = os.path.abspath(out_args['out_dir']) - if not os.path.exists(out_dir): os.mkdir(out_dir) - if out_args['out_name']: - file_prefix = out_args['out_name'] - else: - file_prefix = os.path.basename(os.path.splitext(igblast_output)[0]) - file_prefix = os.path.join(out_dir, file_prefix) + # Parse and write output + with open(aligner_output, 'r') as f: + parse_iter = IgBLASTReader(f, seq_dict, repo_dict, + parse_scores=parse_scores, parse_regions=parse_regions, + parse_igblast_cdr3=parse_igblast_cdr3) + file_prefix = getFilePrefix(aligner_output, out_args) + writeDb(parse_iter, parse_iter.fields, file_prefix, total_count, + no_parse=no_parse, partial=partial, out_args=out_args) - total_count = countSeqFile(seq_file) - - # Create - repo_dict = getRepo(repo) - igblast_dict = readIgBlast(igblast_output, seq_dict, repo_dict, - score_fields=score_fields, region_fields=region_fields) - writeDb(igblast_dict, file_prefix, total_count, no_parse=no_parse, - score_fields=score_fields, region_fields=region_fields, out_args=out_args) + return None -# TODO: may be able to merge with parseIgBlast -def parseIMGT(imgt_output, seq_file=None, no_parse=True, score_fields=False, - region_fields=False, out_args=default_out_args): +# TODO: may be able to merge with other mains +def parseIHMM(aligner_output, seq_file, repo, no_parse=True, partial=False, + parse_scores=False, parse_regions=False, out_args=default_out_args): """ - Main for IMGT aligned sample sequences + Main for iHMMuneAlign aligned sample sequences. Arguments: - imgt_output = zipped file or unzipped folder output by IMGT - seq_file = FASTA file input to IMGT (from which to get seqID) - no_parse = if ID is to be parsed for pRESTO output with default delimiters - score_fields = if True add alignment score fields to output file - region_fields = if True add FWR and CDR region fields to output file - out_args = common output argument dictionary from parseCommonArgs - - Returns: - None + aligner_output : iHMMune-Align output file to process. + seq_file : fasta file input to iHMMuneAlign (from which to get sequence). + repo : folder with germline repertoire files. + no_parse : if ID is to be parsed for pRESTO output with default delimiters. + partial : If True put incomplete alignments in the pass file. + parse_scores : if True parse alignment scores. + parse_regions : if True add FWR and CDR region fields. + out_args : common output argument dictionary from parseCommonArgs. + + Returns: + None """ # Print parameter info log = OrderedDict() - log['START'] = 'MakeDb' - log['ALIGNER'] = 'IMGT' - log['ALIGN_RESULTS'] = imgt_output - log['SEQ_FILE'] = os.path.basename(seq_file) if seq_file else '' + log['START'] = 'MakeDB' + log['ALIGNER'] = 'iHMMune-Align' + log['ALIGNER_OUTPUT'] = os.path.basename(aligner_output) + log['SEQ_FILE'] = os.path.basename(seq_file) log['NO_PARSE'] = no_parse - log['SCORE_FIELDS'] = score_fields - log['REGION_FIELDS'] = region_fields + log['PARTIAL'] = partial + log['SCORES'] = parse_scores + log['REGIONS'] = parse_regions printLog(log) - - # Get individual IMGT result files - temp_dir, imgt_files = extractIMGT(imgt_output) - - # Formalize out_dir and file-prefix - if not out_args['out_dir']: - out_dir = os.path.dirname(os.path.abspath(imgt_output)) - else: - out_dir = os.path.abspath(out_args['out_dir']) - if not os.path.exists(out_dir): os.mkdir(out_dir) - if out_args['out_name']: - file_prefix = out_args['out_name'] - else: - file_prefix = os.path.splitext(os.path.split(os.path.abspath(imgt_output))[1])[0] - file_prefix = os.path.join(out_dir, file_prefix) - total_count = countDbFile(imgt_files[0]) - - # Get (parsed) IDs from fasta file submitted to IMGT - id_dict = getIDforIMGT(seq_file) if seq_file else {} - - # Create - imgt_dict = readIMGT(imgt_files, score_fields=score_fields, - region_fields=region_fields) - writeDb(imgt_dict, file_prefix, total_count, id_dict=id_dict, no_parse=no_parse, - score_fields=score_fields, region_fields=region_fields, out_args=out_args) + start_time = time() + printMessage('Loading sequence files', start_time=start_time, width=25) + # Count records in sequence file + total_count = countSeqFile(seq_file) + # Get input sequence dictionary + seq_dict = getSeqDict(seq_file) + # Create germline repo dictionary + repo_dict = readRepo(repo) + printMessage('Done', start_time=start_time, end=True, width=25) - # Delete temp directory - rmtree(temp_dir) + # Parse and write output + with open(aligner_output, 'r') as f: + parse_iter = IHMMuneReader(f, seq_dict, repo_dict, + parse_scores=parse_scores, parse_regions=parse_regions) + file_prefix = getFilePrefix(aligner_output, out_args) + writeDb(parse_iter, parse_iter.fields, file_prefix, total_count, + no_parse=no_parse, partial=partial, out_args=out_args) + + return None def getArgParser(): """ - Defines the ArgumentParser + Defines the ArgumentParser. - Arguments: - None - Returns: - an ArgumentParser object + argparse.ArgumentParser """ fields = dedent( ''' output files: db-pass - database of parsed alignment records. + database of alignment records with functionality information, + V and J calls, and a junction region. db-fail - database with records failing alignment. + database with records that fail due to no functionality information + (did not pass IMGT), no V call, no J call, or no junction region. - output fields: - SEQUENCE_ID, SEQUENCE_INPUT, FUNCTIONAL, IN_FRAME, STOP, MUTATED_INVARIANT, - INDELS, V_CALL, D_CALL, J_CALL, SEQUENCE_VDJ and/or SEQUENCE_IMGT, - V_SEQ_START, V_SEQ_LENGTH, V_GERM_START_VDJ and/or V_GERM_START_IMGT, - V_GERM_LENGTH_VDJ and/or V_GERM_LENGTH_IMGT, N1_LENGTH, - D_SEQ_START, D_SEQ_LENGTH, D_GERM_START, D_GERM_LENGTH, N2_LENGTH, + universal output fields: + SEQUENCE_ID, SEQUENCE_INPUT, SEQUENCE_VDJ, SEQUENCE_IMGT, + FUNCTIONAL, IN_FRAME, STOP, MUTATED_INVARIANT, INDELS, + V_CALL, D_CALL, J_CALL, + V_SEQ_START, V_SEQ_LENGTH, + D_SEQ_START, D_SEQ_LENGTH, D_GERM_START, D_GERM_LENGTH, J_SEQ_START, J_SEQ_LENGTH, J_GERM_START, J_GERM_LENGTH, - JUNCTION_LENGTH, JUNCTION, V_SCORE, V_IDENTITY, V_EVALUE, V_BTOP, - J_SCORE, J_IDENTITY, J_EVALUE, J_BTOP, FWR1_IMGT, FWR2_IMGT, FWR3_IMGT, - FWR4_IMGT, CDR1_IMGT, CDR2_IMGT, CDR3_IMGT + JUNCTION_LENGTH, JUNCTION, NP1_LENGTH, NP2_LENGTH, + FWR1_IMGT, FWR2_IMGT, FWR3_IMGT, FWR4_IMGT, + CDR1_IMGT, CDR2_IMGT, CDR3_IMGT + + imgt specific output fields: + V_GERM_START_IMGT, V_GERM_LENGTH_IMGT, + N1_LENGTH, N2_LENGTH, P3V_LENGTH, P5D_LENGTH, P3D_LENGTH, P5J_LENGTH, + D_FRAME, V_SCORE, V_IDENTITY, J_SCORE, J_IDENTITY, + + igblast specific output fields: + V_GERM_START_VDJ, V_GERM_LENGTH_VDJ, + V_EVALUE, V_SCORE, V_IDENTITY, V_BTOP, + J_EVALUE, J_SCORE, J_IDENTITY, J_BTOP. + CDR3_IGBLAST_NT, CDR3_IGBLAST_AA + + ihmm specific output fields: + V_GERM_START_VDJ, V_GERM_LENGTH_VDJ, + HMM_SCORE ''') # Define ArgumentParser @@ -933,11 +401,11 @@ parser_parent = getCommonArgParser(seq_in=False, seq_out=False, log=False) # IgBlast Aligner - parser_igblast = subparsers.add_parser('igblast', help='Process IgBlast output', - parents=[parser_parent], - formatter_class=CommonHelpFormatter) - parser_igblast.set_defaults(func=parseIgBlast) - parser_igblast.add_argument('-i', nargs='+', action='store', dest='aligner_files', + parser_igblast = subparsers.add_parser('igblast', parents=[parser_parent], + formatter_class=CommonHelpFormatter, + help='Process IgBLAST output.', + description='Process IgBLAST output.') + parser_igblast.add_argument('-i', nargs='+', action='store', dest='aligner_outputs', required=True, help='''IgBLAST output files in format 7 with query sequence (IgBLAST argument \'-outfmt "7 std qseq sseq btop"\').''') @@ -947,50 +415,112 @@ set of germlines used in the IgBLAST alignment.''') parser_igblast.add_argument('-s', action='store', nargs='+', dest='seq_files', required=True, - help='List of input FASTA files containing sequences') + help='''List of input FASTA files (with .fasta, .fna or .fa + extension), containing sequences.''') parser_igblast.add_argument('--noparse', action='store_true', dest='no_parse', - help='''Specify if input IDs should not be parsed to add - new columns to database.''') - parser_igblast.add_argument('--scores', action='store_true', dest='score_fields', + help='''Specify to prevent input sequence headers from being parsed + to add new columns to database. Parsing of sequence headers requires + headers to be in the pRESTO annotation format, so this should be specified + when sequence headers are incompatible with the pRESTO annotation scheme. + Note, unrecognized header formats will default to this behavior.''') + parser_igblast.add_argument('--partial', action='store_true', dest='partial', + help='''If specified, include incomplete V(D)J alignments in + the pass file instead of the fail file.''') + parser_igblast.add_argument('--scores', action='store_true', dest='parse_scores', help='''Specify if alignment score metrics should be included in the output. Adds the V_SCORE, V_IDENTITY, V_EVALUE, V_BTOP, J_SCORE, J_IDENTITY, J_BTOP, and J_EVALUE columns.''') - parser_igblast.add_argument('--regions', action='store_true', dest='region_fields', - help='''Specify if IMGT framework and CDR regions should be + parser_igblast.add_argument('--regions', action='store_true', dest='parse_regions', + help='''Specify if IMGT FWR and CDRs should be included in the output. Adds the FWR1_IMGT, FWR2_IMGT, FWR3_IMGT, FWR4_IMGT, CDR1_IMGT, CDR2_IMGT, and CDR3_IMGT columns.''') - + parser_igblast.add_argument('--cdr3', action='store_true', + dest='parse_igblast_cdr3', + help='''Specify if the CDR3 sequences generated by IgBLAST + should be included in the output. Adds the columns + CDR3_IGBLAST_NT and CDR3_IGBLAST_AA. Requires IgBLAST + version 1.5 or greater.''') + parser_igblast.set_defaults(func=parseIgBLAST) + # IMGT aligner - parser_imgt = subparsers.add_parser('imgt', help='Process IMGT/HighV-Quest output', - parents=[parser_parent], - formatter_class=CommonHelpFormatter) - imgt_arg_group = parser_imgt.add_mutually_exclusive_group(required=True) - imgt_arg_group.add_argument('-i', nargs='+', action='store', dest='aligner_files', - help='''Either zipped IMGT output files (.zip) or a folder - containing unzipped IMGT output files (which must - include 1_Summary, 2_IMGT-gapped, 3_Nt-sequences, - and 6_Junction).''') + parser_imgt = subparsers.add_parser('imgt', parents=[parser_parent], + formatter_class=CommonHelpFormatter, + help='''Process IMGT/HighV-Quest output + (does not work with V-QUEST).''', + description='''Process IMGT/HighV-Quest output + (does not work with V-QUEST).''') + parser_imgt.add_argument('-i', nargs='+', action='store', dest='aligner_outputs', + help='''Either zipped IMGT output files (.zip or .txz) or a + folder containing unzipped IMGT output files (which must + include 1_Summary, 2_IMGT-gapped, 3_Nt-sequences, + and 6_Junction).''') parser_imgt.add_argument('-s', nargs='*', action='store', dest='seq_files', required=False, - help='List of input FASTA files containing sequences') + help='''List of input FASTA files (with .fasta, .fna or .fa + extension) containing sequences.''') parser_imgt.add_argument('--noparse', action='store_true', dest='no_parse', - help='''Specify if input IDs should not be parsed to add new - columns to database.''') - parser_imgt.add_argument('--scores', action='store_true', dest='score_fields', + help='''Specify to prevent input sequence headers from being parsed + to add new columns to database. Parsing of sequence headers requires + headers to be in the pRESTO annotation format, so this should be specified + when sequence headers are incompatible with the pRESTO annotation scheme. + Note, unrecognized header formats will default to this behavior.''') + parser_imgt.add_argument('--partial', action='store_true', dest='partial', + help='''If specified, include incomplete V(D)J alignments in + the pass file instead of the fail file.''') + parser_imgt.add_argument('--scores', action='store_true', dest='parse_scores', help='''Specify if alignment score metrics should be included in the output. Adds the V_SCORE, V_IDENTITY, - J_SCORE and J_IDENTITY. Note, this will also add - the columns V_EVALUE, V_BTOP, J_EVALUE and J_BTOP, - but they will be empty for IMGT output.''') - parser_imgt.add_argument('--regions', action='store_true', dest='region_fields', - help='''Specify if IMGT framework and CDR regions should be + J_SCORE and J_IDENTITY.''') + parser_imgt.add_argument('--regions', action='store_true', dest='parse_regions', + help='''Specify if IMGT FWRs and CDRs should be included in the output. Adds the FWR1_IMGT, FWR2_IMGT, FWR3_IMGT, FWR4_IMGT, CDR1_IMGT, CDR2_IMGT, and CDR3_IMGT columns.''') + parser_imgt.add_argument('--junction', action='store_true', dest='parse_junction', + help='''Specify if detailed junction fields should be + included in the output. Adds the columns + N1_LENGTH, N2_LENGTH, P3V_LENGTH, P5D_LENGTH, P3D_LENGTH, + P5J_LENGTH, D_FRAME.''') parser_imgt.set_defaults(func=parseIMGT) + # iHMMuneAlign Aligner + parser_ihmm = subparsers.add_parser('ihmm', parents=[parser_parent], + formatter_class=CommonHelpFormatter, + help='Process iHMMune-Align output.', + description='Process iHMMune-Align output.') + parser_ihmm.add_argument('-i', nargs='+', action='store', dest='aligner_outputs', + required=True, + help='''iHMMune-Align output file.''') + parser_ihmm.add_argument('-r', nargs='+', action='store', dest='repo', required=True, + help='''List of folders and/or FASTA files containing + IMGT-gapped germline sequences corresponding to the + set of germlines used in the IgBLAST alignment.''') + parser_ihmm.add_argument('-s', action='store', nargs='+', dest='seq_files', + required=True, + help='''List of input FASTA files (with .fasta, .fna or .fa + extension) containing sequences.''') + parser_ihmm.add_argument('--noparse', action='store_true', dest='no_parse', + help='''Specify to prevent input sequence headers from being parsed + to add new columns to database. Parsing of sequence headers requires + headers to be in the pRESTO annotation format, so this should be specified + when sequence headers are incompatible with the pRESTO annotation scheme. + Note, unrecognized header formats will default to this behavior.''') + parser_ihmm.add_argument('--partial', action='store_true', dest='partial', + help='''If specified, include incomplete V(D)J alignments in + the pass file instead of the fail file.''') + parser_ihmm.add_argument('--scores', action='store_true', dest='parse_scores', + help='''Specify if alignment score metrics should be + included in the output. Adds the path score of the + iHMMune-Align hidden Markov model to HMM_SCORE.''') + parser_ihmm.add_argument('--regions', action='store_true', dest='parse_regions', + help='''Specify if IMGT FWRs and CDRs should be + included in the output. Adds the FWR1_IMGT, FWR2_IMGT, + FWR3_IMGT, FWR4_IMGT, CDR1_IMGT, CDR2_IMGT, and + CDR3_IMGT columns.''') + parser_ihmm.set_defaults(func=parseIHMM) + return parser @@ -1000,7 +530,7 @@ """ parser = getArgParser() args = parser.parse_args() - args_dict = parseCommonArgs(args, in_arg='aligner_files') + args_dict = parseCommonArgs(args, in_arg='aligner_outputs') # Set no ID parsing if sequence files are not provided if 'seq_files' in args_dict and not args_dict['seq_files']: @@ -1008,18 +538,18 @@ # Delete if 'seq_files' in args_dict: del args_dict['seq_files'] - if 'aligner_files' in args_dict: del args_dict['aligner_files'] + if 'aligner_outputs' in args_dict: del args_dict['aligner_outputs'] if 'command' in args_dict: del args_dict['command'] if 'func' in args_dict: del args_dict['func'] if args.command == 'imgt': - for i in range(len(args.__dict__['aligner_files'])): - args_dict['imgt_output'] = args.__dict__['aligner_files'][i] + for i in range(len(args.__dict__['aligner_outputs'])): + args_dict['aligner_output'] = args.__dict__['aligner_outputs'][i] args_dict['seq_file'] = args.__dict__['seq_files'][i] \ if args.__dict__['seq_files'] else None args.func(**args_dict) - elif args.command == 'igblast': - for i in range(len(args.__dict__['aligner_files'])): - args_dict['igblast_output'] = args.__dict__['aligner_files'][i] + elif args.command == 'igblast' or args.command == 'ihmm': + for i in range(len(args.__dict__['aligner_outputs'])): + args_dict['aligner_output'] = args.__dict__['aligner_outputs'][i] args_dict['seq_file'] = args.__dict__['seq_files'][i] args.func(**args_dict)