Leukemia (ALL vs AML) marker-gene classification — 100% linaldb, embedded

This notebook uses the real linaldb package installed from PyPI in linal-hub/.venv (same embedded PyO3 bindings as 01_linaldb_quickstart.ipynb), against a real, non-trivial dataset: the classic Golub et al. (1999) leukemia gene-expression microarray dataset (notebooks/data/examples/cancer_cell_detection/) — 7,129 genes measured across 72 patients, split into a 38-patient training set and a 34-patient independent test set, with each patient labeled ALL (Acute Lymphoblastic Leukemia) or AML (Acute Myeloid Leukemia).

Goal: push linaldb's own DSL as far as it will go for a real data-science workflow — marker-gene discovery, per-class centroids, similarity-based classification, accuracy scoring — instead of reaching for pandas/numpy/scikit-learn. The rule followed throughout:

• Python's stdlib csv module reads the raw files and reshapes rows into DSL literals — pure I/O glue, no statistics.
• Every actual computation (correlation, centroids, similarity, classification, accuracy) runs as linaldb DSL via db.execute(...), read back through ExecuteResult.rows / TensorResult.data — never .to_pandas() / .to_numpy() / .query().
• matplotlib appears exactly once, at the very end, to render two charts from already-computed linaldb results — nothing upstream of that touches it.

Along the way, building this notebook found (and got fixed) two real bugs in the engine itself — see Section 0 for the sanity check, and the closing section for the full story.

import csv
import time

import linaldb

print("linaldb", linaldb.__version__)

db = linaldb.Db(data_dir="./data")
db.active_db(), db.data_dir()

linaldb 0.1.9

('default', './data')

0. Sanity check: CORRELATE and centering a vector by its own mean

Earlier drafts of this notebook found two real bugs while exercising the engine this hard: CORRELATE a WITH b is documented (docs/DSL_REFERENCE.md) as Pearson correlation but was actually computing a raw, unnormalized dot product, and SUM/MEAN/STDEV returned the wrong tensor shape internally, which silently corrupted vector - MEAN(vector) (and similar ops) past the first element. Both are now fixed — engine v0.1.79, PyPI linaldb 0.1.2, which requirements.txt now requires. Full story in the closing section.

The cell below is a quick sanity check confirming both now behave correctly against a hand-computed ground truth, before trusting either at the scale of 7,129 genes in Section 3.

def mean(v):
    return sum(v) / len(v)


def center(v):
    m = mean(v)
    return [x - m for x in v]


def pearson_by_hand(x, y):
    xc, yc = center(x), center(y)
    num = sum(a * b for a, b in zip(xc, yc))
    denx = sum(a * a for a in xc) ** 0.5
    deny = sum(b * b for b in yc) ** 0.5
    return num / (denx * deny)


a = [1.0, 2.0, 3.0, 4.0]
b = [2.0, 4.0, 5.0, 9.0]
print("hand-computed Pearson correlation:", pearson_by_hand(a, b))

hand-computed Pearson correlation: 0.9647638212377322

db.execute(f"VECTOR proba = {a}")
db.execute(f"VECTOR probb = {b}")

db.execute("LET corr = CORRELATE proba WITH probb")
corr = db.execute("SHOW corr").data[0]
print("CORRELATE proba WITH probb ->", corr)
print("matches hand-computed Pearson correlation:", abs(corr - pearson_by_hand(a, b)) < 1e-5)

db.execute("LET mean_a = MEAN proba")
db.execute("LET proba_centered = proba - mean_a")
proba_centered = db.execute("SHOW proba_centered").data
print("proba - MEAN(proba) ->", proba_centered)
print("every element centered correctly:",
      all(abs(x - (v - mean(a))) < 1e-5 for x, v in zip(proba_centered, a)))

CORRELATE proba WITH probb -> 0.9647638201713562
matches hand-computed Pearson correlation: True
proba - MEAN(proba) -> [-1.5, -0.5, 0.5, 1.5]
every element centered correctly: True

Confirmed: CORRELATE now returns the true Pearson correlation directly, and vector - MEAN(vector) now centers every element (not just the first). Section 3 below calls CORRELATE natively on each gene — no Python-side mean-centering workaround needed anymore.

