Makedb number

bin/MakeDb.py

@@ -146,27 +146,32 @@ def correctIMGTFields(receptor, references):
     imgt_dict = {'sequence_imgt': None,
                  'v_germ_start_imgt': None,
                  'v_germ_length_imgt': None,
+                 'v_call': None,
                  'germline_imgt': None}
 
     # Check for necessary fields
     try:
         if not all([receptor.sequence_imgt,
                     receptor.v_germ_start_imgt,
                     receptor.v_germ_length_imgt,
-                    receptor.v_call]):
+                    receptor.v_call,
+                    receptor.germline_imgt]):
             raise AttributeError
     except AttributeError:
         return None
 
     # Gap V region
     try:
-        gapped = gapV(receptor.sequence_imgt,
+        indels,gapped = gapV(receptor.sequence_imgt,
                       receptor.v_germ_start_imgt,
                       receptor.v_germ_length_imgt,
                       receptor.v_call,
-                      references)
+                      references,
+                      remove_indels=True,
+                      germline_alignment = receptor.germline_imgt)
     except KeyError as e:
-        printWarning(e)
+        raise KeyError(e)
+        #printWarning(e)
         return None
 
     # Verify IMGT-gapped sequence and junction concur
@@ -179,18 +184,17 @@ def correctIMGTFields(receptor, references):
 
     # Rebuild germline sequence
     receptor.setDict(gapped, parse=False)
-    __, germlines, __ = buildGermline(receptor, references)
-    # log, germlines, genes = buildGermline(receptor, references)
-    # print(log)
+    #__, germlines, __ = buildGermline(receptor, references)
+    log, germlines, __ = buildGermline(receptor, references)
     if germlines is not None:
         gapped['germline_imgt'] = germlines['full']
     else:
-        return None
+        return None, log, None
 
     # Update return object
     imgt_dict.update(gapped)
 
-    return imgt_dict
+    return imgt_dict, log, indels
 
 
 def getSeqDict(seq_file):
