From 3ebb93199abe698bc15fa58fea0483c0d84aaedc Mon Sep 17 00:00:00 2001 From: John Orgera <65687576+johnoooh@users.noreply.github.com> Date: Wed, 26 Aug 2026 16:38:07 -0400 Subject: [PATCH 1/2] feat: add ANNOTATE_HLAHD module Post-processes HLA-HD class I output (final.result.txt + per-locus est.txt) into a P-group-annotated TSV, with an optional self-contained HTML report (pass --skip_html via task.ext.args to omit it). - Container: ghcr.io/mskcc-omics-workflows/hlahd-tools:1.0.0 (companion PR: mskcc-omics-workflows/containers#85) - Script vendored from mskcc/HLA_HD_workflow's scripts/annotate_hlahd.py + scripts/hlahd_annotate/ into resources/usr/bin/ (private repo, so it can't be pulled at container-build time -- same pattern as neoantigen-utils-base before neoantigen-utils had its own repo) - Test data: mskcc-omics-workflows/test-datasets@feature/annotate_hlahd (synthetic fixtures; pending Review Team promotion to an official 'annotate_hlahd' branch) - Tests: real-data run, --skip_html run (asserts the report output is empty), and stub -- all passing locally against the built image --- modules/msk/annotate_hlahd/environment.yml | 7 + modules/msk/annotate_hlahd/main.nf | 50 ++++++ modules/msk/annotate_hlahd/meta.yml | 64 +++++++ .../resources/usr/bin/annotate_hlahd.py | 67 +++++++ .../usr/bin/hlahd_annotate/__init__.py | 0 .../usr/bin/hlahd_annotate/annotate.py | 107 ++++++++++++ .../usr/bin/hlahd_annotate/parse_est.py | 163 ++++++++++++++++++ .../usr/bin/hlahd_annotate/parse_final.py | 62 +++++++ .../usr/bin/hlahd_annotate/pgroup.py | 71 ++++++++ .../usr/bin/hlahd_annotate/report.py | 107 ++++++++++++ .../usr/bin/templates/report.html.j2 | 135 +++++++++++++++ modules/msk/annotate_hlahd/tests/main.nf.test | 115 ++++++++++++ .../annotate_hlahd/tests/main.nf.test.snap | 62 +++++++ .../msk/annotate_hlahd/tests/nextflow.config | 5 + modules/msk/annotate_hlahd/tests/tags.yml | 2 + tests/config/test_data.config | 8 + 16 files changed, 1025 insertions(+) create mode 100644 modules/msk/annotate_hlahd/environment.yml create mode 100644 modules/msk/annotate_hlahd/main.nf create mode 100644 modules/msk/annotate_hlahd/meta.yml create mode 100755 modules/msk/annotate_hlahd/resources/usr/bin/annotate_hlahd.py create mode 100644 modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/__init__.py create mode 100644 modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/annotate.py create mode 100644 modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/parse_est.py create mode 100644 modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/parse_final.py create mode 100644 modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/pgroup.py create mode 100644 modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/report.py create mode 100644 modules/msk/annotate_hlahd/resources/usr/bin/templates/report.html.j2 create mode 100644 modules/msk/annotate_hlahd/tests/main.nf.test create mode 100644 modules/msk/annotate_hlahd/tests/main.nf.test.snap create mode 100644 modules/msk/annotate_hlahd/tests/nextflow.config create mode 100644 modules/msk/annotate_hlahd/tests/tags.yml diff --git a/modules/msk/annotate_hlahd/environment.yml b/modules/msk/annotate_hlahd/environment.yml new file mode 100644 index 00000000..da913a15 --- /dev/null +++ b/modules/msk/annotate_hlahd/environment.yml @@ -0,0 +1,7 @@ +--- +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/environment-schema.json +channels: + - conda-forge + - bioconda +dependencies: + - "ANNOTATE_HLAHD=HERE" diff --git a/modules/msk/annotate_hlahd/main.nf b/modules/msk/annotate_hlahd/main.nf new file mode 100644 index 00000000..b9abff28 --- /dev/null +++ b/modules/msk/annotate_hlahd/main.nf @@ -0,0 +1,50 @@ +process ANNOTATE_HLAHD { + tag "$meta.id" + label 'process_single' + + conda "${moduleDir}/environment.yml" + container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ? + 'docker://ghcr.io/mskcc-omics-workflows/hlahd-tools:1.0.0': + 'ghcr.io/mskcc-omics-workflows/hlahd-tools:1.0.0' }" + + input: + tuple val(meta), path(result_dir) + path pgroup_file + + output: + tuple val(meta), path("${prefix}_annotated.tsv"), emit: tsv + tuple val(meta), path("${prefix}_report.html"), emit: report, optional: true + path "versions.yml", emit: versions + + when: + task.ext.when == null || task.ext.when + + script: + def args = task.ext.args ?: '' + prefix = task.ext.prefix ?: "${meta.id}" + """ + annotate_hlahd.py \\ + --result_dir ${result_dir} \\ + --sample ${prefix} \\ + --pgroup_file ${pgroup_file} \\ + --outdir . \\ + ${args} + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + python: \$(python3 --version | sed 's/Python //') + END_VERSIONS + """ + + stub: + prefix = task.ext.prefix ?: "${meta.id}" + """ + echo -e "locus\\tallele1\\tallele2\\tp_group" > ${prefix}_annotated.tsv + echo "stub report" > ${prefix}_report.html + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + python: 3.11 + END_VERSIONS + """ +} diff --git a/modules/msk/annotate_hlahd/meta.yml b/modules/msk/annotate_hlahd/meta.yml new file mode 100644 index 00000000..b0f63aef --- /dev/null +++ b/modules/msk/annotate_hlahd/meta.yml @@ -0,0 +1,64 @@ +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/meta-schema.json +name: "annotate_hlahd" +description: Annotate HLA-HD class I output with IMGT P-groups and quality flags +keywords: + - hla + - hlahd + - annotation + - immunogenomics +tools: + - "annotate_hlahd": + description: "Post-processes HLA-HD's per-sample final.result.txt/*.est.txt class I output, mapping each allele call to its IMGT P-group and flagging low-confidence calls." + homepage: "https://github.com/mskcc/HLA_HD_workflow" + documentation: "https://github.com/mskcc/HLA_HD_workflow" + licence: + - "MIT" + identifier: "" +input: + - - meta: + type: map + description: | + Groovy Map containing sample information + e.g. `[ id:'sample1' ]` + - result_dir: + type: directory + description: Directory containing HLA-HD's _final.result.txt and _{A,B,C}.est.txt + pattern: "*" + - - pgroup_file: + type: file + description: IMGT wmda/hla_nom_p.txt P-group reference table + pattern: "*.txt" +output: + - tsv: + - meta: + type: map + description: | + Groovy Map containing sample information + e.g. `[ id:'sample1' ]` + - ${prefix}_annotated.tsv: + type: file + description: Per-allele annotated TSV with P-group and quality-flag columns + pattern: "*_annotated.tsv" + - report: + - meta: + type: map + description: | + Groovy Map containing sample information + e.g. `[ id:'sample1' ]` + - ${prefix}_report.html: + type: file + description: | + Self-contained HTML report for the sample. Optional — pass `--skip_html` + via `task.ext.args` to skip generating it. + pattern: "*_report.html" + - versions: + - versions.yml: + type: file + description: File containing software versions + pattern: "versions.yml" + ontologies: + - edam: http://edamontology.org/format_3750 +authors: + - "@johnoooh" +maintainers: + - "@johnoooh" diff --git a/modules/msk/annotate_hlahd/resources/usr/bin/annotate_hlahd.py b/modules/msk/annotate_hlahd/resources/usr/bin/annotate_hlahd.py new file mode 100755 index 00000000..9faa7e87 --- /dev/null +++ b/modules/msk/annotate_hlahd/resources/usr/bin/annotate_hlahd.py @@ -0,0 +1,67 @@ +#!/usr/bin/env python3 +""" +annotate_hlahd.py — Post-process HLA-HD class I output. + +Produces per-sample annotated TSV and self-contained HTML report. + +Usage: + annotate_hlahd.py \ + --result_dir results/SAMPLE1/ \ + --sample SAMPLE1 \ + --pgroup_file /path/to/hla_nom_p.txt \ + --outdir output/ + +Vendored from mskcc/HLA_HD_workflow's scripts/annotate_hlahd.py + scripts/hlahd_annotate/ +for the ANNOTATE_HLAHD module (see modules/msk/annotate_hlahd/resources/usr/bin/). +""" +import argparse +import sys +from pathlib import Path + +# hlahd_annotate/ is a sibling of this file (both live in resources/usr/bin/, +# which Nextflow adds to PATH for any process including this module). +sys.path.insert(0, str(Path(__file__).parent)) + +from hlahd_annotate.annotate import annotate_sample +from hlahd_annotate.report import render_report + +TEMPLATE_DIR = Path(__file__).parent / "templates" + + +def parse_args(): + p = argparse.ArgumentParser( + description='Annotate HLA-HD class I output with P-groups and quality flags.' + ) + p.add_argument('--result_dir', required=True, type=Path, + help='Directory containing _final.result.txt and _{A,B,C}.est.txt') + p.add_argument('--sample', required=True, + help='Sample ID (used as filename prefix for output files)') + p.add_argument('--pgroup_file', required=True, type=Path, + help='IMGT wmda/hla_nom_p.txt P-group reference table') + p.add_argument('--outdir', required=True, type=Path, + help='Output directory for TSV and HTML report (created if absent)') + p.add_argument('--hlahd_version', default='v1.7.1', + help='HLA-HD version string for the report footer (default: v1.7.1)') + p.add_argument('--skip_html', action='store_true', + help='Skip generating the HTML report; write only the annotated TSV') + return p.parse_args() + + +def main(): + args = parse_args() + args.outdir.mkdir(parents=True, exist_ok=True) + + df = annotate_sample(args.result_dir, args.sample, args.pgroup_file) + + tsv_path = args.outdir / f"{args.sample}_annotated.tsv" + df.to_csv(tsv_path, sep='\t', index=False) + print(f"TSV written: {tsv_path}") + + if not args.skip_html: + html_path = render_report(df, args.sample, str(args.pgroup_file), args.outdir, TEMPLATE_DIR, + hlahd_version=args.hlahd_version) + print(f"HTML written: {html_path}") + + +if __name__ == '__main__': + main() diff --git a/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/__init__.py b/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/annotate.py b/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/annotate.py new file mode 100644 index 00000000..2ec85844 --- /dev/null +++ b/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/annotate.py @@ -0,0 +1,107 @@ +import sys +from pathlib import Path +import pandas as pd + +from .pgroup import load_pgroup_table, lookup_pgroup +from .parse_final import parse_final_result, CLASS_I_LOCI +from .parse_est import parse_est, cross_check_allele + +# Map allele_position to (best_pairs index, coverage key within that pair) +_POSITION_TO_COVERAGE = { + 'allele1': (0, 'allele1_coverage'), + 'allele2': (0, 'allele2_coverage'), + 'allele1_pair2': (1, 'allele1_coverage'), + 'allele2_pair2': (1, 'allele2_coverage'), +} + + +def annotate_sample( + result_dir: Path, + sample: str, + pgroup_file: Path, +) -> pd.DataFrame: + """ + Combine parse_final, parse_est, and P-group lookup into the annotated DataFrame. + Emits a WARNING to stderr for each allele that fails est cross-check. + + Returns DataFrame with columns: + sample, locus, allele, allele_position, resolution, + p_group, p_group_found, multiple_best_pairs, has_ambiguous_pair, + ambiguous_alleles, est_mismatch, est_mismatch_detail, + exon2_depth, exon2_incomp, exon3_depth, exon3_incomp, incomplete_coverage + """ + result_dir = Path(result_dir) + final_file = result_dir / f"{sample}_final.result.txt" + + df = parse_final_result(final_file) + pgroup_lookup = load_pgroup_table(pgroup_file) + + # Load est data for each class I locus (returns empty-flag dict if file missing) + est_data: dict[str, dict] = {} + for locus in CLASS_I_LOCI: + est_file = result_dir / f"{sample}_{locus}.est.txt" + est_data[locus] = parse_est(est_file) + + records = [] + for _, row in df.iterrows(): + locus = row['locus'] + allele = row['allele'] + est = est_data.get(locus, {}) + is_pair2 = 'pair2' in row['allele_position'] + + # P-group lookup + p_group, p_group_found = lookup_pgroup(allele, pgroup_lookup) + + # Cross-check (primary pair alleles only; pair2 alleles not cross-checked) + est_mismatch = False + est_mismatch_detail = None + if allele and not is_pair2: + match, detail = cross_check_allele(allele, est.get('best_pair_alleles', [])) + if not match: + est_mismatch = True + est_mismatch_detail = detail + print( + f"WARNING: est_mismatch for sample={sample} locus={locus} " + f"allele_position={row['allele_position']}: {detail}", + file=sys.stderr, + ) + + # multiple_best_pairs: prefer est.txt (more precise) over final.result.txt count + multiple_best_pairs = est.get('multiple_best_pairs', row['multiple_best_pairs']) + + # Per-allele exon coverage from est best_pairs + coverage = {'exon2_depth': None, 'exon2_incomp': None, + 'exon3_depth': None, 'exon3_incomp': None} + best_pairs = est.get('best_pairs', []) + # Unrecognized allele_position falls through to all-None coverage (intended silent-null) + mapping = _POSITION_TO_COVERAGE.get(row['allele_position']) + if mapping: + idx, key = mapping + if idx < len(best_pairs): + coverage = dict(best_pairs[idx][key]) + incomplete_coverage = ( + (coverage['exon2_incomp'] or 0) > 0 + or (coverage['exon3_incomp'] or 0) > 0 + ) + + records.append({ + 'sample': sample, + 'locus': locus, + 'allele': allele, + 'allele_position': row['allele_position'], + 'resolution': row['resolution'], + 'p_group': p_group, + 'p_group_found': p_group_found, + 'multiple_best_pairs': multiple_best_pairs, + 'has_ambiguous_pair': est.get('has_ambiguous_pair', False), + 'ambiguous_alleles': est.get('ambiguous_alleles'), + 'est_mismatch': est_mismatch, + 'est_mismatch_detail': est_mismatch_detail, + 'exon2_depth': coverage['exon2_depth'], + 'exon2_incomp': coverage['exon2_incomp'], + 'exon3_depth': coverage['exon3_depth'], + 'exon3_incomp': coverage['exon3_incomp'], + 'incomplete_coverage': incomplete_coverage, + }) + + return pd.DataFrame(records) diff --git a/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/parse_est.py b/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/parse_est.py new file mode 100644 index 00000000..5dd81946 --- /dev/null +++ b/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/parse_est.py @@ -0,0 +1,163 @@ +from pathlib import Path + +_NULL_ALLELE = {'-', ''} + +_EXONS = ('exon2', 'exon3') + + +def _parse_coverage(field: str) -> dict: + """ + Parse a coverage field like 'exon2:262.426:comp.0,exon3:303.069:incomp.7' + into {exon2_depth, exon2_incomp, exon3_depth, exon3_incomp}. + comp.N -> incomp 0; incomp.N -> incomp N. Missing/null -> None values. + """ + cov = {} + for e in _EXONS: + cov[f'{e}_depth'] = None + cov[f'{e}_incomp'] = None + + if not field or field.strip() in _NULL_ALLELE: + return cov + + for part in field.split(','): + bits = part.split(':') + if len(bits) != 3: + continue + exon, depth, status = bits[0].strip(), bits[1].strip(), bits[2].strip() + if exon not in _EXONS: + continue + try: + cov[f'{exon}_depth'] = float(depth) + except ValueError: + cov[f'{exon}_depth'] = None + if status.startswith('incomp.'): + try: + cov[f'{exon}_incomp'] = int(status.split('.', 1)[1]) + except ValueError: + cov[f'{exon}_incomp'] = None + elif status.startswith('comp'): + cov[f'{exon}_incomp'] = 0 + + return cov + + +def _first_allele(field: str) -> str | None: + """First allele in a comma-separated est column, HLA- stripped; None if null.""" + field = field.strip() + if field in _NULL_ALLELE: + return None + first = field.split(',')[0].strip().removeprefix('HLA-') + return first or None + + +def parse_est(est_file: Path) -> dict: + """ + Parse *_{A,B,C}.est.txt. Used for flag extraction and cross-validation only. + Never used as the authoritative source for allele calls. + + Returns dict with: + multiple_best_pairs (bool): True if #Best allele pair count > 1 + has_ambiguous_pair (bool): True if #Other ambiguous pair section present + ambiguous_alleles (str | None): 'allele1 / allele2' from ambiguous section + best_pair_alleles (list[str]): all alleles from best-pair data rows, HLA- stripped + best_pairs (list[dict]): per best pair, with representative names + allele1/allele2 (str | None) and allele1_coverage/allele2_coverage dicts + """ + result = { + 'multiple_best_pairs': False, + 'has_ambiguous_pair': False, + 'ambiguous_alleles': None, + 'best_pair_alleles': [], + 'best_pairs': [], + } + + if not est_file.exists(): + return result + + with open(est_file) as f: + lines = [line.rstrip('\n') for line in f] + + in_ambiguous_section = False + best_pair_alleles = [] + best_pairs = [] + + for line in lines: + if not line: + continue + + if line.startswith('#Best allele pair'): + # Format: "#Best allele pair\t" + parts = line.split('\t') + count = int(parts[1]) if len(parts) > 1 else 1 + result['multiple_best_pairs'] = count > 1 + in_ambiguous_section = False + + elif line.startswith('#Other ambiguous pair'): + result['has_ambiguous_pair'] = True + in_ambiguous_section = True + + elif line.startswith('#'): + in_ambiguous_section = False + + else: + cols = line.split('\t') + if in_ambiguous_section: + # Ambiguous data lines: single allele per column (not lists) + if len(cols) >= 2: + a1 = cols[0].strip().removeprefix('HLA-') + a2 = cols[1].strip().removeprefix('HLA-') + result['ambiguous_alleles'] = f"{a1} / {a2}" + else: + # Best pair data lines: col0=allele1_comma_list, col1=allele2_comma_list + for col in cols[:2]: + col = col.strip() + if col in _NULL_ALLELE: + continue + for a in col.split(','): + a = a.strip().removeprefix('HLA-') + if a: + best_pair_alleles.append(a) + allele1_cov = _parse_coverage(cols[2]) if len(cols) > 2 else _parse_coverage('') + allele2_cov = _parse_coverage(cols[3]) if len(cols) > 3 else _parse_coverage('') + best_pairs.append({ + 'allele1': _first_allele(cols[0]) if len(cols) > 0 else None, + 'allele2': _first_allele(cols[1]) if len(cols) > 1 else None, + 'allele1_coverage': allele1_cov, + 'allele2_coverage': allele2_cov, + }) + + result['best_pair_alleles'] = best_pair_alleles + result['best_pairs'] = best_pairs + return result + + +def cross_check_allele( + allele: str | None, + best_pair_alleles: list[str], +) -> tuple[bool, str | None]: + """ + Check that an allele from final.result.txt appears in the est.txt best pair allele lists. + Matching is performed at the field depth of the reported allele (2 or 3 fields). + + Returns (match_found, detail_string). + If match_found is True, detail_string is None. + Hard warning to stderr is the caller's responsibility. + """ + if not allele or not best_pair_alleles: + return True, None + + a = allele.removeprefix('HLA-') + fields = a.split(':') + n_fields = len(fields) + + truncated = set() + for bp in best_pair_alleles: + bp_fields = bp.split(':') + truncated.add(':'.join(bp_fields[:n_fields])) + + if a in truncated: + return True, None + + shown = sorted(truncated)[:10] + detail = f"final={allele}; est_best_alleles={','.join(shown)}" + return False, detail diff --git a/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/parse_final.py b/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/parse_final.py new file mode 100644 index 00000000..6564a448 --- /dev/null +++ b/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/parse_final.py @@ -0,0 +1,62 @@ +from pathlib import Path +import warnings +import pandas as pd + +CLASS_I_LOCI = {'A', 'B', 'C'} +_NULL_VALUES = {'Not typed', '-', ''} +# Positions for alleles: primary pair then optional second pair +_PAIR_POSITIONS = ['allele1', 'allele2', 'allele1_pair2', 'allele2_pair2'] + + +def parse_final_result(result_file: Path) -> pd.DataFrame: + """ + Parse *_final.result.txt. Returns only class I loci (A, B, C). + + Returns DataFrame with columns: + locus, allele, allele_position, resolution, multiple_best_pairs + where: + - allele: str with HLA- prefix, or None/NaN if Not typed / - + - allele_position: 'allele1', 'allele2', 'allele1_pair2', 'allele2_pair2' + - resolution: number of colon-separated fields (0 if null) + - multiple_best_pairs: True if row has >2 allele columns + """ + rows = [] + with open(result_file) as f: + for line in f: + line = line.strip() + if not line: + continue + parts = line.split('\t') + locus = parts[0] + if locus not in CLASS_I_LOCI: + continue + + alleles = parts[1:] + multiple_best_pairs = len(alleles) > 2 + + if len(alleles) > len(_PAIR_POSITIONS): + warnings.warn( + f"Row for locus {locus!r} has {len(alleles)} allele columns " + f"(max supported: {len(_PAIR_POSITIONS)}); extra alleles ignored.", + UserWarning, + stacklevel=2, + ) + + for pos, allele in zip(_PAIR_POSITIONS, alleles): + allele_val = None if allele in _NULL_VALUES else allele + if allele_val is not None: + fields = allele_val.removeprefix('HLA-').split(':') + resolution = len(fields) + else: + resolution = 0 + + rows.append({ + 'locus': locus, + 'allele': allele_val, + 'allele_position': pos, + 'resolution': resolution, + 'multiple_best_pairs': multiple_best_pairs, + }) + + return pd.DataFrame(rows, columns=['locus', 'allele', 'allele_position', + 'resolution', 'multiple_best_pairs']) diff --git a/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/pgroup.py b/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/pgroup.py new file mode 100644 index 00000000..a0079729 --- /dev/null +++ b/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/pgroup.py @@ -0,0 +1,71 @@ +from pathlib import Path + + +def load_pgroup_table(pgroup_file: Path) -> dict[str, str]: + """ + Parse IMGT wmda hla_nom_p.txt. + Returns lookup dict: allele at 2/3/4-field resolution -> P-group designation. + e.g., 'A*02:01' -> 'A*02:01P', 'A*02:01:01' -> 'A*02:01P' + """ + lookup: dict[str, str] = {} + with open(pgroup_file) as f: + for line in f: + line = line.strip() + if not line or line.startswith('#'): + continue + parts = line.split(';') + if len(parts) != 3: + continue + locus_prefix = parts[0] # e.g., 'A*' + alleles_str = parts[1] # e.g., '02:01:01:01/02:01:01:02' + pg_raw = parts[2].strip() + if not pg_raw: + continue # skip alleles not assigned to any P-group + p_group = locus_prefix + pg_raw # e.g., 'A*02:01P' + + for allele_fields in alleles_str.split('/'): + allele_fields = allele_fields.strip() + if not allele_fields: + continue + # Strip trailing IMGT expression suffixes (N=null, L=low, S=secreted, + # Q=questionable, C=aberrant cytoplasm, A=aberrant, G=null genomic) + # from the last colon-field. IMGT uses uppercase only. + fields = allele_fields.split(':') + fields[-1] = fields[-1].rstrip('NLSQCAG') + # MAX_HLA_FIELDS = 4; build keys at 2-, 3-, and 4-field resolution + for n in range(2, min(5, len(fields) + 1)): + key = locus_prefix + ':'.join(fields[:n]) + if key not in lookup: + lookup[key] = p_group + return lookup + + +def lookup_pgroup(allele: str | None, lookup: dict[str, str]) -> tuple[str | None, bool]: + """ + Look up P-group for a reported allele. + + Strips HLA- prefix before lookup. Tries match at reported field depth, + then falls back to shorter fields (minimum 2 fields). + + Sentinel values treated as "not found" (returns (None, False)): + - None + - 'Not typed' (HLA-HD value when locus has insufficient reads) + - '-' (HLA-HD value when only one allele is identified) + + Returns: + (p_group, found): tuple of the P-group string and a bool indicating success. + If not found, p_group is None and found is False. + """ + if not allele or allele in ('Not typed', '-'): + return None, False + + a = allele.removeprefix('HLA-') + fields = a.split(':') + + # Try from full depth down to 2 fields + for n in range(len(fields), 1, -1): + key = ':'.join(fields[:n]) + if key in lookup: + return lookup[key], True + + return None, False diff --git a/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/report.py b/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/report.py new file mode 100644 index 00000000..6b533cba --- /dev/null +++ b/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/report.py @@ -0,0 +1,107 @@ +from pathlib import Path +from datetime import date +import pandas as pd +from jinja2 import Environment, FileSystemLoader + + +def _cov_cell(a1, a2, exon): + """Display string for one exon across both alleles, e.g. '262x / 12x(!)'.""" + parts = [] + for a in (a1, a2): + depth = a.get(f'{exon}_depth') if a else None + if a is None or pd.isna(depth): + parts.append('—') + continue + mark = '(!)' if (a.get(f'{exon}_incomp') or 0) > 0 else '' + parts.append(f"{depth:.0f}x{mark}") + return ' / '.join(parts) + + +def _cov_detail(label, a): + """Per-allele breakdown line, or None if no coverage present.""" + if a is None or (pd.isna(a.get('exon2_depth')) and pd.isna(a.get('exon3_depth'))): + return None + def fmt(exon): + depth = a.get(f'{exon}_depth') + incomp = a.get(f'{exon}_incomp') + if pd.isna(depth): + return 'n/a' + return f"{depth:.1f}x incomp.{int(incomp) if pd.notna(incomp) else '?'}" + return f"{label}: exon2 {fmt('exon2')} | exon3 {fmt('exon3')}" + + +def render_report( + df: pd.DataFrame, + sample: str, + pgroup_file: str, + outdir: Path, + template_dir: Path, + hlahd_version: str = 'v1.7.1', + pileup: dict | None = None, +) -> Path: + """Render Jinja2 HTML report for one sample. Returns path to written file.""" + env = Environment(loader=FileSystemLoader(str(template_dir)), autoescape=True) + template = env.get_template('report.html.j2') + + loci_data = [] + for locus in ['A', 'B', 'C']: + locus_df = df[df['locus'] == locus] + if locus_df.empty: + continue + + def _get(pos, _ldf=locus_df): + rows = _ldf[_ldf['allele_position'] == pos] + return rows.iloc[0].to_dict() if not rows.empty else None + + a1 = _get('allele1') + a2 = _get('allele2') + a1p2 = _get('allele1_pair2') + a2p2 = _get('allele2_pair2') + + flags = [] + if a1 and a1['multiple_best_pairs']: + flags.append('MULTIPLE_BEST_PAIRS') + if a1 and a1['has_ambiguous_pair']: + flags.append('AMBIGUOUS_PAIR') + if (a1 and a1['est_mismatch']) or (a2 and a2['est_mismatch']): + flags.append('EST_MISMATCH') + + if (a1 and a1.get('incomplete_coverage')) or (a2 and a2.get('incomplete_coverage')): + flags.append('INCOMPLETE_COVERAGE') + + coverage_detail = [ + line for line in (_cov_detail('Allele 1', a1), _cov_detail('Allele 2', a2)) + if line + ] + + loci_data.append({ + 'locus': locus, + 'allele1': a1['allele'] if a1 else None, + 'p_group1': a1['p_group'] if a1 and pd.notna(a1['p_group']) else '—', + 'allele2': a2['allele'] if a2 else None, + 'p_group2': a2['p_group'] if a2 and pd.notna(a2['p_group']) else '—', + 'flags': flags, + 'ambiguous_alleles': a1['ambiguous_alleles'] if a1 else None, + 'est_mismatch_detail': ( + (a1['est_mismatch_detail'] if a1 and a1['est_mismatch'] else None) + or (a2['est_mismatch_detail'] if a2 and a2['est_mismatch'] else None) + ), + 'pair2_allele1': a1p2['allele'] if a1p2 else None, + 'pair2_allele2': a2p2['allele'] if a2p2 else None, + 'exon2_cov': _cov_cell(a1, a2, 'exon2'), + 'exon3_cov': _cov_cell(a1, a2, 'exon3'), + 'coverage_detail': coverage_detail, + }) + + html = template.render( + sample=sample, + loci=loci_data, + pgroup_file=pgroup_file, + run_date=date.today().isoformat(), + hlahd_version=hlahd_version, + pileup=pileup or {}, + ) + + out_file = Path(outdir) / f"{sample}_report.html" + out_file.write_text(html, encoding='utf-8') + return out_file diff --git a/modules/msk/annotate_hlahd/resources/usr/bin/templates/report.html.j2 b/modules/msk/annotate_hlahd/resources/usr/bin/templates/report.html.j2 new file mode 100644 index 00000000..f0c6da32 --- /dev/null +++ b/modules/msk/annotate_hlahd/resources/usr/bin/templates/report.html.j2 @@ -0,0 +1,135 @@ + + + + +HLA-HD Report — {{ sample }} + + + + +

HLA Class I Report

+

Sample: {{ sample }}  + Date: {{ run_date }}  + HLA-HD: {{ hlahd_version }}

+ +

HLA Class I Summary

+ + + + + + + + + + + + {% for locus in loci %} + + + + + + + + + + + {% endfor %} + +
LocusAllele 1P-group 1Allele 2P-group 2Exon 2 covExon 3 covFlags
HLA-{{ locus.locus }}{{ locus.allele1 or 'Not typed' }}{{ locus.p_group1 }}{{ locus.allele2 or 'Not typed' }}{{ locus.p_group2 }}{{ locus.exon2_cov }}{{ locus.exon3_cov }} + {% for flag in locus.flags %} + {% if flag == 'AMBIGUOUS_PAIR' %} + AMBIGUOUS PAIR + {% elif flag == 'MULTIPLE_BEST_PAIRS' %} + MULTIPLE BEST PAIRS + {% elif flag == 'EST_MISMATCH' %} + EST MISMATCH + {% elif flag == 'INCOMPLETE_COVERAGE' %} + INCOMPLETE COVERAGE + {% endif %} + {% endfor %} +
+ +{% set flagged = loci | selectattr('flags') | list %} +{% if flagged %} +

Flag Details

+{% for locus in flagged %} +
+ HLA-{{ locus.locus }}{% for f in locus.flags %} [{{ f }}]{% endfor %} + {% if 'AMBIGUOUS_PAIR' in locus.flags %} +

Ambiguous pair (not selected): {{ locus.ambiguous_alleles }}

+ {% endif %} + {% if 'MULTIPLE_BEST_PAIRS' in locus.flags and locus.pair2_allele1 %} +

Second best pair: {{ locus.pair2_allele1 }} / {{ locus.pair2_allele2 }}

+ {% endif %} + {% if 'EST_MISMATCH' in locus.flags %} +

Est mismatch detail: {{ locus.est_mismatch_detail }}

+ {% endif %} + {% if 'INCOMPLETE_COVERAGE' in locus.flags %} +

Coverage (mean depth, incomplete positions):

+
    + {% for line in locus.coverage_detail %} +
  • {{ line }}
    • + {% endfor %} +
+ {% endif %} +
+{% endfor %} +{% endif %} + +
+ Flag definitions:
+ AMBIGUOUS PAIR: HLA-HD identified an equivalent allele combination. The displayed alleles were selected based on population frequency.
+ MULTIPLE BEST PAIRS: HLA-HD could not disambiguate between equally-scored allele pairs.
+ EST MISMATCH: Unexpected discrepancy between the final result allele and the alleles in the .est.txt file. Review manually.
+
+ INCOMPLETE COVERAGE: One or more best-pair alleles had positions in exon 2 or exon 3 with zero read coverage (incomp > 0). A "(!)" marks the affected allele. Lower-confidence call — review manually. +
+ +{% if pileup %} +

Read Evidence at Discriminating / Mismatch Positions

+

Stacked read counts (A/C/G/T) at the positions that distinguish the called alleles. + Hatched (top) = allele 1's expected base; dotted = allele 2's; solid = any other base; + red columns are low coverage. At est-mismatch / ambiguous positions the dotted segment + shows the alternative allele's expected base rather than allele 2's.

+{% for locus in ['A', 'B', 'C'] if locus in pileup %} +
+

HLA-{{ locus }} + Open in IGV

+ HLA-{{ locus }} read pileup +

{{ pileup[locus].caption }}

+
+{% endfor %} +{% endif %} + +
+ HLA-HD {{ hlahd_version }} • P-group reference: {{ pgroup_file }} +
+ + diff --git a/modules/msk/annotate_hlahd/tests/main.nf.test b/modules/msk/annotate_hlahd/tests/main.nf.test new file mode 100644 index 00000000..9ef53ac0 --- /dev/null +++ b/modules/msk/annotate_hlahd/tests/main.nf.test @@ -0,0 +1,115 @@ +nextflow_process { + + name "Test Process ANNOTATE_HLAHD" + script "../main.nf" + process "ANNOTATE_HLAHD" + + tag "modules" + tag "modules_nfcore" + tag "modules_msk" + tag "annotate_hlahd" + + // Test 1: full run with real (synthetic) data. + // The module's result_dir input expects a single directory containing + // _final.result.txt and _{A,B,C}.est.txt; test_data.config + // registers those as individual files, so this test stages them into one + // directory before invoking the process. + test("annotate_hlahd - result_dir - annotated tsv and html") { + + when { + process { + """ + def resultDir = File.createTempDir() + new File(resultDir, "test_sample_final.result.txt") + .bytes = file(params.test_data_mskcc['annotate_hlahd']['final_result_txt'], checkIfExists: true).bytes + new File(resultDir, "test_sample_A.est.txt") + .bytes = file(params.test_data_mskcc['annotate_hlahd']['est_a_txt'], checkIfExists: true).bytes + new File(resultDir, "test_sample_B.est.txt") + .bytes = file(params.test_data_mskcc['annotate_hlahd']['est_b_txt'], checkIfExists: true).bytes + new File(resultDir, "test_sample_C.est.txt") + .bytes = file(params.test_data_mskcc['annotate_hlahd']['est_c_txt'], checkIfExists: true).bytes + + input[0] = [ + [ id:'test_sample' ], + file(resultDir) + ] + input[1] = file(params.test_data_mskcc['annotate_hlahd']['pgroup_file'], checkIfExists: true) + """ + } + } + + then { + assertAll( + { assert process.success }, + { assert snapshot( + process.out.tsv, + process.out.report, + process.out.versions + ).match() } + ) + } + } + + // Test 2: --skip_html should suppress the (optional) report output. + test("annotate_hlahd - skip_html - no report emitted") { + + config "./nextflow.config" + + when { + process { + """ + def resultDir = File.createTempDir() + new File(resultDir, "test_sample_final.result.txt") + .bytes = file(params.test_data_mskcc['annotate_hlahd']['final_result_txt'], checkIfExists: true).bytes + new File(resultDir, "test_sample_A.est.txt") + .bytes = file(params.test_data_mskcc['annotate_hlahd']['est_a_txt'], checkIfExists: true).bytes + new File(resultDir, "test_sample_B.est.txt") + .bytes = file(params.test_data_mskcc['annotate_hlahd']['est_b_txt'], checkIfExists: true).bytes + new File(resultDir, "test_sample_C.est.txt") + .bytes = file(params.test_data_mskcc['annotate_hlahd']['est_c_txt'], checkIfExists: true).bytes + + input[0] = [ + [ id:'test_sample' ], + file(resultDir) + ] + input[1] = file(params.test_data_mskcc['annotate_hlahd']['pgroup_file'], checkIfExists: true) + """ + } + } + + then { + assertAll( + { assert process.success }, + { assert process.out.report == [] }, + { assert snapshot(process.out.tsv, process.out.versions).match() } + ) + } + } + + // Test 3: stub run (no container needed) + test("annotate_hlahd - stub") { + + options "-stub" + + when { + process { + """ + input[0] = [ + [ id:'test_sample' ], + file('result_dir') + ] + input[1] = file('hla_nom_p.txt') + """ + } + } + + then { + assertAll( + { assert process.success }, + { assert path(process.out.tsv.get(0).get(1)).exists() }, + { assert path(process.out.report.get(0).get(1)).exists() }, + { assert snapshot(process.out.versions).match() } + ) + } + } +} diff --git a/modules/msk/annotate_hlahd/tests/main.nf.test.snap b/modules/msk/annotate_hlahd/tests/main.nf.test.snap new file mode 100644 index 00000000..d01962d6 --- /dev/null +++ b/modules/msk/annotate_hlahd/tests/main.nf.test.snap @@ -0,0 +1,62 @@ +{ + "annotate_hlahd - skip_html - no report emitted": { + "content": [ + [ + [ + { + "id": "test_sample" + }, + "test_sample_annotated.tsv:md5,a8636bc09b81afefdd01e388317b1978" + ] + ], + [ + "versions.yml:md5,7303efd9ae753a785695afaae7bca831" + ] + ], + "timestamp": "2026-08-26T16:36:59.315681", + "meta": { + "nf-test": "0.9.5", + "nextflow": "26.04.6" + } + }, + "annotate_hlahd - stub": { + "content": [ + [ + "versions.yml:md5,e1c7ca2bd9524d224402cef26a5fca25" + ] + ], + "timestamp": "2026-08-26T16:37:04.22217", + "meta": { + "nf-test": "0.9.5", + "nextflow": "26.04.6" + } + }, + "annotate_hlahd - result_dir - annotated tsv and html": { + "content": [ + [ + [ + { + "id": "test_sample" + }, + "test_sample_annotated.tsv:md5,a8636bc09b81afefdd01e388317b1978" + ] + ], + [ + [ + { + "id": "test_sample" + }, + "test_sample_report.html:md5,c78df8f69cf678d42b6aa6fad40d223a" + ] + ], + [ + "versions.yml:md5,7303efd9ae753a785695afaae7bca831" + ] + ], + "timestamp": "2026-08-26T16:36:51.100058", + "meta": { + "nf-test": "0.9.5", + "nextflow": "26.04.6" + } + } +} \ No newline at end of file diff --git a/modules/msk/annotate_hlahd/tests/nextflow.config b/modules/msk/annotate_hlahd/tests/nextflow.config new file mode 100644 index 00000000..52cfe423 --- /dev/null +++ b/modules/msk/annotate_hlahd/tests/nextflow.config @@ -0,0 +1,5 @@ +process { + withName: 'ANNOTATE_HLAHD' { + ext.args = '--skip_html' + } +} diff --git a/modules/msk/annotate_hlahd/tests/tags.yml b/modules/msk/annotate_hlahd/tests/tags.yml new file mode 100644 index 00000000..18559792 --- /dev/null +++ b/modules/msk/annotate_hlahd/tests/tags.yml @@ -0,0 +1,2 @@ +annotate_hlahd: + - "modules/msk/annotate_hlahd/**" diff --git a/tests/config/test_data.config b/tests/config/test_data.config index 347b83d6..f4786b2d 100644 --- a/tests/config/test_data.config +++ b/tests/config/test_data.config @@ -816,5 +816,13 @@ params { test_chr22_collapsed_grouped_bam = "${params.test_data_base_msk}/feature/fgbio_collectduplexseqmetrics/testdata/chr22_collapsed_grouped.bam" test_chr22_collapsed_grouped_bam_bai = "${params.test_data_base_msk}/feature/fgbio_collectduplexseqmetrics/testdata/chr22_collapsed_grouped.bam.bai" } + // NOTE: data lives on the feature/annotate_hlahd branch pending Review Team promotion to an official 'annotate_hlahd' branch. + 'annotate_hlahd' { + final_result_txt = "${params.test_data_base_msk}/feature/annotate_hlahd/annotate_hlahd/test_sample_final.result.txt" + est_a_txt = "${params.test_data_base_msk}/feature/annotate_hlahd/annotate_hlahd/test_sample_A.est.txt" + est_b_txt = "${params.test_data_base_msk}/feature/annotate_hlahd/annotate_hlahd/test_sample_B.est.txt" + est_c_txt = "${params.test_data_base_msk}/feature/annotate_hlahd/annotate_hlahd/test_sample_C.est.txt" + pgroup_file = "${params.test_data_base_msk}/feature/annotate_hlahd/annotate_hlahd/hla_nom_p_demo.txt" + } } } From 36a5d61877712b2cb18799f14a6a0181b8ffab5e Mon Sep 17 00:00:00 2001 From: John Orgera <65687576+johnoooh@users.noreply.github.com> Date: Wed, 26 Aug 2026 17:36:53 -0400 Subject: [PATCH 2/2] refactor: use mskcc/hla-annotate package instead of vendored copy The hlahd-tools:1.0.0 container now installs mskcc/hla-annotate (see companion containers PR) and exposes the 'annotate_hlahd' console script, so the module no longer needs to vendor a copy of the script under resources/usr/bin/. Re-verified: nf-test suite (real-data, --skip_html, stub) passes against the updated image -- identical output hashes to the vendored version. --- modules/msk/annotate_hlahd/main.nf | 2 +- modules/msk/annotate_hlahd/meta.yml | 4 +- .../resources/usr/bin/annotate_hlahd.py | 67 ------- .../usr/bin/hlahd_annotate/__init__.py | 0 .../usr/bin/hlahd_annotate/annotate.py | 107 ------------ .../usr/bin/hlahd_annotate/parse_est.py | 163 ------------------ .../usr/bin/hlahd_annotate/parse_final.py | 62 ------- .../usr/bin/hlahd_annotate/pgroup.py | 71 -------- .../usr/bin/hlahd_annotate/report.py | 107 ------------ .../usr/bin/templates/report.html.j2 | 135 --------------- .../annotate_hlahd/tests/main.nf.test.snap | 6 +- 11 files changed, 6 insertions(+), 718 deletions(-) delete mode 100755 modules/msk/annotate_hlahd/resources/usr/bin/annotate_hlahd.py delete mode 100644 modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/__init__.py delete mode 100644 modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/annotate.py delete mode 100644 modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/parse_est.py delete mode 100644 modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/parse_final.py delete mode 100644 modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/pgroup.py delete mode 100644 modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/report.py delete mode 100644 modules/msk/annotate_hlahd/resources/usr/bin/templates/report.html.j2 diff --git a/modules/msk/annotate_hlahd/main.nf b/modules/msk/annotate_hlahd/main.nf index b9abff28..1b2c2898 100644 --- a/modules/msk/annotate_hlahd/main.nf +++ b/modules/msk/annotate_hlahd/main.nf @@ -23,7 +23,7 @@ process ANNOTATE_HLAHD { def args = task.ext.args ?: '' prefix = task.ext.prefix ?: "${meta.id}" """ - annotate_hlahd.py \\ + annotate_hlahd \\ --result_dir ${result_dir} \\ --sample ${prefix} \\ --pgroup_file ${pgroup_file} \\ diff --git a/modules/msk/annotate_hlahd/meta.yml b/modules/msk/annotate_hlahd/meta.yml index b0f63aef..f0a4f2da 100644 --- a/modules/msk/annotate_hlahd/meta.yml +++ b/modules/msk/annotate_hlahd/meta.yml @@ -9,8 +9,8 @@ keywords: tools: - "annotate_hlahd": description: "Post-processes HLA-HD's per-sample final.result.txt/*.est.txt class I output, mapping each allele call to its IMGT P-group and flagging low-confidence calls." - homepage: "https://github.com/mskcc/HLA_HD_workflow" - documentation: "https://github.com/mskcc/HLA_HD_workflow" + homepage: "https://github.com/mskcc/hla-annotate" + documentation: "https://github.com/mskcc/hla-annotate" licence: - "MIT" identifier: "" diff --git a/modules/msk/annotate_hlahd/resources/usr/bin/annotate_hlahd.py b/modules/msk/annotate_hlahd/resources/usr/bin/annotate_hlahd.py deleted file mode 100755 index 9faa7e87..00000000 --- a/modules/msk/annotate_hlahd/resources/usr/bin/annotate_hlahd.py +++ /dev/null @@ -1,67 +0,0 @@ -#!/usr/bin/env python3 -""" -annotate_hlahd.py — Post-process HLA-HD class I output. - -Produces per-sample annotated TSV and self-contained HTML report. - -Usage: - annotate_hlahd.py \ - --result_dir results/SAMPLE1/ \ - --sample SAMPLE1 \ - --pgroup_file /path/to/hla_nom_p.txt \ - --outdir output/ - -Vendored from mskcc/HLA_HD_workflow's scripts/annotate_hlahd.py + scripts/hlahd_annotate/ -for the ANNOTATE_HLAHD module (see modules/msk/annotate_hlahd/resources/usr/bin/). -""" -import argparse -import sys -from pathlib import Path - -# hlahd_annotate/ is a sibling of this file (both live in resources/usr/bin/, -# which Nextflow adds to PATH for any process including this module). -sys.path.insert(0, str(Path(__file__).parent)) - -from hlahd_annotate.annotate import annotate_sample -from hlahd_annotate.report import render_report - -TEMPLATE_DIR = Path(__file__).parent / "templates" - - -def parse_args(): - p = argparse.ArgumentParser( - description='Annotate HLA-HD class I output with P-groups and quality flags.' - ) - p.add_argument('--result_dir', required=True, type=Path, - help='Directory containing _final.result.txt and _{A,B,C}.est.txt') - p.add_argument('--sample', required=True, - help='Sample ID (used as filename prefix for output files)') - p.add_argument('--pgroup_file', required=True, type=Path, - help='IMGT wmda/hla_nom_p.txt P-group reference table') - p.add_argument('--outdir', required=True, type=Path, - help='Output directory for TSV and HTML report (created if absent)') - p.add_argument('--hlahd_version', default='v1.7.1', - help='HLA-HD version string for the report footer (default: v1.7.1)') - p.add_argument('--skip_html', action='store_true', - help='Skip generating the HTML report; write only the annotated TSV') - return p.parse_args() - - -def main(): - args = parse_args() - args.outdir.mkdir(parents=True, exist_ok=True) - - df = annotate_sample(args.result_dir, args.sample, args.pgroup_file) - - tsv_path = args.outdir / f"{args.sample}_annotated.tsv" - df.to_csv(tsv_path, sep='\t', index=False) - print(f"TSV written: {tsv_path}") - - if not args.skip_html: - html_path = render_report(df, args.sample, str(args.pgroup_file), args.outdir, TEMPLATE_DIR, - hlahd_version=args.hlahd_version) - print(f"HTML written: {html_path}") - - -if __name__ == '__main__': - main() diff --git a/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/__init__.py b/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/__init__.py deleted file mode 100644 index e69de29b..00000000 diff --git a/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/annotate.py b/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/annotate.py deleted file mode 100644 index 2ec85844..00000000 --- a/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/annotate.py +++ /dev/null @@ -1,107 +0,0 @@ -import sys -from pathlib import Path -import pandas as pd - -from .pgroup import load_pgroup_table, lookup_pgroup -from .parse_final import parse_final_result, CLASS_I_LOCI -from .parse_est import parse_est, cross_check_allele - -# Map allele_position to (best_pairs index, coverage key within that pair) -_POSITION_TO_COVERAGE = { - 'allele1': (0, 'allele1_coverage'), - 'allele2': (0, 'allele2_coverage'), - 'allele1_pair2': (1, 'allele1_coverage'), - 'allele2_pair2': (1, 'allele2_coverage'), -} - - -def annotate_sample( - result_dir: Path, - sample: str, - pgroup_file: Path, -) -> pd.DataFrame: - """ - Combine parse_final, parse_est, and P-group lookup into the annotated DataFrame. - Emits a WARNING to stderr for each allele that fails est cross-check. - - Returns DataFrame with columns: - sample, locus, allele, allele_position, resolution, - p_group, p_group_found, multiple_best_pairs, has_ambiguous_pair, - ambiguous_alleles, est_mismatch, est_mismatch_detail, - exon2_depth, exon2_incomp, exon3_depth, exon3_incomp, incomplete_coverage - """ - result_dir = Path(result_dir) - final_file = result_dir / f"{sample}_final.result.txt" - - df = parse_final_result(final_file) - pgroup_lookup = load_pgroup_table(pgroup_file) - - # Load est data for each class I locus (returns empty-flag dict if file missing) - est_data: dict[str, dict] = {} - for locus in CLASS_I_LOCI: - est_file = result_dir / f"{sample}_{locus}.est.txt" - est_data[locus] = parse_est(est_file) - - records = [] - for _, row in df.iterrows(): - locus = row['locus'] - allele = row['allele'] - est = est_data.get(locus, {}) - is_pair2 = 'pair2' in row['allele_position'] - - # P-group lookup - p_group, p_group_found = lookup_pgroup(allele, pgroup_lookup) - - # Cross-check (primary pair alleles only; pair2 alleles not cross-checked) - est_mismatch = False - est_mismatch_detail = None - if allele and not is_pair2: - match, detail = cross_check_allele(allele, est.get('best_pair_alleles', [])) - if not match: - est_mismatch = True - est_mismatch_detail = detail - print( - f"WARNING: est_mismatch for sample={sample} locus={locus} " - f"allele_position={row['allele_position']}: {detail}", - file=sys.stderr, - ) - - # multiple_best_pairs: prefer est.txt (more precise) over final.result.txt count - multiple_best_pairs = est.get('multiple_best_pairs', row['multiple_best_pairs']) - - # Per-allele exon coverage from est best_pairs - coverage = {'exon2_depth': None, 'exon2_incomp': None, - 'exon3_depth': None, 'exon3_incomp': None} - best_pairs = est.get('best_pairs', []) - # Unrecognized allele_position falls through to all-None coverage (intended silent-null) - mapping = _POSITION_TO_COVERAGE.get(row['allele_position']) - if mapping: - idx, key = mapping - if idx < len(best_pairs): - coverage = dict(best_pairs[idx][key]) - incomplete_coverage = ( - (coverage['exon2_incomp'] or 0) > 0 - or (coverage['exon3_incomp'] or 0) > 0 - ) - - records.append({ - 'sample': sample, - 'locus': locus, - 'allele': allele, - 'allele_position': row['allele_position'], - 'resolution': row['resolution'], - 'p_group': p_group, - 'p_group_found': p_group_found, - 'multiple_best_pairs': multiple_best_pairs, - 'has_ambiguous_pair': est.get('has_ambiguous_pair', False), - 'ambiguous_alleles': est.get('ambiguous_alleles'), - 'est_mismatch': est_mismatch, - 'est_mismatch_detail': est_mismatch_detail, - 'exon2_depth': coverage['exon2_depth'], - 'exon2_incomp': coverage['exon2_incomp'], - 'exon3_depth': coverage['exon3_depth'], - 'exon3_incomp': coverage['exon3_incomp'], - 'incomplete_coverage': incomplete_coverage, - }) - - return pd.DataFrame(records) diff --git a/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/parse_est.py b/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/parse_est.py deleted file mode 100644 index 5dd81946..00000000 --- a/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/parse_est.py +++ /dev/null @@ -1,163 +0,0 @@ -from pathlib import Path - -_NULL_ALLELE = {'-', ''} - -_EXONS = ('exon2', 'exon3') - - -def _parse_coverage(field: str) -> dict: - """ - Parse a coverage field like 'exon2:262.426:comp.0,exon3:303.069:incomp.7' - into {exon2_depth, exon2_incomp, exon3_depth, exon3_incomp}. - comp.N -> incomp 0; incomp.N -> incomp N. Missing/null -> None values. - """ - cov = {} - for e in _EXONS: - cov[f'{e}_depth'] = None - cov[f'{e}_incomp'] = None - - if not field or field.strip() in _NULL_ALLELE: - return cov - - for part in field.split(','): - bits = part.split(':') - if len(bits) != 3: - continue - exon, depth, status = bits[0].strip(), bits[1].strip(), bits[2].strip() - if exon not in _EXONS: - continue - try: - cov[f'{exon}_depth'] = float(depth) - except ValueError: - cov[f'{exon}_depth'] = None - if status.startswith('incomp.'): - try: - cov[f'{exon}_incomp'] = int(status.split('.', 1)[1]) - except ValueError: - cov[f'{exon}_incomp'] = None - elif status.startswith('comp'): - cov[f'{exon}_incomp'] = 0 - - return cov - - -def _first_allele(field: str) -> str | None: - """First allele in a comma-separated est column, HLA- stripped; None if null.""" - field = field.strip() - if field in _NULL_ALLELE: - return None - first = field.split(',')[0].strip().removeprefix('HLA-') - return first or None - - -def parse_est(est_file: Path) -> dict: - """ - Parse *_{A,B,C}.est.txt. Used for flag extraction and cross-validation only. - Never used as the authoritative source for allele calls. - - Returns dict with: - multiple_best_pairs (bool): True if #Best allele pair count > 1 - has_ambiguous_pair (bool): True if #Other ambiguous pair section present - ambiguous_alleles (str | None): 'allele1 / allele2' from ambiguous section - best_pair_alleles (list[str]): all alleles from best-pair data rows, HLA- stripped - best_pairs (list[dict]): per best pair, with representative names - allele1/allele2 (str | None) and allele1_coverage/allele2_coverage dicts - """ - result = { - 'multiple_best_pairs': False, - 'has_ambiguous_pair': False, - 'ambiguous_alleles': None, - 'best_pair_alleles': [], - 'best_pairs': [], - } - - if not est_file.exists(): - return result - - with open(est_file) as f: - lines = [line.rstrip('\n') for line in f] - - in_ambiguous_section = False - best_pair_alleles = [] - best_pairs = [] - - for line in lines: - if not line: - continue - - if line.startswith('#Best allele pair'): - # Format: "#Best allele pair\t" - parts = line.split('\t') - count = int(parts[1]) if len(parts) > 1 else 1 - result['multiple_best_pairs'] = count > 1 - in_ambiguous_section = False - - elif line.startswith('#Other ambiguous pair'): - result['has_ambiguous_pair'] = True - in_ambiguous_section = True - - elif line.startswith('#'): - in_ambiguous_section = False - - else: - cols = line.split('\t') - if in_ambiguous_section: - # Ambiguous data lines: single allele per column (not lists) - if len(cols) >= 2: - a1 = cols[0].strip().removeprefix('HLA-') - a2 = cols[1].strip().removeprefix('HLA-') - result['ambiguous_alleles'] = f"{a1} / {a2}" - else: - # Best pair data lines: col0=allele1_comma_list, col1=allele2_comma_list - for col in cols[:2]: - col = col.strip() - if col in _NULL_ALLELE: - continue - for a in col.split(','): - a = a.strip().removeprefix('HLA-') - if a: - best_pair_alleles.append(a) - allele1_cov = _parse_coverage(cols[2]) if len(cols) > 2 else _parse_coverage('') - allele2_cov = _parse_coverage(cols[3]) if len(cols) > 3 else _parse_coverage('') - best_pairs.append({ - 'allele1': _first_allele(cols[0]) if len(cols) > 0 else None, - 'allele2': _first_allele(cols[1]) if len(cols) > 1 else None, - 'allele1_coverage': allele1_cov, - 'allele2_coverage': allele2_cov, - }) - - result['best_pair_alleles'] = best_pair_alleles - result['best_pairs'] = best_pairs - return result - - -def cross_check_allele( - allele: str | None, - best_pair_alleles: list[str], -) -> tuple[bool, str | None]: - """ - Check that an allele from final.result.txt appears in the est.txt best pair allele lists. - Matching is performed at the field depth of the reported allele (2 or 3 fields). - - Returns (match_found, detail_string). - If match_found is True, detail_string is None. - Hard warning to stderr is the caller's responsibility. - """ - if not allele or not best_pair_alleles: - return True, None - - a = allele.removeprefix('HLA-') - fields = a.split(':') - n_fields = len(fields) - - truncated = set() - for bp in best_pair_alleles: - bp_fields = bp.split(':') - truncated.add(':'.join(bp_fields[:n_fields])) - - if a in truncated: - return True, None - - shown = sorted(truncated)[:10] - detail = f"final={allele}; est_best_alleles={','.join(shown)}" - return False, detail diff --git a/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/parse_final.py b/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/parse_final.py deleted file mode 100644 index 6564a448..00000000 --- a/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/parse_final.py +++ /dev/null @@ -1,62 +0,0 @@ -from pathlib import Path -import warnings -import pandas as pd - -CLASS_I_LOCI = {'A', 'B', 'C'} -_NULL_VALUES = {'Not typed', '-', ''} -# Positions for alleles: primary pair then optional second pair -_PAIR_POSITIONS = ['allele1', 'allele2', 'allele1_pair2', 'allele2_pair2'] - - -def parse_final_result(result_file: Path) -> pd.DataFrame: - """ - Parse *_final.result.txt. Returns only class I loci (A, B, C). - - Returns DataFrame with columns: - locus, allele, allele_position, resolution, multiple_best_pairs - where: - - allele: str with HLA- prefix, or None/NaN if Not typed / - - - allele_position: 'allele1', 'allele2', 'allele1_pair2', 'allele2_pair2' - - resolution: number of colon-separated fields (0 if null) - - multiple_best_pairs: True if row has >2 allele columns - """ - rows = [] - with open(result_file) as f: - for line in f: - line = line.strip() - if not line: - continue - parts = line.split('\t') - locus = parts[0] - if locus not in CLASS_I_LOCI: - continue - - alleles = parts[1:] - multiple_best_pairs = len(alleles) > 2 - - if len(alleles) > len(_PAIR_POSITIONS): - warnings.warn( - f"Row for locus {locus!r} has {len(alleles)} allele columns " - f"(max supported: {len(_PAIR_POSITIONS)}); extra alleles ignored.", - UserWarning, - stacklevel=2, - ) - - for pos, allele in zip(_PAIR_POSITIONS, alleles): - allele_val = None if allele in _NULL_VALUES else allele - if allele_val is not None: - fields = allele_val.removeprefix('HLA-').split(':') - resolution = len(fields) - else: - resolution = 0 - - rows.append({ - 'locus': locus, - 'allele': allele_val, - 'allele_position': pos, - 'resolution': resolution, - 'multiple_best_pairs': multiple_best_pairs, - }) - - return pd.DataFrame(rows, columns=['locus', 'allele', 'allele_position', - 'resolution', 'multiple_best_pairs']) diff --git a/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/pgroup.py b/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/pgroup.py deleted file mode 100644 index a0079729..00000000 --- a/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/pgroup.py +++ /dev/null @@ -1,71 +0,0 @@ -from pathlib import Path - - -def load_pgroup_table(pgroup_file: Path) -> dict[str, str]: - """ - Parse IMGT wmda hla_nom_p.txt. - Returns lookup dict: allele at 2/3/4-field resolution -> P-group designation. - e.g., 'A*02:01' -> 'A*02:01P', 'A*02:01:01' -> 'A*02:01P' - """ - lookup: dict[str, str] = {} - with open(pgroup_file) as f: - for line in f: - line = line.strip() - if not line or line.startswith('#'): - continue - parts = line.split(';') - if len(parts) != 3: - continue - locus_prefix = parts[0] # e.g., 'A*' - alleles_str = parts[1] # e.g., '02:01:01:01/02:01:01:02' - pg_raw = parts[2].strip() - if not pg_raw: - continue # skip alleles not assigned to any P-group - p_group = locus_prefix + pg_raw # e.g., 'A*02:01P' - - for allele_fields in alleles_str.split('/'): - allele_fields = allele_fields.strip() - if not allele_fields: - continue - # Strip trailing IMGT expression suffixes (N=null, L=low, S=secreted, - # Q=questionable, C=aberrant cytoplasm, A=aberrant, G=null genomic) - # from the last colon-field. IMGT uses uppercase only. - fields = allele_fields.split(':') - fields[-1] = fields[-1].rstrip('NLSQCAG') - # MAX_HLA_FIELDS = 4; build keys at 2-, 3-, and 4-field resolution - for n in range(2, min(5, len(fields) + 1)): - key = locus_prefix + ':'.join(fields[:n]) - if key not in lookup: - lookup[key] = p_group - return lookup - - -def lookup_pgroup(allele: str | None, lookup: dict[str, str]) -> tuple[str | None, bool]: - """ - Look up P-group for a reported allele. - - Strips HLA- prefix before lookup. Tries match at reported field depth, - then falls back to shorter fields (minimum 2 fields). - - Sentinel values treated as "not found" (returns (None, False)): - - None - - 'Not typed' (HLA-HD value when locus has insufficient reads) - - '-' (HLA-HD value when only one allele is identified) - - Returns: - (p_group, found): tuple of the P-group string and a bool indicating success. - If not found, p_group is None and found is False. - """ - if not allele or allele in ('Not typed', '-'): - return None, False - - a = allele.removeprefix('HLA-') - fields = a.split(':') - - # Try from full depth down to 2 fields - for n in range(len(fields), 1, -1): - key = ':'.join(fields[:n]) - if key in lookup: - return lookup[key], True - - return None, False diff --git a/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/report.py b/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/report.py deleted file mode 100644 index 6b533cba..00000000 --- a/modules/msk/annotate_hlahd/resources/usr/bin/hlahd_annotate/report.py +++ /dev/null @@ -1,107 +0,0 @@ -from pathlib import Path -from datetime import date -import pandas as pd -from jinja2 import Environment, FileSystemLoader - - -def _cov_cell(a1, a2, exon): - """Display string for one exon across both alleles, e.g. '262x / 12x(!)'.""" - parts = [] - for a in (a1, a2): - depth = a.get(f'{exon}_depth') if a else None - if a is None or pd.isna(depth): - parts.append('—') - continue - mark = '(!)' if (a.get(f'{exon}_incomp') or 0) > 0 else '' - parts.append(f"{depth:.0f}x{mark}") - return ' / '.join(parts) - - -def _cov_detail(label, a): - """Per-allele breakdown line, or None if no coverage present.""" - if a is None or (pd.isna(a.get('exon2_depth')) and pd.isna(a.get('exon3_depth'))): - return None - def fmt(exon): - depth = a.get(f'{exon}_depth') - incomp = a.get(f'{exon}_incomp') - if pd.isna(depth): - return 'n/a' - return f"{depth:.1f}x incomp.{int(incomp) if pd.notna(incomp) else '?'}" - return f"{label}: exon2 {fmt('exon2')} | exon3 {fmt('exon3')}" - - -def render_report( - df: pd.DataFrame, - sample: str, - pgroup_file: str, - outdir: Path, - template_dir: Path, - hlahd_version: str = 'v1.7.1', - pileup: dict | None = None, -) -> Path: - """Render Jinja2 HTML report for one sample. Returns path to written file.""" - env = Environment(loader=FileSystemLoader(str(template_dir)), autoescape=True) - template = env.get_template('report.html.j2') - - loci_data = [] - for locus in ['A', 'B', 'C']: - locus_df = df[df['locus'] == locus] - if locus_df.empty: - continue - - def _get(pos, _ldf=locus_df): - rows = _ldf[_ldf['allele_position'] == pos] - return rows.iloc[0].to_dict() if not rows.empty else None - - a1 = _get('allele1') - a2 = _get('allele2') - a1p2 = _get('allele1_pair2') - a2p2 = _get('allele2_pair2') - - flags = [] - if a1 and a1['multiple_best_pairs']: - flags.append('MULTIPLE_BEST_PAIRS') - if a1 and a1['has_ambiguous_pair']: - flags.append('AMBIGUOUS_PAIR') - if (a1 and a1['est_mismatch']) or (a2 and a2['est_mismatch']): - flags.append('EST_MISMATCH') - - if (a1 and a1.get('incomplete_coverage')) or (a2 and a2.get('incomplete_coverage')): - flags.append('INCOMPLETE_COVERAGE') - - coverage_detail = [ - line for line in (_cov_detail('Allele 1', a1), _cov_detail('Allele 2', a2)) - if line - ] - - loci_data.append({ - 'locus': locus, - 'allele1': a1['allele'] if a1 else None, - 'p_group1': a1['p_group'] if a1 and pd.notna(a1['p_group']) else '—', - 'allele2': a2['allele'] if a2 else None, - 'p_group2': a2['p_group'] if a2 and pd.notna(a2['p_group']) else '—', - 'flags': flags, - 'ambiguous_alleles': a1['ambiguous_alleles'] if a1 else None, - 'est_mismatch_detail': ( - (a1['est_mismatch_detail'] if a1 and a1['est_mismatch'] else None) - or (a2['est_mismatch_detail'] if a2 and a2['est_mismatch'] else None) - ), - 'pair2_allele1': a1p2['allele'] if a1p2 else None, - 'pair2_allele2': a2p2['allele'] if a2p2 else None, - 'exon2_cov': _cov_cell(a1, a2, 'exon2'), - 'exon3_cov': _cov_cell(a1, a2, 'exon3'), - 'coverage_detail': coverage_detail, - }) - - html = template.render( - sample=sample, - loci=loci_data, - pgroup_file=pgroup_file, - run_date=date.today().isoformat(), - hlahd_version=hlahd_version, - pileup=pileup or {}, - ) - - out_file = Path(outdir) / f"{sample}_report.html" - out_file.write_text(html, encoding='utf-8') - return out_file diff --git a/modules/msk/annotate_hlahd/resources/usr/bin/templates/report.html.j2 b/modules/msk/annotate_hlahd/resources/usr/bin/templates/report.html.j2 deleted file mode 100644 index f0c6da32..00000000 --- a/modules/msk/annotate_hlahd/resources/usr/bin/templates/report.html.j2 +++ /dev/null @@ -1,135 +0,0 @@ - - - - -HLA-HD Report — {{ sample }} - - - - -

HLA Class I Report

-

Sample: {{ sample }}  - Date: {{ run_date }}  - HLA-HD: {{ hlahd_version }}

- -

HLA Class I Summary

- - - - - - - - - - - - {% for locus in loci %} - - - - - - - - - - - {% endfor %} - -
LocusAllele 1P-group 1Allele 2P-group 2Exon 2 covExon 3 covFlags
HLA-{{ locus.locus }}{{ locus.allele1 or 'Not typed' }}{{ locus.p_group1 }}{{ locus.allele2 or 'Not typed' }}{{ locus.p_group2 }}{{ locus.exon2_cov }}{{ locus.exon3_cov }} - {% for flag in locus.flags %} - {% if flag == 'AMBIGUOUS_PAIR' %} - AMBIGUOUS PAIR - {% elif flag == 'MULTIPLE_BEST_PAIRS' %} - MULTIPLE BEST PAIRS - {% elif flag == 'EST_MISMATCH' %} - EST MISMATCH - {% elif flag == 'INCOMPLETE_COVERAGE' %} - INCOMPLETE COVERAGE - {% endif %} - {% endfor %} -
- -{% set flagged = loci | selectattr('flags') | list %} -{% if flagged %} -

Flag Details

-{% for locus in flagged %} -
- HLA-{{ locus.locus }}{% for f in locus.flags %} [{{ f }}]{% endfor %} - {% if 'AMBIGUOUS_PAIR' in locus.flags %} -

Ambiguous pair (not selected): {{ locus.ambiguous_alleles }}

- {% endif %} - {% if 'MULTIPLE_BEST_PAIRS' in locus.flags and locus.pair2_allele1 %} -

Second best pair: {{ locus.pair2_allele1 }} / {{ locus.pair2_allele2 }}

- {% endif %} - {% if 'EST_MISMATCH' in locus.flags %} -

Est mismatch detail: {{ locus.est_mismatch_detail }}

- {% endif %} - {% if 'INCOMPLETE_COVERAGE' in locus.flags %} -

Coverage (mean depth, incomplete positions):

-
    - {% for line in locus.coverage_detail %} -
  • {{ line }}
    • - {% endfor %} -
- {% endif %} -
-{% endfor %} -{% endif %} - -
- Flag definitions:
- AMBIGUOUS PAIR: HLA-HD identified an equivalent allele combination. The displayed alleles were selected based on population frequency.
- MULTIPLE BEST PAIRS: HLA-HD could not disambiguate between equally-scored allele pairs.
- EST MISMATCH: Unexpected discrepancy between the final result allele and the alleles in the .est.txt file. Review manually.
-
- INCOMPLETE COVERAGE: One or more best-pair alleles had positions in exon 2 or exon 3 with zero read coverage (incomp > 0). A "(!)" marks the affected allele. Lower-confidence call — review manually. -
- -{% if pileup %} -

Read Evidence at Discriminating / Mismatch Positions

-

Stacked read counts (A/C/G/T) at the positions that distinguish the called alleles. - Hatched (top) = allele 1's expected base; dotted = allele 2's; solid = any other base; - red columns are low coverage. At est-mismatch / ambiguous positions the dotted segment - shows the alternative allele's expected base rather than allele 2's.

-{% for locus in ['A', 'B', 'C'] if locus in pileup %} -
-

HLA-{{ locus }} - Open in IGV

- HLA-{{ locus }} read pileup -

{{ pileup[locus].caption }}

-
-{% endfor %} -{% endif %} - - - - diff --git a/modules/msk/annotate_hlahd/tests/main.nf.test.snap b/modules/msk/annotate_hlahd/tests/main.nf.test.snap index d01962d6..e734df81 100644 --- a/modules/msk/annotate_hlahd/tests/main.nf.test.snap +++ b/modules/msk/annotate_hlahd/tests/main.nf.test.snap @@ -13,7 +13,7 @@ "versions.yml:md5,7303efd9ae753a785695afaae7bca831" ] ], - "timestamp": "2026-08-26T16:36:59.315681", + "timestamp": "2026-08-26T17:30:52.702765", "meta": { "nf-test": "0.9.5", "nextflow": "26.04.6" @@ -25,7 +25,7 @@ "versions.yml:md5,e1c7ca2bd9524d224402cef26a5fca25" ] ], - "timestamp": "2026-08-26T16:37:04.22217", + "timestamp": "2026-08-26T17:30:57.636612", "meta": { "nf-test": "0.9.5", "nextflow": "26.04.6" @@ -53,7 +53,7 @@ "versions.yml:md5,7303efd9ae753a785695afaae7bca831" ] ], - "timestamp": "2026-08-26T16:36:51.100058", + "timestamp": "2026-08-26T17:30:44.124301", "meta": { "nf-test": "0.9.5", "nextflow": "26.04.6"