Document and validate that CreateGermlines only supports nucleotide sequences

bin/CreateGermlines.py

@@ -1,6 +1,8 @@
 #!/usr/bin/env python3
 """
-Reconstructs germline sequences from alignment data
+Reconstructs germline sequences from nucleotide alignment data
+
+Note: Amino acid sequences are not currently supported.
 """
 
 # Info
@@ -40,7 +42,7 @@ def createGermlines(db_file, references, seq_field=default_seq_field, v_field=de
     Arguments:
       db_file : input tab-delimited database file.
       references : folders and/or files containing germline repertoire data in FASTA format.
-      seq_field : field in which to look for sequence.
+      seq_field : field in which to look for sequence. Must be a nucleotide sequence field.
       v_field : field in which to look for V call.
       d_field : field in which to look for D call.
       j_field : field in which to look for J call.
@@ -53,6 +55,9 @@ def createGermlines(db_file, references, seq_field=default_seq_field, v_field=de
 
     Returns:
       dict: names of the 'pass' and 'fail' output files.
+
+    Note:
+      Amino acid sequences are not currently supported. The seq_field must contain nucleotide sequences.
     """
     # Print parameter info
     log = OrderedDict()
@@ -115,6 +120,32 @@ def createGermlines(db_file, references, seq_field=default_seq_field, v_field=de
         if f not in db_iter.fields:
             printError('%s field does not exist in input database file.' % f)
 
+    # Check for amino acid sequences
+    aa_fields = ['sequence_aa', 'sequence_alignment_aa', 'SEQUENCE_AA', 'SEQUENCE_ALIGNMENT_AA']
+    if seq_field in aa_fields:
+        printError('Reconstruction of germlines from amino acid sequences is not currently supported. '
+                   'Please use a nucleotide sequence field (e.g., sequence_alignment or SEQUENCE_IMGT).')
+
+    # Also check actual sequence content to detect amino acids
+    # Sample first few sequences to detect amino acid-specific characters
+    aa_chars = set('EFIKLMPQWY')  # Amino acid characters not valid in nucleotide sequences
+    sample_count = 0
+    max_samples = 10
+    for rec in db_iter:
+        seq = rec.getField(seq_field)
+        if seq and any(c in aa_chars for c in seq.upper().replace('.', '').replace('-', '')):
+            printError('Amino acid sequence detected in %s field. '
+                       'Reconstruction of germlines from amino acid sequences is not currently supported. '
+                       'Please use a nucleotide sequence field.' % seq_field)
+        sample_count += 1
+        if sample_count >= max_samples:
+            break
+
+    # Reset the file handle after sampling
+    db_handle.close()
+    db_handle = open(db_file, 'rt')
+    db_iter = reader(db_handle)
+
     # Translate to Receptor attribute names
     v_field = schema.toReceptor(v_field)
     d_field = schema.toReceptor(d_field)
@@ -299,7 +330,8 @@ def getArgParser():
                              germline call used for the entire clone within the
                              germline_v_call, germline_d_call and germline_j_call fields.''')
     group.add_argument('--sf', action='store', dest='seq_field', default=None,
-                        help='''Field containing the aligned sequence.
+                        help='''Field containing the aligned sequence. Must be a nucleotide sequence field.
+                             Amino acid sequences are not currently supported.
                              Defaults to sequence_alignment (airr) or SEQUENCE_IMGT (changeo).''')
     group.add_argument('--vf', action='store', dest='v_field', default=None,
                         help='''Field containing the germline V segment call.

docs/examples/germlines.rst

@@ -60,9 +60,16 @@ argument is recommended, as this will generate a single germline of consensus le
     there is insufficient information to successfully reconstruct germline sequences
     for all cases using :ref:`CreateGermlines`.
 
+.. warning::
+
+    :ref:`CreateGermlines` only supports reconstruction from **nucleotide sequences**.
+    Amino acid sequences are not currently supported. Ensure that the sequence field 
+    (e.g., ``sequence_alignment`` or ``SEQUENCE_IMGT``) contains nucleotide data, not 
+    amino acid data (e.g., ``sequence_alignment_aa`` or ``SEQUENCE_AA``).
+
 .. seealso::
 
-    The `TIgGER <https://github.com/immcantation/tigger>`__ R package provided tools for
+    The `TIgGER <https://github.com/immcantation/tigger>`__ R package provides tools for
     identifying novel polymorphisms and building a personalized germline database. To
     use the germline corrections provided by `TIgGER <https://github.com/immcantation/tigger>`__
     you would replace the V-segment germline file with the one generated by

docs/methods/germlines.rst

@@ -24,3 +24,11 @@ single V and J allele to use as the germline from all the assigned
 annotations in each clone. The selection is made by simple majority rule of all 
 the allele calls in the clone. After the germline reconstruction process, all 
 sequences belonging to the same clone have been assigned the same germline.
+
+Limitations
+--------------------------------------------------------------------------------
+
+:ref:`CreateGermlines` currently only supports reconstruction from **nucleotide sequences**.
+Amino acid sequences are not supported at this time. The tool requires nucleotide-based 
+positional fields (e.g., ``v_germ_start_imgt``, ``d_germ_start``, ``j_germ_start``) and 
+nucleotide sequence data to perform accurate germline reconstruction.