@@ -698,17 +702,19 @@ def parseIHMM(aligner_file, seq_file, repo, cellranger_file=None, validate='stri
 
 
 
-def numberAIRR(aligner_file, repo=None, format=default_format,
-               out_file=None, out_args=default_out_args, nproc=None):
+def numberAIRR(aligner_file, germline_reference=None, format=default_format,
+               out_file=None, out_args=default_out_args,
+               debug_mode=False, nproc=None):
     """
     Inserts IMGT numbering into V fields
 
     Arguments:
       aligner_file (str): AIRR Rearrangement file from the alignment tool.
-      repo (str): folder with germline repertoire files. If None, do not updated alignment columns with IMGT gaps.
+      germline_reference (str): folder with germline repertoire files. If None, do not updated alignment columns with IMGT gaps.
       format (str): output format.
       out_file (str): output file name. Automatically generated from the input file if None.
       out_args (dict): common output argument dictionary from parseCommonArgs.
+      debug_mode (bool): if True, output debug information to log file.
       nproc (int): number of processes to use for parallel processing.
 
     Returns:
@@ -745,7 +751,7 @@ def numberAIRR(aligner_file, repo=None, format=default_format,
         printError(e)
 
     # Load references
-    reference_dict = readGermlines(repo)
+    reference_dict = readGermlines(germline_reference)
 
     # Check for IMGT-gaps in germlines
     if all('...' not in x for x in reference_dict.values()):
@@ -775,7 +781,8 @@ def numberAIRR(aligner_file, repo=None, format=default_format,
         printProgress(rec_count, result_count, 0.05, start_time=start_time)
         rec_count += 1
         # Update IMGT fields
-        imgt_dict = correctIMGTFields(rec, reference_dict)
+        imgt_dict, error_log, indels = correctIMGTFields(rec, reference_dict)
+
         # Write records
         if imgt_dict is not None:
             pass_count += 1
@@ -793,10 +800,19 @@ def numberAIRR(aligner_file, repo=None, format=default_format,
                                 ('D_CALL', rec.d_call),
                                 ('J_CALL', rec.j_call),
                                 ('PRODUCTIVE', rec.functional)])
+            if indels:
+                log['INDELS'] = 'Found %s indels (I:insertion, D:deletion). The insertions were removed from the sequence alignment to be consistent with IMGT positions:' % len(indels) + ','.join(['(%s%s:%s)' % (p,b,t) for p,b,t in indels])
+            if not imgt_dict:
+                log['ERROR'] = error_log['ERROR']
+            if (not imgt_dict) or debug_mode:
+                log['SEQUENCE'] = error_log['SEQUENCE']
+                log['GERMLINE'] = error_log['GERMLINE']
+                log['REGIONS'] = error_log['REGIONS']
+
             printLog(log, log_handle)
 
     # Print counts
-    printProgress(rec_count, result_count, 0.05, start_time=start_time)
+    #printProgress(rec_count, result_count, 0.05, start_time=start_time)
     log = OrderedDict()
     log['OUTPUT'] = os.path.basename(pass_handle.name)
     log['RECORDS'] = rec_count
@@ -1070,23 +1086,23 @@ def getArgParser():
     parser_ihmm.set_defaults(func=parseIHMM, validate='strict')
 
     # Subparser to normalize AIRR file with IMGT-numbering
-    # desc_number = dedent('''
-    #                      Inserts IMGT numbering spacers into sequence_alignment, rebuilds the germline sequence
-    #                      in germline_alignment, and adjusts the values in the coordinate fields v_germline_start
-    #                      and v_germline_end accordingly.
-    #                      ''')
-    # parser_number = subparsers.add_parser('number', parents=[parser_parent],
-    #                                       formatter_class=CommonHelpFormatter, add_help=False,
-    #                                       help='Add IMGT-numbering to an AIRR Rearrangement TSV.',
-    #                                       description=desc_number)
-    # group_number = parser_number.add_argument_group('aligner parsing arguments')
-    # group_number.add_argument('-i', nargs='+', action='store', dest='aligner_files', required=True,
-    #                         help='''AIRR Rearrangement TSV files.''')
-    # group_number.add_argument('-r', nargs='+', action='store', dest='repo', required=False,
-    #                         help='''List of folders and/or fasta files containing
-    #                              IMGT-numbered germline sequences corresponding to the
-    #                              set of germlines used for the alignment.''')
-    # parser_number.set_defaults(func=numberAIRR)
+    desc_number = dedent('''
+                         Inserts IMGT numbering spacers into sequence_alignment, rebuilds the germline sequence
+                         in germline_alignment, and adjusts the values in the coordinate fields v_germline_start
+                         and v_germline_end accordingly.
+                         ''')
+    parser_number = subparsers.add_parser('number', parents=[parser_parent],
+                                          formatter_class=CommonHelpFormatter, add_help=False,
+                                          help='Add IMGT-numbering to an AIRR Rearrangement TSV.',
+                                          description=desc_number)
+    group_number = parser_number.add_argument_group('aligner parsing arguments')
+    group_number.add_argument('-i', nargs='+', action='store', dest='aligner_files', required=True,
+                            help='''AIRR Rearrangement TSV files.''')
+    group_number.add_argument('-g', nargs='+', action='store', dest='germline_reference', required=False,
+                            help='''List of folders and/or fasta files containing
+                                 IMGT-numbered germline sequences corresponding to the
+                                 set of germlines used for the alignment.''')
+    parser_number.set_defaults(func=numberAIRR)
 
     return parser
 

changeo/Alignment.py

@@ -8,9 +8,11 @@
 # Imports
 import re
 from Bio.Seq import Seq
+from Bio import Align
 
 # Presto and changeo imports
 from changeo.Gene import getVAllele, getJAllele
+from presto.Sequence import getDNAScoreDict
 
 # Load regions
 # import yaml
@@ -250,8 +252,49 @@ def alignmentPositions(alignment):
 
     return result
 
+def find_indels(sequence_alignment, germline_alignment):
+    """
+    Identify indels between a sequence alignment and a germline alignment.
 
-def gapV(seq, v_germ_start, v_germ_length, v_call, references, asis_calls=False):
+    Arguments:
+      sequence_alignment (str): Aligned sequence.
+      germline_alignment (str): Aligned germline sequence.
+    Returns:
+      list: List of indels as tuples (position, base, type).
+    """
+    indels = []
+    ref_pos = 1
+    query_pos = 1
+    for ref_base, query_base in zip(germline_alignment, sequence_alignment):
+        if ref_base == '-' and query_base != '-':
+            indels.append((ref_pos, query_base, "I"))
+        elif ref_base != '-' and query_base == '-':
+            indels.append((ref_pos, ref_base, "D"))
+        ref_pos += 1
+        query_pos += 1
+
+    return indels
+
+def remove_insertions(sequence_alignment, indels):
+    """
+    Remove insertions from a sequence alignment based on identified indels.
+
+    Arguments:
+      sequence_alignment (str): Aligned sequence.
+      indels (list): List of indels as tuples (position, base, type). The indel type can be I (insertion) or D (deletion).
+
+    Returns:
+      str: Sequence alignment with insertions removed.
+    """
+    seq_imgt = sequence_alignment
+    offset = -1 # Offset starts with -1 since positions are 1-based, and python indexing is 0-based
+    for indel_pos, indel_base, indel_type in indels:
+        if indel_type == "I":
+            seq_imgt = seq_imgt[:indel_pos + offset] + seq_imgt[indel_pos + offset + 1:]
+            offset -= 1
+    return seq_imgt
+
+def gapV(seq, v_germ_start, v_germ_length, v_call, references, asis_calls=False, remove_indels=False, germline_alignment=None):
     """
     Construction IMGT-gapped V segment sequences.
 
@@ -272,10 +315,12 @@ def gapV(seq, v_germ_start, v_germ_length, v_call, references, asis_calls=False)
     # Initialize return object
     imgt_dict = {'sequence_imgt': None,
                  'v_germ_start_imgt': None,
-                 'v_germ_length_imgt': None}
+                 'v_germ_length_imgt': None,
+                 'v_call': None}
 
     # Initialize imgt gapped sequence
