Support disjoint assembling features for mid-region germline imputation (FR3 gap)

.changeset/disjoint-assembling-feature-fr3-imputation.md

@@ -0,0 +1,12 @@
+---
+'@platforma-open/milaboratories.mixcr-amplicon-alignment.workflow': minor
+'@platforma-open/milaboratories.mixcr-amplicon-alignment.ui': minor
+'@platforma-open/milaboratories.mixcr-amplicon-alignment.model': minor
+'@platforma-open/milaboratories.mixcr-amplicon-alignment': minor
+---
+
+Support disjoint assembling features for germline imputation across a mid-region (FR3) gap.
+
+When 2×150 reads don't overlap for long-CDR3 clones, a short uncovered window is left in FR3 and those clones are dropped at assembly. The assembling feature can now be a disjoint, comma-separated list of pieces (e.g. `FR1Begin:FR3Begin(+40),FR3Begin(+46):FR4End`) that brackets the gap: each mate fully covers one piece so the clone survives, fully-covered regions are exported as-is, and the skipped FR3 window plus the full VDJRegion are germline-imputed from the assigned V/J germline on export.
+
+The "Assembling feature" dropdown gains a "Custom (advanced)" option that reveals a free-text field for entering an arbitrary MiXCR gene feature, including such a disjoint one.

model/src/index.ts