1. Load & reshape the raw CSVs (stdlib only, no computation)

The source files are gene-major: each row is one gene, each pair of columns is one patient's (expression value, present/absent/marginal call). We drop the call columns (a detection QC flag we don't need) and keep the numeric expression values. Patient IDs in the CSV header are not in numeric order (e.g. patients 34-38 appear before 28-33 in the training file), so we read them off the header rather than assuming a sequence. This section only re-indexes numbers already sitting in the files on disk — nothing here is a statistical computation.

DATA_DIR = "data/examples/cancer_cell_detection"


def load_labels(path):
    labels = {}
    with open(path, newline="") as f:
        reader = csv.reader(f)
        next(reader)  # header: patient,cancer
        for row in reader:
            if row:
                labels[int(row[0])] = row[1]
    return labels


def load_expression(path):
    """Returns (patient_ids_in_column_order, [(gene_id, description, [expr per patient]), ...])."""
    with open(path, newline="") as f:
        reader = csv.reader(f)
        header = next(reader)
        patient_cols = [(int(header[i]), i) for i in range(2, len(header), 2)]
        genes = []
        for row in reader:
            if not row or not row[0]:
                continue
            desc, gene_id = row[0], row[1]
            values = [float(row[i]) for (_pid, i) in patient_cols]
            genes.append((gene_id, desc, values))
    patient_ids = [pid for (pid, _i) in patient_cols]
    return patient_ids, genes

labels = load_labels(f"{DATA_DIR}/actual.csv")
train_patients, train_genes = load_expression(f"{DATA_DIR}/data_set_ALL_AML_train.csv")
test_patients, test_genes = load_expression(f"{DATA_DIR}/data_set_ALL_AML_independent.csv")

n_genes, n_train, n_test = len(train_genes), len(train_patients), len(test_patients)
print(f"{n_genes} genes x {n_train} train patients + {n_test} independent test patients")
print("train label split:", {lbl: sum(1 for p in train_patients if labels[p] == lbl) for lbl in ("ALL", "AML")})
print("test label split: ", {lbl: sum(1 for p in test_patients if labels[p] == lbl) for lbl in ("ALL", "AML")})

7129 genes x 38 train patients + 34 independent test patients
train label split: {'ALL': 27, 'AML': 11}
test label split:  {'ALL': 20, 'AML': 14}

2. Load into linaldb as real datasets

Two orientations of the same data, both as hybrid scalar+Vector tables:

• Patient-major (train_patients / test_patients): one row per patient, a Vector(7129) column holding that patient's full expression profile — used for classification in Section 5.
• Gene-major (train_genes): one row per gene, a Vector(38) column holding that gene's expression across the 38 training patients — the orientation marker-gene scoring in Section 3 needs.

def to_patient_major(patient_ids, genes):
    out = {pid: [0.0] * len(genes) for pid in patient_ids}
    for gi, (_gid, _desc, values) in enumerate(genes):
        for pi, pid in enumerate(patient_ids):
            out[pid][gi] = values[pi]
    return out


train_patient_major = to_patient_major(train_patients, train_genes)
test_patient_major = to_patient_major(test_patients, test_genes)

db.execute(f"DATASET train_patients COLUMNS (patient_id: Int, label: String, expression: Vector({n_genes}))")
db.execute(f"DATASET test_patients COLUMNS (patient_id: Int, label: String, expression: Vector({n_genes}))")

for pid in train_patients:
    db.execute(f'INSERT INTO train_patients VALUES ({pid}, "{labels[pid]}", {train_patient_major[pid]})')
for pid in test_patients:
    db.execute(f'INSERT INTO test_patients VALUES ({pid}, "{labels[pid]}", {test_patient_major[pid]})')

print(f"inserted {n_train} train + {n_test} test patient rows, each Vector({n_genes})")

inserted 38 train + 34 test patient rows, each Vector(7129)

db.execute(f"DATASET train_genes COLUMNS (gene_id: String, description: String, profile: Vector({n_train}))")

for gid, desc, values in train_genes:
    desc_escaped = desc.replace('"', "'")
    db.execute(f'INSERT INTO train_genes VALUES ("{gid}", "{desc_escaped}", {values})')