-    seq_imgt = '.' * (int(v_germ_start) - 1) + seq
+    seq_imgt = seq
+    #print("\n" + seq_imgt)
 
     # Extract first V call
     if not asis_calls:
@@ -287,30 +332,50 @@ def gapV(seq, v_germ_start, v_germ_length, v_call, references, asis_calls=False)
     try:
     #if vgene in references:
         vgap = references[vgene]
-        # Iterate over gaps in the germline segment
-        gaps = re.finditer(r'\.', vgap)
-        gapcount = int(v_germ_start) - 1
-        for gap in gaps:
-            i = gap.start()
-            # Break if gap begins after V region
-            if i >= v_germ_length + gapcount:
-                break
-            # Insert gap into IMGT sequence
-            seq_imgt = seq_imgt[:i] + '.' + seq_imgt[i:]
-            # Update gap counter
-            gapcount += 1
-
-        imgt_dict['sequence_imgt'] = seq_imgt
-        # Update IMGT positioning information for V
-        imgt_dict['v_germ_start_imgt'] = 1
-        imgt_dict['v_germ_length_imgt'] = v_germ_length + gapcount
     except KeyError as e:
         raise KeyError('%s was not found in the germline repository.' % vgene)
-    #else:
-    #    printWarning('%s was not found in the germline repository. IMGT-gapped sequence cannot be determined.' % vgene)
-
-    return imgt_dict
 