@@ -260,6 +260,11 @@ export const platforma = BlockModel.create("Heavy")
   })
 
   .argsValid((ctx) => {
+    // A "Custom (advanced)" assembling feature left empty would reach parseAssemblingFeature and
+    // panic; block the run until it is filled.
+    if (ctx.args.assemblingFeature !== undefined && ctx.args.assemblingFeature.trim() === "") {
+      return false;
+    }
     const mode = ctx.uiState.referenceInputMode ?? "fastaSequence";
     const hasDataset = ctx.args.datasetRef !== undefined;
     if (mode === "libraryFile") {

test/src/exportSpecs.test.ts

@@ -14,6 +14,10 @@ import { test, expect, describe } from "vitest";
 // Mirrors formatAssemblingFeature in calculate-export-specs.lib.tengo
 function formatAssemblingFeature(fstr: string): string {
   if (fstr === "VDJRegion" || fstr === "CDR3") return fstr;
+  // Disjoint feature (comma-separated): the composite is NOT exported as a single column (MiXCR
+  // renames such headers unpredictably); productiveFeature is overridden to CDR3, so this value is
+  // unused for disjoint.
+  if (fstr.includes(",")) return fstr;
   const parts = fstr.split(":");
   if (parts.length === 1) return `{${parts[0]}Begin:${parts[0]}End}`;
   return `{${parts[0]}Begin:${parts[1]}End}`;
@@ -23,13 +27,16 @@ function formatAssemblingFeature(fstr: string): string {
 // MiXCR export args (-isProductive/-isOOF/-hasStops/-nFeature). FR1:FR4 is the full
 // VDJRegion, so it is normalized to "VDJRegion" to match MiXCR's column naming.
 function productiveFeature(assemblingFeature: string): string {
+  // Disjoint feature: productivity is driven off CDR3 (predictable, covered), not the composite.
+  if (assemblingFeature.includes(",")) return "CDR3";
   if (assemblingFeature === "FR1:FR4") return "VDJRegion";
   return formatAssemblingFeature(assemblingFeature);
 }
 
 // Mirrors outputProductiveFeature logic
 // MiXCR has named aliases for ranges ending at FR4; other ranges use {XBegin:YEnd}
 function outputProductiveFeature(assemblingFeature: string): string {
+  if (assemblingFeature.includes(",")) return "CDR3";
   const productive = formatAssemblingFeature(assemblingFeature);
   if (assemblingFeature !== "VDJRegion" && assemblingFeature !== "CDR3") {
     const parts = assemblingFeature.split(":");
@@ -42,6 +49,59 @@ function outputProductiveFeature(assemblingFeature: string): string {
   return productive;
 }
 
+// Region features in 5'->3' order, used to classify a disjoint assembling feature.
+const FEATURE_ORDER = ["FR1", "CDR1", "FR2", "CDR2", "FR3", "CDR3", "FR4"];
+
+// Mirrors parseRefPoint: "FR3Begin", "FR3Begin(+40)", "FR4End", "FR3End(-20)" -> parts
+function parseRefPoint(rp: string): { region: string; edge: string; offset: number } {
+  const [base, offsetPart] = rp.split("(");
+  let offset = 0;
+  if (offsetPart !== undefined) {
+    const n = parseInt(offsetPart.replace(")", ""), 10);
+    if (!Number.isNaN(n)) offset = n;
+  }
+  let region = "";
+  let edge = "";
+  if (base.endsWith("Begin")) {
+    edge = "Begin";
+    region = base.slice(0, -5);
+  } else if (base.endsWith("End")) {
+    edge = "End";
+    region = base.slice(0, -3);
+  }
+  return { region, edge, offset };
+}
+
+// Mirrors coveredRegions: whole regions fully spanned by some piece of a disjoint feature.
+function coveredRegions(disjointStr: string): Set<string> {
+  const covered = new Set<string>();
+  for (const p of disjointStr.split(",")) {
+    const be = p.split(":");
+    if (be.length !== 2) continue;
+    const s = parseRefPoint(be[0]);
+    const e = parseRefPoint(be[1]);
+    const si = FEATURE_ORDER.indexOf(s.region);
+    const ei = FEATURE_ORDER.indexOf(e.region);
+    if (si === -1 || ei === -1) continue;
+    let first = si;
+    if (s.edge === "Begin") {
+      if (s.offset > 0) first = si + 1;
+    } else {
+      first = si + 1;
+    }
+    let last = ei;
+    if (e.edge === "End") {
+      if (e.offset < 0) last = ei - 1;
+    } else {
+      last = ei - 1;
+    }
+    for (let i = first; i <= last; i++) {
+      if (i >= 0 && i < FEATURE_ORDER.length) covered.add(FEATURE_ORDER[i]);
+    }
+  }
+  return covered;
+}
+
 // Mirrors parseAssemblingFeature
 function parseAssemblingFeature(assemblingFeature: string) {
   if (assemblingFeature === "VDJRegion" || assemblingFeature === "CDR3") {
@@ -59,6 +119,20 @@ function parseAssemblingFeature(assemblingFeature: string) {
     };
   }
 
+  // Disjoint feature (comma-separated pieces): fully-covered regions are exported as-is,
+  // the partially-covered gap region + flanks + VDJRegion are germline-imputed.
+  if (assemblingFeature.includes(",")) {
+    const covered = coveredRegions(assemblingFeature);
+    const imputed: string[] = [];
+    const nonImputed: string[] = [];
+    for (const f of FEATURE_ORDER) {
+      if (covered.has(f)) nonImputed.push(f);
+      else imputed.push(f);
+    }
+    imputed.push("VDJRegion");
+    return { imputed, nonImputed };
+  }
+
   const features = ["FR1", "CDR1", "FR2", "CDR2", "FR3", "CDR3", "FR4"];
   const [begin, end] = assemblingFeature.split(":");
   const iBegin = features.indexOf(begin);
@@ -100,6 +174,13 @@ function computeClonotypeKeyAndExport(
 
   if (assemblingFeature === "CDR3") {
     clonotypeKeyColumns = ["nSeqCDR3", "bestVGene", "bestJGene"];
+  } else if (assemblingFeature.includes(",")) {
+    // Disjoint feature: key on the covered whole regions (predictable names) + V + J.
+    clonotypeKeyColumns = [
+      ...parsed.nonImputed.filter((f) => f !== "VDJRegion").map((f) => `nSeq${f}`),
+      "bestVGene",
+      "bestJGene",
+    ];
   } else {
     // VDJRegion is the assembling feature itself only when it's NOT in the imputed list
     const vdjIsAssemblingFeature = imputedFeaturesMap["VDJRegion"] === undefined;
@@ -116,7 +197,9 @@ function computeClonotypeKeyAndExport(
 
   const isRangeFeature = assemblingFeature !== "CDR3" && assemblingFeature !== "VDJRegion";
   const vdjIsImputed = imputedFeaturesMap["VDJRegion"] === true;
-  const needsAssemblingFeatureExport = isRangeFeature && vdjIsImputed;
+  // Disjoint features never export a combined column (composite has no predictable name).
+  const needsAssemblingFeatureExport =
+    isRangeFeature && vdjIsImputed && !assemblingFeature.includes(",");
 
   let assemblingFeatureColumn: string | undefined;
   if (needsAssemblingFeatureExport) {
@@ -323,3 +406,63 @@ describe("export-report flag column naming (productiveFeature)", () => {
     expect(productiveFeature("CDR1:CDR3")).toBe("{CDR1Begin:CDR3End}");
   });
 });
+
+// Disjoint assembling feature — recovers long-CDR3 clones whose 2x150 reads leave a ~6 nt
+// FR3 gap (Valerio "germline-imputation-fr3-gap"). The feature brackets the uncovered window;
+// each mate fully covers one piece, the clone survives assembly, and the gap is imputed from
+// germline on export. The two pieces reach into FR3 via offsets, so FR3 is only partially
+// covered and must be germline-imputed (never exported as a real, non-imputed column).
+describe("disjoint assembling feature (FR3-gap germline imputation)", () => {
+  // Valerio's real feature (excludes all of FR3). The composite covered sequence is NOT exported
+  // as a single column (MiXCR renames such headers unpredictably); instead the block keys on the
+  // covered regions and drives productivity off CDR3 — all predictable, verbatim column names.
+  const DISJOINT = "CDR1Begin:FR3Begin,CDR3Begin:FR4End";
+
+  test("covered regions are non-imputed; FR3 gap, flanks and VDJRegion are imputed", () => {
+    const r = parseAssemblingFeature(DISJOINT);
+    expect(r.nonImputed).toEqual(["CDR1", "FR2", "CDR2", "CDR3", "FR4"]);
+    expect(r.imputed).toContain("FR1");
+    expect(r.imputed).toContain("FR3");
+    expect(r.imputed).toContain("VDJRegion");
+    // FR3 is not covered — it must NOT be exported as a real (non-imputed) column, otherwise the
+    // "region_not_covered" placeholder would shadow the germline-imputed one.
+    expect(r.nonImputed).not.toContain("FR3");
+  });
+
+  test("clonotype key is the covered whole regions plus V and J (predictable names)", () => {
+    const r = computeClonotypeKeyAndExport(DISJOINT, true);
+    expect(r.clonotypeKeyColumns).toEqual([
+      "nSeqCDR1",
+      "nSeqFR2",
+      "nSeqCDR2",
+      "nSeqCDR3",
+      "nSeqFR4",
+      "bestVGene",
+      "bestJGene",
+    ]);
+    // No composite assembling-feature column is exported for disjoint features.
+    expect(r.needsAssemblingFeatureExport).toBe(false);
+  });
+
+  test("productivity is driven off CDR3 (predictable column name)", () => {
+    expect(`isProductive${outputProductiveFeature(DISJOINT)}`).toBe("isProductiveCDR3");
+  });
+
+  test("a wider gap leaves the same whole regions covered", () => {
+    // Offsets only change the imputed FR3 window width; whole-region coverage is unchanged.
+    const r = parseAssemblingFeature("FR1Begin:FR3Begin(+30),FR3Begin(+60):FR4End");
+    expect(r.nonImputed).toEqual(["FR1", "CDR1", "FR2", "CDR2", "CDR3", "FR4"]);
+    expect(r.imputed).toContain("FR3");
+  });
+
+  test("a 3-piece disjoint feature (two internal gaps) is classified by coverage", () => {
+    // Gaps in both FR2 and FR3: FR2 and FR3 are partial -> imputed; the rest stays covered.
+    const r = parseAssemblingFeature(
+      "FR1Begin:FR2Begin(+10),FR2Begin(+16):FR3Begin(+40),FR3Begin(+46):FR4End",
+    );
+    expect(r.nonImputed).toEqual(["FR1", "CDR1", "CDR2", "CDR3", "FR4"]);
+    expect(r.imputed).toContain("FR2");
+    expect(r.imputed).toContain("FR3");
+    expect(r.imputed).toContain("VDJRegion");
+  });
+});

test/src/wf.test.ts

@@ -264,6 +264,112 @@ blockTest(
   },
 );
 
+blockTest(
+  "disjoint FR3-gap feature with imputation",
+  { timeout: 300000 },
+  async ({ rawPrj: project, ml, helpers, expect }) => {
+    const sndBlockId = await project.addBlock("Samples & Data", samplesAndDataBlockSpec);
+    const alignBlockId = await project.addBlock("MiXCR Amplicon Alignment", myBlockSpec);
+
+    const sample1Id = uniquePlId();
+    const dataset1Id = uniquePlId();
+
+    const r1Handle = await helpers.getLocalFileHandle("./assets/s1_R1.fastq.gz");
+    const r2Handle = await helpers.getLocalFileHandle("./assets/s1_R2.fastq.gz");
+
+    await project.setBlockArgs(sndBlockId, {
+      metadata: [],
+      sampleIds: [sample1Id],
+      sampleLabelColumnLabel: "Sample Name",
+      sampleLabels: { [sample1Id]: "Sample 1" },
+      datasets: [
+        {
+          id: dataset1Id,
+          label: "Dataset 1",
+          content: {
+            type: "Fastq",
+            readIndices: ["R1", "R2"],
+            gzipped: true,
+            data: {
+              [sample1Id]: {
+                R1: r1Handle,
+                R2: r2Handle,
+              },
+            },
+          },
+        },
+      ],
+    } satisfies SamplesAndDataBlockArgs);
+    await project.runBlock(sndBlockId);
+
+    await helpers.awaitBlockDoneAndGetStableBlockState(sndBlockId, 8000);
+
+    // Wait for input options to propagate
+    const alignStableState1 = (await awaitStableState(
+      project.getBlockState(alignBlockId),
+      25000,
+    )) as InferBlockState<typeof platforma>;
+
+    const alignOutputs1 = wrapOutputs(alignStableState1.outputs);
+
+    // Configure the amplicon alignment block with a DISJOINT assembling feature: two pieces
+    // bracketing a mid-FR3 window. The window is excluded from the clonal sequence (so a
+    // long-CDR3 clone whose mates don't overlap there still assembles) and is germline-imputed
+    // on export. Run on s1 to verify the mechanism + MiXCR's disjoint column naming independently
+    // of a natural read gap (the disjoint feature excludes the window regardless of coverage).
+    const vGenesFasta = `>ref_heavy\n${referenceSequence}`;
+    const jGenesFasta = `>ref_heavy_j\n${referenceSequence.slice(-80)}`;
+
+    await project.setBlockArgs(alignBlockId, {
+      datasetRef: alignOutputs1.inputOptions[0].ref,
+      chains: "IGHeavy",
+      tagPattern: "",
+      vGenes: vGenesFasta,
+      jGenes: jGenesFasta,
+      assemblingFeature: "FR1Begin:FR3Begin(+30),FR3Begin(+36):FR4End",
+      imputeGermline: true,
+      cloneClusteringMode: "relaxed",
+    } satisfies BlockArgs);
+
+    await project.runBlock(alignBlockId);
+    const alignStableState3 = await helpers.awaitBlockDoneAndGetStableBlockState(
+      alignBlockId,
+      250000,
+    );
+    const outputs3 = wrapOutputs<BlockOutputs>(
+      alignStableState3.outputs as unknown as BlockOutputs,
+    );
+
+    // Reaching "done" with complete reports implicitly verifies that (a) the disjoint feature
+    // assembled clones across mates and (b) MiXCR's exported column names for the disjoint
+    // feature match what the workflow requests (a naming mismatch would error the export).
+    expect(outputs3.reports.isComplete).toEqual(true);
+
+    const reportEntries = outputs3.reports.data;
+    const alignJsonReportEntry = reportEntries.find(
+      (entry) => entry.key[1] === "align" && entry.key[2] === "json",
+    );
+    expect(alignJsonReportEntry).toBeDefined();
+
+    const alignReport = AlignReport.parse(
+      JSON.parse(
+        Buffer.from(
+          await ml.driverKit.blobDriver.getContent(
+            alignJsonReportEntry!.value!.handle as Parameters<
+              typeof ml.driverKit.blobDriver.getContent
+            >[0],
+          ),
+        ).toString("utf8"),
+      ),
+    );
+    expect(alignReport).toBeDefined();
+    expect(alignReport.totalReadsProcessed).greaterThan(0);
+
+    const qcEntry = outputs3.qc!.data[0];
+    expect(qcEntry).toBeDefined();
+  },
+);
+
 blockTest(
   "CDR1:CDR3 without imputation",
   { timeout: 300000 },