print(f"inserted {n_genes} gene rows, each Vector({n_train})")
db.execute("SELECT COUNT(*) AS n FROM train_genes").rows

inserted 7129 gene rows, each Vector(38)

[[7129]]

3. Marker-gene discovery via linaldb's native CORRELATE

For each of the 7,129 genes: hand its training-set expression profile and the ALL/AML label signal directly to linaldb's CORRELATE, and read back the genuine Pearson correlation between that gene's expression and the disease subtype — no manual centering needed, CORRELATE does it internally. This directly reproduces Golub et al.'s own "neighborhood analysis" marker-selection technique — genes whose expression correlates strongly (positive or negative) with the label are the useful discriminators.

train_label_signal = [1.0 if labels[p] == "AML" else -1.0 for p in train_patients]
db.execute(f"VECTOR label_signal = {train_label_signal}")

gene_profiles = db.execute("SELECT gene_id, description, profile FROM train_genes").rows

t0 = time.time()
scores = []
for gid, desc, profile in gene_profiles:
    db.execute(f"VECTOR gene_vec = {profile}")
    db.execute("LET gene_corr = CORRELATE gene_vec WITH label_signal")
    corr = db.execute("SHOW gene_corr").data[0]
    scores.append((gid, desc, corr))

print(f"scored {len(scores)} genes via {len(scores) * 3} linaldb calls in {time.time() - t0:.2f}s")

scored 7129 genes via 21387 linaldb calls in 0.33s

# Push the scores back into linaldb and let it do the ranking, rather than sorting client-side.
db.execute("DATASET marker_scores COLUMNS (gene_id: String, description: String, correlation: Float)")
for gid, desc, corr in scores:
    desc_escaped = desc.replace('"', "'")
    db.execute(f'INSERT INTO marker_scores VALUES ("{gid}", "{desc_escaped}", {corr})')

TOP_K = 25  # per direction -> up to 50 marker genes total
top_aml = db.execute(f"SELECT gene_id, description, correlation FROM marker_scores ORDER BY correlation DESC LIMIT {TOP_K}").rows
top_all = db.execute(f"SELECT gene_id, description, correlation FROM marker_scores ORDER BY correlation ASC LIMIT {TOP_K}").rows

print("Top AML-associated genes (highest positive correlation):")
for gid, desc, corr in top_aml[:5]:
    print(f"  {gid:20s} {corr:+.3f}  {desc}")
print("Top ALL-associated genes (most negative correlation):")
for gid, desc, corr in top_all[:5]:
    print(f"  {gid:20s} {corr:+.3f}  {desc}")

marker_rows = {(g, d, c) for g, d, c in top_aml} | {(g, d, c) for g, d, c in top_all}
marker_ids = sorted({g for g, _d, _c in marker_rows})
print(f"\nmarker gene set size: {len(marker_ids)}")

Top AML-associated genes (highest positive correlation):
  U50136_rna1_at       +0.828  Leukotriene C4 synthase (LTC4S) gene
  X95735_at            +0.822  Zyxin
  M55150_at            +0.811  FAH Fumarylacetoacetate
  M16038_at            +0.777  LYN V-yes-1 Yamaguchi sarcoma viral related oncogene homolog
  Y12670_at            +0.776  LEPR Leptin receptor
Top ALL-associated genes (most negative correlation):
  U22376_cds2_s_at     -0.653  C-myb gene extracted from Human (c-myb) gene; complete primary cds; and five complete alternatively spliced cds
  M86406_at            -0.638  ACTN2 Actinin alpha 2 
  U37055_rna1_s_at     -0.622  Hepatocyte growth factor-like protein gene
  D38128_at            -0.619  PTGIR Prostaglandin I2 (prostacyclin) receptor (IP)
  X15414_at            -0.612  ALDR1 Aldehyde reductase 1 (low Km aldose reductase)

marker gene set size: 50

4. Reduced feature vectors from the marker genes

Re-slice the same Python-side patient-major data down to just the marker genes (plain list indexing — again, no computation) and load a second pair of datasets with a much smaller Vector column, to compare against the full 7,129-dim classifier in Section 5.