+    # We need to check for indels first
+    # If aligned germline is provided use that to identify indels
+    # Otherwise perform pairwise alignment with reference
+    indels = []
+    if remove_indels:
+        if not germline_alignment:
+            raise ValueError('Aligned germline sequence must be provided when remove_indels is True.')
+
+        ungapped_reference = germline_alignment
+        #print(ungapped_reference)
+
+        if seq_imgt.count('-') > 0 or ungapped_reference.count('-') > 0:
+            indels = find_indels(seq, ungapped_reference)
+            seq_imgt = remove_insertions(seq, indels)
+
+    # Iterate over gaps in the germline segment
+    seq_imgt = '.' * (int(v_germ_start) - 1) + seq_imgt
+
+    gaps = re.finditer(r'\.', vgap)
+    gapcount = int(v_germ_start) - 1
+    for gap in gaps:
+        i = gap.start()
+        # Break if gap begins after V region
+        if i >= v_germ_length + gapcount:
+            break
+        # Insert gap into IMGT sequence
+        seq_imgt = seq_imgt[:i] + '.' + seq_imgt[i:]
+        # Update gap counter
+        gapcount += 1
+
+
+
+    # Assign IMGT gapped sequence
+    imgt_dict['sequence_imgt'] = seq_imgt
+    # Update IMGT positioning information for V
+    imgt_dict['v_germ_start_imgt'] = 1
+    imgt_dict['v_germ_length_imgt'] = v_germ_length + gapcount
+    imgt_dict['v_call'] = v_call
+
+    # Return indels if present
+    return (indels if indels else None), imgt_dict
 
 def inferJunction(seq, j_germ_start, j_germ_length, j_call, references, asis_calls=False, regions='default'):
     """

changeo/Applications.py

@@ -181,6 +181,7 @@ def runIgBLASTN(fasta, igdata, loci='ig', organism='human', vdb=None, ddb=None,
            '-germline_db_D', str(d_germ),
            '-germline_db_J', str(j_germ),
            '-outfmt', outfmt,
+           '-show_translation',
            '-domain_system', 'imgt']
 
     # Add C-region arguments for igblastn v1.18.0

changeo/Gene.py

@@ -8,6 +8,7 @@
 # Imports
 import re
 from collections import OrderedDict
+from Bio.Seq import Seq
 
 # Presto and changeo imports
 from changeo.Defaults import v_attr, d_attr, j_attr, seq_attr
@@ -301,6 +302,13 @@ def getDGermline(receptor, references, d_field=d_attr, amino_acid=False):
     elif dgene in references:
         # Define D germline sequence
         dseq = references[dgene]
+        # Handle D gene inversions
+        if receptor.d_germ_start > receptor.d_germ_end:
+            dlen = receptor.d_germ_start - receptor.d_germ_end + 1
+            reverse_complement = Seq(dseq).reverse_complement()
+            dseq = str(reverse_complement)
+            dstart = receptor.d_germ_end -1
+
         germ_dseq = dseq[dstart:(dstart + dlen)]
     else:
         germ_dseq = None
@@ -509,6 +517,7 @@ def buildGermline(receptor, references, seq_field=seq_attr, v_field=v_attr,
     if germ_dseq is None:
         log['ERROR'] = 'Allele %s is not in the provided germline database.' % dgene
         return log, None, None
+    log['D_GENE'] = germ_dseq
 
     # Build J segment germline sequence
     jgene, germ_jseq = getJGermline(receptor, references, j_field=j_field, amino_acid=amino_acid)

changeo/IO.py

@@ -1462,7 +1462,7 @@ def parseSections(self, sections):
         # Create IMGT-gapped sequence
         if ('v_call' in db and db['v_call']) and ('sequence_trim' in db and db['sequence_trim']):
             try:
-                imgt_dict = gapV(db['sequence_trim'],
+                _, imgt_dict = gapV(db['sequence_trim'],
                                  v_germ_start=db['v_germ_start_vdj'],
                                  v_germ_length=db['v_germ_length_vdj'],
                                  v_call=db['v_call'],
@@ -1676,7 +1676,7 @@ def parseSections(self, sections):
         # Create IMGT-gapped sequence
         if ('v_call' in db and db['v_call']) and ('sequence_aa_trim' in db and db['sequence_aa_trim']):
             try:
-                gap = gapV(db['sequence_aa_trim'],
+                _, gap = gapV(db['sequence_aa_trim'],
                            v_germ_start=db['v_germ_aa_start_vdj'],
                            v_germ_length=db['v_germ_aa_length_vdj'],
                            v_call=db['v_call'],
@@ -2132,7 +2132,7 @@ def parseRecord(self, record):
         # Create IMGT-gapped sequence
         if 'v_call' in db and db['v_call'] and 'sequence_vdj' in db and db['sequence_vdj']:
             try:
-                imgt_dict = gapV(db['sequence_vdj'],
+                _, imgt_dict = gapV(db['sequence_vdj'],
                                  v_germ_start=db['v_germ_start_vdj'],
                                  v_germ_length=db['v_germ_length_vdj'],
                                  v_call=db['v_call'],

changeo/Receptor.py

@@ -1390,7 +1390,8 @@ def d_germ_end(self):
         """
         Position of the last nucleotide in the D germline sequence alignment
         """
-        try:  return self.d_germ_start + self.d_germ_length - 1
+        try:  
+            return self.d_germ_start + self.d_germ_length - 1
         except TypeError:  return None
 
     @property