diff --git a/bin/CreateGermlines.py b/bin/CreateGermlines.py index 8fd730b..2bf10cf 100755 --- a/bin/CreateGermlines.py +++ b/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. diff --git a/docs/examples/germlines.rst b/docs/examples/germlines.rst index 455a465..ea974db 100644 --- a/docs/examples/germlines.rst +++ b/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 `__ R package provided tools for + The `TIgGER `__ R package provides tools for identifying novel polymorphisms and building a personalized germline database. To use the germline corrections provided by `TIgGER `__ you would replace the V-segment germline file with the one generated by diff --git a/docs/methods/germlines.rst b/docs/methods/germlines.rst index f46acfa..ab36240 100644 --- a/docs/methods/germlines.rst +++ b/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.