gene_index = {gid: i for i, (gid, _desc, _v) in enumerate(train_genes)}
marker_indices = [gene_index[g] for g in marker_ids]
K = len(marker_indices)

db.execute(f"DATASET train_patients_reduced COLUMNS (patient_id: Int, label: String, expression: Vector({K}))")
db.execute(f"DATASET test_patients_reduced COLUMNS (patient_id: Int, label: String, expression: Vector({K}))")

for pid in train_patients:
    reduced = [train_patient_major[pid][i] for i in marker_indices]
    db.execute(f'INSERT INTO train_patients_reduced VALUES ({pid}, "{labels[pid]}", {reduced})')
for pid in test_patients:
    reduced = [test_patient_major[pid][i] for i in marker_indices]
    db.execute(f'INSERT INTO test_patients_reduced VALUES ({pid}, "{labels[pid]}", {reduced})')

print(f"built train/test patient datasets with Vector({K}) reduced expression")

built train/test patient datasets with Vector(50) reduced expression

5. Nearest-centroid classification: full 7,129-dim vs. 50-marker-gene reduced

Same technique linal-db-rs/examples/hdf5_digit_classification.lnl uses for handwritten digits, applied here to leukemia subtypes: per-class centroids via AVG_VEC + GROUP BY, a cosine-similarity JOIN against the held-out test patients, ROW_NUMBER() to keep each patient's single best-matching centroid, and CASE/SUM/COUNT for accuracy.

def classify(train_table, test_table, threshold=0.0):
    cte = (
        f"WITH centroids AS (SELECT label, AVG_VEC(expression) AS centroid FROM {train_table} GROUP BY label), "
        f"classified AS (SELECT t.patient_id AS patient_id, t.label AS true_label, c.label AS predicted_label, "
        f"COSINE_SIM(t.expression, c.centroid) AS similarity, "
        f"ROW_NUMBER() OVER (PARTITION BY patient_id ORDER BY similarity DESC) AS rn "
        f"FROM {test_table} t JOIN centroids c ON COSINE_SIM(t.expression, c.centroid) > {threshold}) "
    )
    preds = db.execute(cte + "SELECT patient_id, true_label, predicted_label, similarity "
                              "FROM classified WHERE rn = 1 ORDER BY patient_id").rows
    overall = db.execute(cte + "SELECT COUNT(*) AS total, "
                                "SUM(CASE WHEN predicted_label = true_label THEN 1 ELSE 0 END) AS correct "
                                "FROM classified WHERE rn = 1").rows[0]
    per_class = db.execute(cte + "SELECT true_label, COUNT(*) AS n, "
                                  "SUM(CASE WHEN predicted_label = true_label THEN 1 ELSE 0 END) AS correct "
                                  "FROM classified WHERE rn = 1 GROUP BY true_label ORDER BY true_label").rows
    return preds, overall, per_class


preds_full, overall_full, per_class_full = classify("train_patients", "test_patients")
total_f, correct_f = overall_full
print(f"FULL {n_genes}-dim:     {correct_f}/{total_f} correct ({100 * correct_f / total_f:.1f}%)")
for lbl, n, c in per_class_full:
    print(f"  {lbl}: {c}/{n}")

FULL 7129-dim:     32/34 correct (94.1%)
  ALL: 20/20
  AML: 12/14

preds_reduced, overall_reduced, per_class_reduced = classify("train_patients_reduced", "test_patients_reduced")
total_r, correct_r = overall_reduced
print(f"REDUCED {K}-dim (markers): {correct_r}/{total_r} correct ({100 * correct_r / total_r:.1f}%)")
for lbl, n, c in per_class_reduced:
    print(f"  {lbl}: {c}/{n}")

majority_baseline = max(sum(1 for _p, t, _pr, _s in preds_full if t == lbl) for lbl in ("ALL", "AML")) / total_f
print(f"\nmajority-class baseline on this test split: {100 * majority_baseline:.1f}%")

REDUCED 50-dim (markers): 31/34 correct (91.2%)
  ALL: 18/20
  AML: 13/14

majority-class baseline on this test split: 58.8%

