Add skeleton of workflow, stages, and jobs for dummy proband analysis

pyproject.toml

@@ -117,6 +117,7 @@ hail = ["hail"]
 # suppress the ARG002 "Unused method argument" warning in the stages.py file
 ## - we don't need generic cpg-flow arguments for every Stage, but need to fit the required method signature
 "src/cpg_seqr_loader/stages.py" = ["ARG002"]
+"src/cpg_seqr_loader/synthetic_proband_stages.py" = ["ARG002"]
 "src/cpg_seqr_loader/scripts/mt_to_es.py" = ["ANN202"]
 "src/cpg_seqr_loader/scripts/annotate_cohort.py" = ["E501"]
 

src/cpg_seqr_loader/jobs/CombineGvcfsIntoVds.py

@@ -127,15 +127,20 @@ def create_combiner_jobs(
         loguru.logger.info(f'Checking if VDS exists: {output_vds}: {output_vds.exists()}')
         return None
 
-    gvcf_add_arg = ''
-    if new_sg_gvcfs:
-        # write the gVCF paths into a temporary file
+    # If the synthetic-proband pre-workflow produced a pre-baked manifest (see the Stage 2 output
+    # `all_gvcf_paths_including_synthetic_gvcfs.txt`), use it verbatim and skip the normal
+    # metamist-driven SG discovery. Otherwise fall back to writing a manifest from new_sg_gvcfs.
+    if synthetic_gvcf_text_file := config.config_retrieve(['combiner', 'synthetic_gvcf_text_file'], None):
+        loguru.logger.info(f'Using synthetic gVCF text file: {synthetic_gvcf_text_file}')
+        localised_version = hail_batch.get_batch().read_input(str(synthetic_gvcf_text_file))
+        gvcf_add_arg = f'--gvcf_add_file {localised_version}'
+    elif new_sg_gvcfs:
         gvcf_path_file = temp_dir / 'gvcfs_to_combine.txt'
-        with gvcf_path_file.open('w') as write_handle:
-            for gvcf_path in new_sg_gvcfs:
-                write_handle.write(f'{gvcf_path!s}\n')
+        utils.write_gvcf_manifest([str(p) for p in new_sg_gvcfs], str(gvcf_path_file))
         localised_version = hail_batch.get_batch().read_input(str(gvcf_path_file))
         gvcf_add_arg = f'--gvcf_add_file {localised_version}'
+    else:
+        gvcf_add_arg = ''
 
     job = hail_batch.get_batch().new_bash_job(
         'CombineGvcfsIntoVds',

src/cpg_seqr_loader/jobs/GenerateSyntheticProbandCombinerInputs.py

@@ -0,0 +1,59 @@
+"""
+Driver-side writer for the synthetic proband combiner-inputs stage.
+
+The single public entry point `write_combiner_inputs` produces both of the stage's output
+artifacts from the driver, no Batch jobs required:
+
+  - Pedigree PED file: 6-column TSV with three rows per qualifying duo family (mother, father,
+    synthetic proband), plus a row for every real SG in the multicohort.
+
+  - gVCF manifest: newline-separated list of every gVCF that will go into the synthetic-trio
+    combiner run. It includes every real SG in the multicohort that has a gVCF (qualifying or
+    not, so we don't silently drop samples from the seqr load), plus every synthetic gVCF
+    from Stage 1's outputs.
+
+Both files are tiny text blobs, so the previous approach of spinning up a Batch job per file
+(only to run `cat > $output <<HEREDOC`) added VM + container startup cost and a heredoc-quoting
+hazard for no benefit. Writing directly to `gs://` via `cpg_utils.Path.open('w')` is the same
+pattern used elsewhere in the repo (see jobs/SubsetMtToDatasetWithHail.py).
+
+The write is skipped in dry-run mode - unlike Batch jobs, which cpg_flow describes rather than
+executes under `workflow.dry_run`, a driver-side `.open('w')` executes unconditionally, so the
+callers must gate on that config themselves.
+"""
+
+from cpg_flow import targets
+from cpg_utils import Path, config
+
+from cpg_seqr_loader.utils import (
+    SyntheticProbandFamily,
+    build_gvcf_manifest_content,
+    build_synthetic_pedigree_content,
+)
+
+
+def write_combiner_inputs(
+    families: list[SyntheticProbandFamily],
+    multicohort: targets.MultiCohort,
+    synthetic_gvcf_paths: dict[str, Path],
+    output_ped: Path,
+    output_manifest: Path,
+) -> None:
+    """Write the PED and gVCF manifest for the synthetic-trio combiner run.
+
+    No Batch jobs are queued - both files are written directly from the driver. The stage's
+    `queue_jobs` still returns an empty job list so cpg_flow's REUSE check keys off the output
+    paths existing on disk, exactly as it would with jobs.
+    """
+    if config.config_retrieve(['workflow', 'dry_run'], False):
+        return
+
+    pedigree_content = build_synthetic_pedigree_content(multicohort, families)
+    with output_ped.open('w') as f:
+        f.write(pedigree_content)
+
+    real_paths = [str(sg.gvcf) for sg in multicohort.get_sequencing_groups() if sg.gvcf is not None]
+    synthetic_paths = [str(synthetic_gvcf_paths[family.family_id]) for family in families]
+    manifest_content = build_gvcf_manifest_content(real_paths + synthetic_paths)
+    with output_manifest.open('w') as f:
+        f.write(manifest_content)

src/cpg_seqr_loader/jobs/GenerateSyntheticProbandGvcfs.py

@@ -0,0 +1,186 @@
+"""
+Job factory for the synthetic proband gVCF stage.
+
+The single public entry point `create_synthetic_gvcf_jobs` builds every Batch job the stage
+needs and returns them in queue order. Internally it uses two private per-artifact factories:
+
+  - _create_synthesis_jobs: one BashJob per duo family that invokes
+    scripts/create_synthetic_proband_gvcf.py with the family's parental gVCFs, writing the
+    output gVCF to the durable per-family path supplied by the stage.
+
+  - _create_analysis_registration_jobs: one BashJob per duo family that invokes
+    scripts/register_synthetic_gvcf_analysis.py to record (or refresh) the synthetic gVCF as a
+    metamist Analysis of type SYNTHETIC_GVCF_ANALYSIS_TYPE. Each registration job depends_on
+    the matching synthesis job so registration only fires once the gVCF actually exists.
+
+The synthesis script accepts gs:// paths directly (it localises inputs itself), so we don't
+need Batch's read_input_group here.
+"""
+
+from typing import TYPE_CHECKING
+
+from cpg_utils import Path, config, hail_batch, to_path
+
+from cpg_seqr_loader.utils import SyntheticProbandFamily
+
+if TYPE_CHECKING:
+    from hailtop.batch.job import BashJob
+
+
+def _create_synthesis_jobs(
+    families: list[SyntheticProbandFamily],
+    output_paths: dict[str, Path],
+    job_attrs: dict,
+) -> dict[str, 'BashJob']:
+    """Build one Batch synthesis job per qualifying duo family.
+
+    Args:
+        families: enumerated by utils.get_families_for_synthetic_probands.
+        output_paths: keyed by family_id; the durable output gVCF path per family.
+        job_attrs: base job attributes from the calling stage (dataset, workflow, etc.).
+
+    Returns:
+        Dict keyed by family_id -> BashJob, so downstream registration jobs can depend on the
+        matching synthesis job. Empty dict if no families qualify.
+    """
+    jobs = {}
+    for family in families:
+        output = output_paths[family.family_id]
+
+        # Per-family skip: cpg_flow's stage-level REUSE check is bypassed here (e.g. when the
+        # registration marker for a sibling family is missing and forces the whole stage to
+        # QUEUE). Avoid re-running expensive synthesis when this family's gVCF already exists.
+        if to_path(output).exists():
+            continue
+
+        job = hail_batch.get_batch().new_bash_job(
+            f'SynthesizeProband_{family.external_family_id}',
+            attributes=job_attrs | {'tool': 'pysam', 'family_id': family.external_family_id},
+        )
+        job.image(config.config_retrieve(['workflow', 'driver_image']))
+
+        # Sized for a WGS trio synthesis: 2 parental gVCFs (~20 GB each) + one output (~20 GB)
+        # + slack. Overridable via config so exome / smaller-genome runs can shrink this.
+        job.storage(config.config_retrieve(['synthetic_proband', 'storage'], '80G'))
+        job.memory(config.config_retrieve(['synthetic_proband', 'memory'], 'standard'))
+
+        job.command(
+            f"""
+            python -m cpg_seqr_loader.scripts.create_synthetic_proband_gvcf \\
+                --mother_gvcf {family.mother_sg.gvcf!s} \\
+                --father_gvcf {family.father_sg.gvcf!s} \\
+                --sample_name {family.synthetic_sample_name} \\
+                --output {output!s}
+            """,
+        )
+
+        jobs[family.family_id] = job
+
+    return jobs
+
+
+def _create_analysis_registration_jobs(
+    families: list[SyntheticProbandFamily],
+    gvcf_paths: dict[str, Path],
+    marker_paths: dict[str, Path],
+    synthesis_jobs: dict[str, 'BashJob'],
+    project: str,
+    job_attrs: dict,
+) -> list['BashJob']:
+    """Build one Batch job per family that registers (or refreshes) the metamist Analysis.
+
+    Each returned job depends on the matching synthesis job so registration only fires once the
+    synthetic gVCF actually exists. If the framework has skipped a synthesis job because the
+    output already exists, the corresponding key won't be present in synthesis_jobs; in that
+    case the registration job simply has no upstream dependency and runs immediately (the
+    script itself is idempotent).
+
+    Per-family skip: if this family's registration marker already exists, we don't queue a
+    registration job for it. Same rationale as the skip in create_synthesis_jobs.
+
+    Args:
+        families: same list passed to create_synthesis_jobs.
+        gvcf_paths: per-family output gVCF paths.
+        marker_paths: per-family registration sentinel paths. The script writes to this on
+            successful registration; its presence is what cpg_flow uses to decide REUSE.
+        synthesis_jobs: return value of create_synthesis_jobs (family_id -> BashJob).
+        project: metamist project name (the dataset the parental SGs live in).
+        job_attrs: base job attributes from the calling stage.
+    """
+    jobs = []
+    for family in families:
+        marker = marker_paths[family.family_id]
+        if to_path(marker).exists():
+            continue
+
+        gvcf = gvcf_paths[family.family_id]
+
+        job = hail_batch.get_batch().new_bash_job(
+            f'RegisterSyntheticProband_{family.external_family_id}',
+            attributes=job_attrs | {'tool': 'metamist', 'family_id': family.external_family_id},
+        )
+        job.image(config.config_retrieve(['workflow', 'driver_image']))
+
+        upstream = synthesis_jobs.get(family.family_id)
+        if upstream is not None:
+            job.depends_on(upstream)
+
+        job.command(
+            f"""
+            python -m cpg_seqr_loader.scripts.register_synthetic_gvcf_analysis \\
+                --project {project} \\
+                --gvcf_path {gvcf!s} \\
+                --sample_name {family.synthetic_sample_name} \\
+                --family_id {family.external_family_id} \\
+                --mother_sg_id {family.mother_sg.id} \\
+                --father_sg_id {family.father_sg.id} \\
+                --mother_source_gvcf {family.mother_sg.gvcf!s} \\
+                --father_source_gvcf {family.father_sg.gvcf!s} \\
+                --marker_path {marker!s}
+            """,
+        )
+
+        jobs.append(job)
+
+    return jobs
+
+
+def create_synthetic_gvcf_jobs(
+    families: list[SyntheticProbandFamily],
+    gvcf_paths: dict[str, Path],
+    marker_paths: dict[str, Path],
+    job_attrs: dict,
+) -> list['BashJob']:
+    """Build every Batch job the synthetic-gVCF stage needs, in the order they'll be queued.
+
+    Returns the synthesis jobs (one per family) followed by the registration jobs (one per
+    family). Registration jobs depend on their matching synthesis job so they only fire once
+    the gVCF exists.
+
+    Callers (the Stage class) should treat this as the single entry point for the stage - the
+    per-artifact factories are private to this module and shouldn't be invoked directly.
+    """
+    synthesis_jobs = _create_synthesis_jobs(
+        families=families,
+        output_paths=gvcf_paths,
+        job_attrs=job_attrs,
+    )
+
+    # Register per family: each Analysis lives in its parents' metamist project (which is the
+    # dataset owning the mother SG), so we can't batch across families that live in different
+    # datasets. cpg-flow's get_metamist().create_analysis handles the access-level suffix
+    # internally.
+    registration_jobs = []
+    for family in families:
+        registration_jobs.extend(
+            _create_analysis_registration_jobs(
+                families=[family],
+                gvcf_paths=gvcf_paths,
+                marker_paths=marker_paths,
+                synthesis_jobs=synthesis_jobs,
+                project=family.mother_sg.dataset.name,
+                job_attrs=job_attrs,
+            ),
+        )
+
+    return list(synthesis_jobs.values()) + registration_jobs