Skip to content

Support AA abundance plot #25

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
mzueva wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Support AA abundance plot #25

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 8 additions & 0 deletions .changeset/tangy-snails-slide.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,8 @@
---
"@platforma-open/milaboratories.redefine-clonotypes.anarci-numbering": minor
"@platforma-open/milaboratories.redefine-clonotypes.workflow": minor
"@platforma-open/milaboratories.redefine-clonotypes.model": minor
"@platforma-open/milaboratories.redefine-clonotypes.ui": minor
---

Introduce per position AA export to support AA abundance plot
1 change: 1 addition & 0 deletions model/src/index.ts
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,7 @@ export type BlockArgs = {
numberingScheme?: 'imgt' | 'kabat' | 'chothia';
mem?: number;
cpu?: number;
exportCdr3AaPositions?: boolean;
};

export const model = BlockModel.create()
Expand Down
26 changes: 26 additions & 0 deletions software/anarci-numbering/src/main.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -280,6 +280,7 @@ def main() -> None:
p.add_argument("--cdr_mapping_kl", required=False, help="CDR annotation mapping JSON for light chain")
p.add_argument("--out_tsv", required=True, help="Output TSV path")
p.add_argument("--stats_tsv", required=False, help="Output numbering stats TSV path")
p.add_argument("--positions_tsv", required=False, help="Output per-position CDR3 amino acid TSV path")
args = p.parse_args()

keys, seqs, fields = read_input_tsv(args.input_tsv)
Expand Down Expand Up @@ -361,6 +362,31 @@ def main() -> None:
}
pl.DataFrame(stats_data).write_csv(args.stats_tsv, separator="\t")

# Write per-position CDR3 amino acid data
if args.positions_tsv:
chain_domain = {"H": "IGHeavy", "KL": "IGLight"}
pos_rows: List[List[str]] = []
unique_keys = list(dict.fromkeys(keys))
for key in unique_keys:
for chain in chains:
anarci_rows, pos_labels = anarci_by_chain[chain]
residues = anarci_rows.get(key) if anarci_rows and pos_labels else None
if residues is None:
continue
ranges = REGION_RANGES[args.scheme][chain]
cdr3_start, cdr3_end = ranges["CDR3"]
for pos_label, residue in zip(pos_labels, residues):
num = position_number(pos_label)
if num is None or num < cdr3_start or num > cdr3_end:
continue
residue = (residue or "").strip()
if residue in {"", "-", "."}:
continue
pos_rows.append([key, chain_domain[chain], pos_label, residue, "1"])
pl.DataFrame(
pos_rows, schema=["clonotypeKey", "chain", "position", "aminoacid", "count"], orient="row"
).write_csv(args.positions_tsv, separator="\t")
Comment on lines +366 to +388

@gemini-code-assist gemini-code-assist Bot Mar 27, 2026

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

high

The current implementation for generating the per-position CDR3 amino acid data builds a list of lists in memory, which can be inefficient and consume a large amount of memory for datasets with many clonotypes. Additionally, the count column is being written as a string "1", but the downstream workflow expects a numeric type (Long).

I suggest refactoring this section to use a more idiomatic polars approach. This will improve performance and memory efficiency by leveraging polars' lazy evaluation and optimized backend operations. The suggested change also corrects the data type of the count column.

    if args.positions_tsv:
        chain_domain = {"H": "IGHeavy", "KL": "IGLight"}
        all_pos_dfs = []
        for chain in chains:
            anarci_rows, pos_labels = anarci_by_chain[chain]
            if not anarci_rows or not pos_labels:
                continue

            ranges = REGION_RANGES[args.scheme][chain]
            cdr3_start, cdr3_end = ranges["CDR3"]

            df = pl.DataFrame(
                list(anarci_rows.items()),
                schema=["clonotypeKey", "residues"],
            )
            df = df.with_columns(pl.lit(pos_labels).alias("pos_labels"))
            df = df.explode(["residues", "pos_labels"])
            df = df.rename({"residues": "aminoacid", "pos_labels": "position"})

            df = df.with_columns(
                pl.col("position").str.extract(r"^(\d+)", 1).cast(pl.Int64).alias("pos_num")
            )
            df = df.filter(
                (pl.col("pos_num") >= cdr3_start)
                & (pl.col("pos_num") <= cdr3_end)
                & (~pl.col("aminoacid").str.strip().is_in(["", "-", "."]))
            )

            df = df.with_columns(
                pl.lit(chain_domain[chain]).alias("chain"),
                pl.lit(1, dtype=pl.Int64).alias("count"),
            )

            all_pos_dfs.append(df.select("clonotypeKey", "chain", "position", "aminoacid", "count"))

        if all_pos_dfs:
            final_df = pl.concat(all_pos_dfs)
            final_df.write_csv(args.positions_tsv, separator="\t")
        else:
            pl.DataFrame(
                schema={
                    "clonotypeKey": pl.Utf8,
                    "chain": pl.Utf8,
                    "position": pl.Utf8,
                    "aminoacid": pl.Utf8,
                    "count": pl.Int64,
                }
            ).write_csv(args.positions_tsv, separator="\t")