6. Persistence

Same SAVE DATASET / SHOW DATASET METADATA conventions as 01_linaldb_quickstart.ipynb — here applied to real, non-trivial datasets (7,129-element vectors, 7,129 marker-score rows) rather than a 6-row synthetic table.

db.execute("SAVE DATASET train_patients")
db.execute("SAVE DATASET test_patients")
db.execute("SAVE DATASET marker_scores")

print("on-disk package:", db.dataset_dir("train_patients"))
db.execute("SHOW DATASET METADATA marker_scores")

on-disk package: ./data/default/datasets/train_patients

'=== Dataset Metadata: marker_scores (In-Memory/Legacy) ===\nVersion: 1\nOrigin: Created\nCreated: 2026-09-16T13:45:58.954890Z\nUpdated: 2026-09-16T13:45:58.970253Z\nRows: 7129\n================================'

7. Two charts (display only — matplotlib touches nothing but already-computed results)

import matplotlib.pyplot as plt

top_markers = sorted(marker_rows, key=lambda r: abs(r[2]), reverse=True)[:15]
names = [g for g, _d, _c in top_markers][::-1]
corrs = [c for _g, _d, c in top_markers][::-1]
colors = ["#c0392b" if c < 0 else "#2980b9" for c in corrs]

fig, ax = plt.subplots(figsize=(6, 5))
ax.barh(names, corrs, color=colors)
ax.axvline(0, color="black", linewidth=0.8)
ax.set_xlabel("Pearson correlation with AML label (linaldb CORRELATE)")
ax.set_title("Top marker genes by |correlation| with ALL/AML")
plt.tight_layout()
plt.show()

fig, ax = plt.subplots(figsize=(4, 3))
labels_ = [f"full\n({n_genes}-dim)", f"markers\n({K}-dim)", "majority\nbaseline"]
values = [100 * correct_f / total_f, 100 * correct_r / total_r, 100 * majority_baseline]
ax.bar(labels_, values, color=["#7f8c8d", "#27ae60", "#bdc3c7"])
ax.set_ylabel("test-set accuracy (%)")
ax.set_ylim(0, 100)
ax.set_title("ALL vs AML classification accuracy")
for i, v in enumerate(values):
    ax.text(i, v + 1.5, f"{v:.1f}%", ha="center")
plt.tight_layout()
plt.show()

8. What this exercised

• Hybrid Vector-column datasets at real scale: Vector(7129) patient rows, Vector(38) gene-profile rows, bulk INSERT of thousands of real (not synthetic) numeric vectors.
• CORRELATE used natively for marker-gene discovery across all 7,129 genes, and MEAN + vector subtraction for the Section 0 sanity check — both fixed as part of building this exact notebook (see below).
• AVG_VEC for per-class centroids, a cosine-similarity JOIN, ROW_NUMBER() OVER (PARTITION BY ...) for top-1 matching, CASE/SUM/COUNT/GROUP BY for accuracy scoring — the same nearest-centroid idiom as linal-db-rs/examples/hdf5_digit_classification.lnl, applied to a real 2-class biological classification problem instead of 10-class digit recognition.
• SAVE DATASET / SHOW DATASET METADATA persistence on the resulting real datasets.
• No pandas, numpy, or sklearn import anywhere in this notebook; matplotlib appears only in Section 7, rendering results linaldb already computed.

Two real engine bugs, found and fixed: building this notebook surfaced two real, reproducible discrepancies between docs/DSL_REFERENCE.md and the engine's actual behavior — CORRELATE computing a raw, unnormalized dot product instead of the documented Pearson correlation, and SUM/MEAN/STDEV returning the wrong tensor shape internally (a rank-1 Vector(1) instead of a true rank-0 scalar), which silently corrupted vector - MEAN(vector) (and similar ops against a reduction result) past the first element. Both were fixed in gorigami/linaldb#81, shipped as engine v0.1.79 and PyPI linaldb 0.1.2 — the exact version this notebook now runs against. Same recurring project lesson: a real end-to-end workflow against real data finds bugs isolated unit tests miss, the same way linal-hub's pytest suite found the 4 bugs fixed in v0.1.77.