diff --git a/.changeset/disjoint-assembling-feature-fr3-imputation.md b/.changeset/disjoint-assembling-feature-fr3-imputation.md new file mode 100644 index 0000000..8007173 --- /dev/null +++ b/.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. diff --git a/model/src/index.ts b/model/src/index.ts index 1c787d2..5502262 100644 --- a/model/src/index.ts +++ b/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") { diff --git a/test/src/exportSpecs.test.ts b/test/src/exportSpecs.test.ts index be4a158..317641e 100644 --- a/test/src/exportSpecs.test.ts +++ b/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,6 +27,8 @@ 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); } @@ -30,6 +36,7 @@ function productiveFeature(assemblingFeature: string): string { // 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 { + const covered = new Set(); + 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"); + }); +}); diff --git a/test/src/wf.test.ts b/test/src/wf.test.ts index a3f96cb..2b511ea 100644 --- a/test/src/wf.test.ts +++ b/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; + + 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( + 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 }, diff --git a/ui/src/pages/SettingsPanel.vue b/ui/src/pages/SettingsPanel.vue index c223b86..812c173 100644 --- a/ui/src/pages/SettingsPanel.vue +++ b/ui/src/pages/SettingsPanel.vue @@ -241,8 +241,14 @@ const assemblingFeatureOptions = [ { value: "FR2:CDR3", label: "FR2:CDR3" }, { value: "CDR2:CDR3", label: "CDR2:CDR3" }, { value: "FR3:CDR3", label: "FR3:CDR3" }, + { value: "custom", label: "Custom (advanced)" }, ]; +// Real MiXCR feature values (everything except the "custom" sentinel used to reveal free text). +const presetAssemblingFeatures = assemblingFeatureOptions + .map((o) => o.value) + .filter((v) => v !== "custom"); + const assemblingFeature = computed({ get: () => app.model.args.assemblingFeature as AssemblingFeature, set: (value: AssemblingFeature) => { @@ -250,6 +256,42 @@ const assemblingFeature = computed({ }, }); +// "Custom (advanced)" mode: reveals a free-text field so an arbitrary MiXCR gene feature can be +// entered — including a disjoint one (e.g. "FR1Begin:FR3Begin(+40),FR3Begin(+46):FR4End") that +// brackets an uncovered window. Tracked as its own ref so an empty custom field doesn't snap the +// dropdown back to a preset; the watch turns it on when a project loads with a non-preset value. +const isCustomAssemblingFeature = ref(false); +watch( + () => app.model.args.assemblingFeature, + (v) => { + if (v !== undefined && v !== "" && !presetAssemblingFeatures.includes(v)) { + isCustomAssemblingFeature.value = true; + } + }, + { immediate: true }, +); + +const assemblingFeatureSelection = computed({ + get: () => { + if (isCustomAssemblingFeature.value) return "custom"; + const v = app.model.args.assemblingFeature; + if (v === undefined || v === "") return "VDJRegion"; + return presetAssemblingFeatures.includes(v) ? v : "custom"; + }, + set: (value: string) => { + if (value === "custom") { + isCustomAssemblingFeature.value = true; + // Don't carry a preset value into the custom text field. + if (presetAssemblingFeatures.includes(app.model.args.assemblingFeature ?? "")) { + app.model.args.assemblingFeature = ""; + } + } else { + isCustomAssemblingFeature.value = false; + app.model.args.assemblingFeature = value; + } + }, +}); + const imputeGermline = computed({ get: () => app.model.args.imputeGermline ?? false, set: (value: boolean) => { @@ -424,13 +466,36 @@ ATCGATCGATCG..." - + + + + + Impute non-covered parts from germline 3' order, used to classify a disjoint assembling feature. +featureOrder := ["FR1", "CDR1", "FR2", "CDR2", "FR3", "CDR3", "FR4"] + +featureIndex := func(r) { + for i, f in featureOrder { + if f == r { + return i + } + } + return -1 +} + +// Parse a MiXCR reference point ("FR3Begin", "FR3Begin(+40)", "FR4End", "FR3End(-20)") +// into { region, edge ("Begin"/"End"), offset }. +parseRefPoint := func(rp) { + parts := text.split(rp, "(") + base := parts[0] + offset := 0 + if len(parts) >= 2 { + inner := parts[1] + if text.has_suffix(inner, ")") { + inner = inner[0:len(inner) - 1] + } + neg := false + if len(inner) > 0 { + if inner[0] == '+' { + inner = inner[1:] + } else if inner[0] == '-' { + neg = true + inner = inner[1:] + } + } + o := int(inner) + if !is_undefined(o) { + offset = o + } + if neg { + offset = -offset + } + } + region := "" + edge := "" + if text.has_suffix(base, "Begin") { + edge = "Begin" + region = base[0:len(base) - 5] + } else if text.has_suffix(base, "End") { + edge = "End" + region = base[0:len(base) - 3] + } + return { + region: region, + edge: edge, + offset: offset + } +} + +// Given a disjoint assembling feature (comma-separated "begin:end" pieces with explicit +// reference points), return the set of whole region features fully covered by some piece. +// A region is fully covered only when a piece spans it edge-to-edge; a region a piece only +// reaches into via an offset (e.g. FR3 in "...:FR3Begin(+40)") is left uncovered so it gets +// germline-imputed on export. +coveredRegions := func(disjointStr) { + covered := {} + for p in text.split(disjointStr, ",") { + be := text.split(p, ":") + if len(be) != 2 { + continue + } + s := parseRefPoint(be[0]) + e := parseRefPoint(be[1]) + si := featureIndex(s.region) + ei := featureIndex(e.region) + if si == -1 || ei == -1 { + continue + } + first := si + if s.edge == "Begin" { + if s.offset > 0 { + first = si + 1 + } + } else { + first = si + 1 + } + last := ei + if e.edge == "End" { + if e.offset < 0 { + last = ei - 1 + } + } else { + last = ei - 1 + } + for i := first; i <= last; i++ { + if i >= 0 && i < len(featureOrder) { + covered[featureOrder[i]] = true + } + } + } + return covered +} + parseAssemblingFeature := func(assemblingFeature) { if assemblingFeature == "VDJRegion" || assemblingFeature == "CDR3" { return { @@ -67,6 +167,37 @@ parseAssemblingFeature := func(assemblingFeature) { } } + // Disjoint assembling feature: comma-separated pieces bracketing an uncovered mid-region + // window, e.g. "FR1Begin:FR3Begin(+40),FR3Begin(+46):FR4End" for the ~6 nt FR3 gap left by + // non-overlapping 2x150 reads on long-CDR3 clones. Each mate fully covers one piece, so the + // clone survives assembly; fully-covered whole regions are exported as-is (nonImputed), while + // the partially-covered gap region, any flanks, and the full VDJRegion are germline-imputed on + // export (MiXCR fills interior uncovered positions from the assigned V/J germline). + if text.contains(assemblingFeature, ",") { + covered := coveredRegions(assemblingFeature) + imputed := [] + nonImputed := [] + for f in featureOrder { + if !is_undefined(covered[f]) { + nonImputed = append(nonImputed, f) + } else { + imputed = append(imputed, f) + } + } + // The full VDJRegion always spans the gap, so it is never fully covered here. + imputed = append(imputed, "VDJRegion") + return { + imputed: imputed, + nonImputed: nonImputed, + // V/J mutation cores would span the imputed gap, so mutation columns are skipped for + // disjoint features (undefined core => the mutation loops no-op). + coreGeneFeatures: { + V: undefined, + J: undefined + } + } + } + be := text.split(assemblingFeature, ":") if len(be) != 2 { ll.panic("assemblingFeature must be in the format of 'begin:end', got " + assemblingFeature) @@ -159,6 +290,15 @@ calculateExportSpecs := func(presetSpecForBack, blockId) { if fstr == "VDJRegion" || fstr == "CDR3" { return fstr } + // Disjoint feature (comma-separated pieces): the composite covered sequence is NOT exported + // as a single column — MiXCR's GeneFeature.encode rewrites such a column header + // unpredictably (ref-point duals like FR3Begin->CDR2End, and whole-region ranges collapse to + // their name, e.g. {CDR3Begin:CDR3End}->CDR3). Instead we key on the individual covered + // regions and drive productivity off CDR3 (see the disjoint override + key branch below), so + // this return value is unused for disjoint features (productiveFeature is overridden). + if text.contains(fstr, ",") { + return fstr + } parts := text.split(fstr, ":") if len(parts) == 1 { return "{" + parts[0] + "Begin:" + parts[0] + "End}" @@ -206,6 +346,16 @@ calculateExportSpecs := func(presetSpecForBack, blockId) { anchorFeature = outputProductiveFeature } + // Disjoint feature: MiXCR can't export the composite covered sequence as a predictably-named + // single column, so drive productivity off CDR3 (always covered, predictable name) and make the + // full-length imputed VDJRegion the main sequence; the clonotype key is the covered regions + // (see the key branch below). + if text.contains(assemblingFeature, ",") { + productiveFeature = "CDR3" + outputProductiveFeature = "CDR3" + anchorFeature = "VDJRegion" + } + features := parsedFeature.nonImputed if imputeGermline { features = features + parsedFeature.imputed @@ -220,6 +370,20 @@ calculateExportSpecs := func(presetSpecForBack, blockId) { [ "-vGene" ], [ "-jGene" ] ] + } else if text.contains(assemblingFeature, ",") { + // Disjoint feature: key on the covered whole regions (predictable, verbatim column names, + // already exported in the feature loop) plus V and J. The composite covered sequence can't + // be a single predictably-named column; the covered regions together identify the clone + // (for synthetic fixed-framework libraries CDR1/CDR2 are independently diversified, so all + // covered regions are needed to distinguish members — CDR3 alone would collapse them). + for f in parsedFeature.nonImputed { + if f != "VDJRegion" { + clonotypeKeyColumns = append(clonotypeKeyColumns, "nSeq" + f) + clonotypeKeyArgs = append(clonotypeKeyArgs, [ "-nFeature", f ]) + } + } + clonotypeKeyColumns = append(clonotypeKeyColumns, "bestVGene", "bestJGene") + clonotypeKeyArgs = append(clonotypeKeyArgs, [ "-vGene" ], [ "-jGene" ]) } else { // VDJRegion is the assembling feature itself only when it's NOT in the imputed list // (e.g. VDJRegion or FR1:FR4 as the assembling feature) @@ -261,7 +425,9 @@ calculateExportSpecs := func(presetSpecForBack, blockId) { // For range features where VDJRegion is not the assembling feature, we need to export // the combined assembling feature sequence column explicitly (individual features are // exported in the loop below, but the combined feature like {CDR1Begin:CDR3End} is not) - needsAssemblingFeatureExport := assemblingFeature != "CDR3" && assemblingFeature != "VDJRegion" && !is_undefined(imputedFeaturesMap["VDJRegion"]) + // Disjoint features are excluded: their composite sequence has no predictably-named MiXCR + // column, so we key on covered regions instead and never export a combined column. + needsAssemblingFeatureExport := assemblingFeature != "CDR3" && assemblingFeature != "VDJRegion" && !text.contains(assemblingFeature, ",") && !is_undefined(imputedFeaturesMap["VDJRegion"]) if needsAssemblingFeatureExport { featureIdL := text.to_lower(formatId(assemblingFeature)) keyColName := "nSeq" + outputProductiveFeature diff --git a/workflow/src/mixcr-analyze.tpl.tengo b/workflow/src/mixcr-analyze.tpl.tengo index ac0b18c..9d7d6b4 100644 --- a/workflow/src/mixcr-analyze.tpl.tengo +++ b/workflow/src/mixcr-analyze.tpl.tengo @@ -69,10 +69,31 @@ self.body(func(inputs) { mixcrCmdBuilder.arg("generic-amplicon") } + // Format the assembling feature into MiXCR --assemble-clonotypes-by syntax. + // VDJRegion, CDR3 -> unchanged + // "X" / "X:Y" -> {XBegin:XEnd} / {XBegin:YEnd} + // disjoint (has ",") -> [{p1},{p2},...] with explicit reference points, e.g. + // "FR1Begin:FR3Begin(+40),FR3Begin(+46):FR4End" + // -> "[{FR1Begin:FR3Begin(+40)},{FR3Begin(+46):FR4End}]" + // NOTE: --assemble-clonotypes-by takes this LIST form. The export side + // (calculate-export-specs.lib.tengo) does NOT reference the composite feature at all — it keys + // the clonotype on the individual covered regions instead, because MiXCR renames a composite + // export column unpredictably. formatAssemblingFeature := func(fstr) { if fstr == "VDJRegion" || fstr == "CDR3" { return fstr } + if text.contains(fstr, ",") { + res := "[" + for i, p in text.split(fstr, ",") { + pbe := text.split(p, ":") + if i > 0 { + res += "," + } + res += "{" + pbe[0] + ":" + pbe[1] + "}" + } + return res + "]" + } parts := text.split(fstr, ":") if len(parts) == 1 { return "{" + parts[0] + "Begin:" + parts[0] + "End}"