if __name__ == "__main__":
main()
16 changes: 15 additions & 1 deletion ui/src/pages/MainPage.vue
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
<script setup lang="ts">
import type { PlRef } from '@platforma-sdk/model';
import { PlAccordionSection, PlAlert, PlBlockPage, PlBtnGhost, PlDropdown, PlDropdownMulti, PlDropdownRef, PlLogView, PlMaskIcon24, PlNumberField, PlSectionSeparator, PlSlideModal, PlTabs } from '@platforma-sdk/ui-vue';
import { PlAccordionSection, PlAlert, PlBlockPage, PlBtnGhost, PlCheckbox, PlDropdown, PlDropdownMulti, PlDropdownRef, PlLogView, PlMaskIcon24, PlNumberField, PlSectionSeparator, PlSlideModal, PlTabs } from '@platforma-sdk/ui-vue';
import { computed, ref, watch, watchEffect } from 'vue';
import { useApp } from '../app';

Expand Down Expand Up @@ -33,6 +33,12 @@ watchEffect(() => {
}
});

watchEffect(() => {
if (app.model.args.numberingScheme === undefined) {
app.model.args.exportCdr3AaPositions = false;
}
});

// Stable fingerprint of chain options — reduces the option array to a primitive string
// so Vue's watch compares by value, not reference. This prevents false triggers from
// model re-evaluations that produce new array objects with identical contents.
Expand Down Expand Up @@ -170,6 +176,14 @@ function numberingWarningForChain(ns: { total: number; numbered: number } | unde
Apply IMGT, Kabat, or Chothia numbering. Requires datasets with VDJRegion or VDJRegionInFrame sequences or assembled on CDR3 (In this case, only the CDR3 region will be numbered). Transformed features are used for clonotype definition.
</template>
</PlDropdown>
<PlCheckbox
v-if="app.model.args.numberingScheme !== undefined"
:model-value="app.model.args.exportCdr3AaPositions ?? false"
:disabled="numberingDisabled"
@update:model-value="(v: boolean) => app.model.args.exportCdr3AaPositions = v"
>
Export CDR3 AA position table
</PlCheckbox>

<PlSectionSeparator>Resource Allocation</PlSectionSeparator>
<PlNumberField
Expand Down
14 changes: 12 additions & 2 deletions workflow/src/anarci-numbering.tpl.tengo
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ anarciSw := assets.importSoftware("@platforma-open/milaboratories.software-anarc
fastaSw := assets.importSoftware("@platforma-open/milaboratories.redefine-clonotypes.assembling-fasta:main")
numberingSw := assets.importSoftware("@platforma-open/milaboratories.redefine-clonotypes.anarci-numbering:main")

self.defineOutputs("numbered", "numberingStats", "anarciLog")
self.defineOutputs("numbered", "numberingStats", "anarciLog", "cdr3Positions")

self.body(func(inputs) {
inputAaTsv := inputs.inputAaTsv
Expand All @@ -16,6 +16,7 @@ self.body(func(inputs) {
bulkChain := inputs.bulkChain
cdrMappingH := inputs.cdrMappingH
cdrMappingKL := inputs.cdrMappingKL
exportCdr3AaPositions := inputs.exportCdr3AaPositions

// Resource configuration with defaults
fastaMem := "4GiB"
Expand Down Expand Up @@ -105,17 +106,26 @@ self.body(func(inputs) {
arg("--kl_csv").arg("anarci.csv_KL.csv")
}
}
if exportCdr3AaPositions {
numberingExec = numberingExec.
arg("--positions_tsv").arg("positions.tsv")
}
// Always write an empty placeholder so the output is never undefined
numberingExec = numberingExec.
writeFile("positions.tsv", "clonotypeKey\tchain\tposition\taminoacid\tcount\n")
numberingExec = numberingExec.
arg("--out_tsv").arg("numbered.tsv").
saveFile("numbered.tsv").
saveFileContent("numbering_stats.tsv").
saveFile("positions.tsv").
printErrStreamToStdout().
cache(24 * 60 * 60 * 1000).
run()

return {
numbered: numberingExec.getFile("numbered.tsv"),
numberingStats: numberingExec.getFileContent("numbering_stats.tsv"),
anarciLog: anarciBuilder.getStdoutStream()
anarciLog: anarciBuilder.getStdoutStream(),
cdr3Positions: numberingExec.getFile("positions.tsv")
}
})
104 changes: 101 additions & 3 deletions workflow/src/main.tpl.tengo
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -216,6 +216,8 @@ wf.body(func(args) {
}

pfBuilder := pframes.pFrameBuilder()
cdr3PosPfBuilder := pframes.pFrameBuilder()
hasCdr3Positions := false
outputs := {}
chainLabels := []

Expand Down Expand Up @@ -262,6 +264,7 @@ wf.body(func(args) {
numberingTsv := undefined
numberingStatsTsv := undefined
anarciLog := undefined
positionsTsv := undefined
numberingChains := []
cdr3Tsv := undefined
cdr3Chains := []
Expand Down Expand Up @@ -510,13 +513,17 @@ wf.body(func(args) {
clonotypeKeyAxisSpec: clonotypeKeyAxisSpec,
cdrMappingH: cdrMappingH,
cdrMappingKL: cdrMappingKL,
exportCdr3AaPositions: args.exportCdr3AaPositions == true,
mem: mem,
cpu: convCpu
})
numberingTsv = numberingPrep.output("numberingTsv", 24 * 60 * 60 * 1000)
numberingStatsTsv = numberingPrep.output("numberingStatsTsv", 24 * 60 * 60 * 1000)
anarciLog = numberingPrep.output("anarciLog")
cdr3Tsv = numberingPrep.output("cdr3Tsv")
if args.exportCdr3AaPositions == true {
positionsTsv = numberingPrep.output("positionsTsv", 24 * 60 * 60 * 1000)
}
}

anonymizationResult := render.create(anonymizationTpl, {
Expand Down Expand Up @@ -935,6 +942,16 @@ wf.body(func(args) {
df_abundances = df_abundances.withColumns(normExpressions...)
df_abundances.save("abundances.tsv")

// 5b. Expand per-position CDR3 amino acid data across samples
if positionsTsv != undefined {
positionsDf := ptWf.frame(positionsTsv, {xsvType: "tsv"})
positionsWithSample := df_abundances.
select(pt.col("sampleId"), pt.col("newClonotypeKey"), pt.col("clonotypeKey")).
join(positionsDf, {on: "clonotypeKey", how: "inner"}).
withoutColumns("clonotypeKey")
positionsWithSample.save("cdr3_positions.tsv")
}

// 6. Aggregate properties — pick representative from most-abundant old clonotype
df_properties := df_abundances.join(ptWf.frame(propertiesTsv, {xsvType: "tsv"}), {on: "clonotypeKey", how: "left"})
if numberingDf != undefined {
Expand Down Expand Up @@ -1121,13 +1138,89 @@ wf.body(func(args) {
partitionKeyLength: 0
}, { splitDataAndSpec: true, mem: mem, cpu: convCpu })

// Build per-position CDR3 amino acid p-frame for Graph Maker
cdr3PositionsPf := undefined
if positionsTsv != undefined {
cdr3PositionsTsvFile := ptResult.getFile("cdr3_positions.tsv")

// Deanonymize the positions TSV (sampleId axis was anonymized)
positionsDeanonimization := render.createEphemeral(
deanonimizationTpl,
{
mapping: mappingRef,
clusteringResult: cdr3PositionsTsvFile
}
)
cdr3PositionsTsvDeanon := positionsDeanonimization.output("deanonimizedTsv")

posSchemeLabel := numberingScheme == "imgt" ? "IMGT" : numberingScheme == "kabat" ? "Kabat" : "Chothia"
posLabel := posSchemeLabel + " CDR3 position"
posColumnLabel := posSchemeLabel + " CDR3 amino acid"
posDomain := {
"pl7.app/vdj/feature": "CDR3",
"pl7.app/vdj/numberingSchema": numberingScheme
}

cdr3PositionsPf = xsv.importFile(cdr3PositionsTsvDeanon, "tsv", {
axes: [
{ column: "sampleId", spec: sampleIdAxisSpec },
{ column: "newClonotypeKey", spec: newClonotypeKeySpec },
{ column: "chain", spec: {
name: "pl7.app/vdj/chain",
type: "String",
annotations: { "pl7.app/label": "Chain" }
}},
{ column: "position", spec: {
name: "pl7.app/vdj/numberingPosition",
type: "String",
domain: posDomain,
annotations: { "pl7.app/label": posLabel }
}}
],
columns: [
{
column: "aminoacid",
spec: {
name: "pl7.app/aminoacid",
valueType: "String",
domain: { "pl7.app/alphabet": "aminoacid" },
annotations: {
"pl7.app/label": posColumnLabel,
"pl7.app/discreteFilters": "true",
"pl7.app/table/visibility": "hidden"
}
}
},
{
column: "count",
spec: {
name: "pl7.app/count",
valueType: "Long",
annotations: {
"pl7.app/label": "AA Count by position",
"pl7.app/table/visibility": "hidden"
}
}
}
],
storageFormat: "Parquet",
partitionKeyLength: 1
}, { splitDataAndSpec: true, mem: "8GiB", cpu: 1 })
}

// Add this chain's PColumns to the combined PFrame
for k, v in abundancePf {
pfBuilder.add(anchor + "/" + k, trace.inject(v.spec), v.data)
}
for k, v in propertyPf {
pfBuilder.add(anchor + "/" + k, trace.inject(v.spec), v.data)
}
if cdr3PositionsPf != undefined {
hasCdr3Positions = true
for k, v in cdr3PositionsPf {
cdr3PosPfBuilder.add(anchor + "/" + k, trace.inject(v.spec), v.data)
}
}

// --- Collect per chain outputs ---
outputs["statsTsvContent_" + chainIdx] = statsTsvContent
Expand All @@ -1144,10 +1237,15 @@ wf.body(func(args) {
outputs.nChains = len(args.selectedChainRefs)
outputs.chainLabels = chainLabels

allExports := {
pf: pfBuilder.build()
}
if hasCdr3Positions {
allExports.cdr3AaPositions = cdr3PosPfBuilder.build()
}

return {
outputs: outputs,
exports: {
pf: pfBuilder.build()
}
exports: allExports
}
})
11 changes: 9 additions & 2 deletions workflow/src/numbering-prep.tpl.tengo
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ numberingTpl := assets.importTemplate(":anarci-numbering")
// Unused outputs get {} placeholder since defineOutputs requires all to be set.
self.defineOutputs(
"numberingTsv", "numberingStatsTsv", "anarciLog",
"cdr3Tsv"
"cdr3Tsv", "positionsTsv"
)

self.body(func(inputs) {
Expand All @@ -29,11 +29,13 @@ self.body(func(inputs) {
cdrMappingKL := inputs.cdrMappingKL
mem := inputs.mem
cpu := inputs.cpu
exportCdr3AaPositions := inputs.exportCdr3AaPositions

numberingTsv := undefined
numberingStatsTsv := undefined
anarciLog := undefined
cdr3Tsv := undefined
positionsTsv := undefined

// VDJRegion numbering via ANARCI
if useVdjRegionNumbering {
Expand Down Expand Up @@ -79,12 +81,16 @@ self.body(func(inputs) {
bulkChain: bulkChain,
cdrMappingH: cdrMappingH,
cdrMappingKL: cdrMappingKL,
exportCdr3AaPositions: exportCdr3AaPositions,
mem: mem,
cpu: cpu
})
numberingTsv = numbering.output("numbered", 24 * 60 * 60 * 1000)
numberingStatsTsv = numbering.output("numberingStats", 24 * 60 * 60 * 1000)
anarciLog = numbering.output("anarciLog")
if exportCdr3AaPositions {
positionsTsv = numbering.output("cdr3Positions", 24 * 60 * 60 * 1000)
}
}

// CDR3 TSV building (when VDJRegion is not available)
Expand Down Expand Up @@ -112,6 +118,7 @@ self.body(func(inputs) {
numberingTsv: numberingTsv != undefined ? numberingTsv : {},
numberingStatsTsv: numberingStatsTsv != undefined ? numberingStatsTsv : {},
anarciLog: anarciLog != undefined ? anarciLog : {},
cdr3Tsv: cdr3Tsv != undefined ? cdr3Tsv : {}
cdr3Tsv: cdr3Tsv != undefined ? cdr3Tsv : {},
positionsTsv: positionsTsv != undefined ? positionsTsv : {}
}
})
Loading