Altered Cancer Pathways - QUBONotebooks: Quantum Integer Programming

Maintained by the JuliaQUBO organization
SECQUOIA · PSR Energy

An original Julia case study balancing mutation coverage and pairwise co-mutation.

Educational use only. This notebook reproduces an optimization formulation on a fixed public-data aggregate. It does not identify a clinically validated pathway and is not medical guidance.

Setup¶

Local installation¶

From the repository root, instantiate the shared Julia environment before opening this notebook:

julia --project=notebooks_jl -e 'using Pkg; Pkg.instantiate()'

Restart the kernel and run every cell with the targeted command:

make verify-cancer-genomics-julia

Ordinary execution reads only the committed files under notebooks_data. It requires no account, credential, live data request, or quantum service. The optional make refresh-tcga-aml maintenance command is documented in notebooks_data/README.md and is not run here.

Google Colab¶

Open the badge above, select a Julia runtime, and run the setup cells. The bootstrap clones this repository only when Colab does not already have it, then activates the same checked-in notebook project and reads the same committed aggregate.

function load_qubonotebooks_bootstrap()
    candidates = (
        joinpath(pwd(), "scripts", "notebook_bootstrap.jl"),
        joinpath(pwd(), "..", "scripts", "notebook_bootstrap.jl"),
        joinpath(pwd(), "QUBONotebooks", "scripts", "notebook_bootstrap.jl"),
        joinpath("/content", "QUBONotebooks", "scripts", "notebook_bootstrap.jl"),
    )

    for candidate in candidates
        if isfile(candidate)
            include(candidate)
            return nothing
        end
    end

    in_colab = haskey(ENV, "COLAB_RELEASE_TAG") || haskey(ENV, "COLAB_JUPYTER_IP") || isdir(joinpath("/content", "sample_data"))
    if in_colab
        repo_dir = get(ENV, "QUBONOTEBOOKS_REPO_DIR", joinpath(pwd(), "QUBONotebooks"))
        if !isdir(repo_dir)
            println("[bootstrap] Cloning JuliaQUBO/QUBONotebooks into $repo_dir")
            run(Cmd(["git", "clone", "--quiet", "--depth", "1", "https://github.com/JuliaQUBO/QUBONotebooks.git", repo_dir]))
        end
        include(joinpath(repo_dir, "scripts", "notebook_bootstrap.jl"))
        return nothing
    end

    error("Could not locate scripts/notebook_bootstrap.jl from $(pwd()).")
end

load_qubonotebooks_bootstrap()

BOOTSTRAP = Base.invokelatest(QUBONotebooksBootstrap.bootstrap_notebook, "9-CancerGenomics")
QUBONOTEBOOKS_REPO_DIR = BOOTSTRAP.repo_dir
JULIA_NOTEBOOKS_DIR = BOOTSTRAP.notebooks_dir
JULIA_PROJECT_DIR = BOOTSTRAP.project_dir
IN_COLAB = BOOTSTRAP.in_colab;

import Pkg

python_warning_filter = "ignore:invalid escape sequence:SyntaxWarning"
python_warning_filters = String.(filter(!isempty, split(get(ENV, "PYTHONWARNINGS", ""), ",")))
if python_warning_filter ∉ python_warning_filters
    push!(python_warning_filters, python_warning_filter)
    ENV["PYTHONWARNINGS"] = join(python_warning_filters, ",")
end

if @isdefined(JULIA_PROJECT_DIR)
    Pkg.activate(JULIA_PROJECT_DIR; io = devnull)
else
    Pkg.activate(@__DIR__; io = devnull)
end
Pkg.instantiate(; io = devnull, allow_autoprecomp = false)

Learning objectives¶

By the end of this notebook you will be able to:

construct coverage and symmetric co-mutation matrices from a patient-by-gene incidence matrix;
explain the double-counting convention in the matrix QUBO;
exhaustively verify a tiny fixture against an independent analytic energy;
decode coverage, distinct-patient coverage, co-mutation, pathway size, and raw energy safely; and
run and validate a seeded local sampler on a provenance-checked TCGA AML aggregate.

Prerequisites¶

Prior notebooks: Notebook 2 introduces JuMP/QUBO models, and Notebook 7 develops exhaustive verification. This case study defines the helpers it needs so it remains independently runnable.

Mathematical background: Binary incidence matrices, matrix multiplication, quadratic objectives, and finite enumeration.

Data background: The committed TCGA AML files are gene-level reproducibility aggregates, not raw clinical records. Patient and sample identifiers are intentionally absent.

Software and accounts: Julia 1.10+ with the shared notebook project instantiated. No external account is required.

# QUBONOTEBOOKS_COLAB_IMPORT_CELL
Base.invokelatest(
    QUBONotebooksBootstrap.warm_notebook_packages!,
    "9-CancerGenomics",
);

From incidence data to coverage and co-mutation¶

Let rows of the binary matrix $B$ represent patients or samples and columns represent genes. An entry $B_{p,g}=1$ means gene $g$ is observed as mutated for patient $p$ . The cited formulation uses the transposed orientation (genes by patients); either convention produces the same gene-by-gene Gram matrix when used consistently.

We start with seven synthetic patients and four synthetic genes. The fixture is deliberately small enough to inspect without biological interpretation.

tiny_patients = ["P1", "P2", "P3", "P4", "P5", "P6", "P7"]
tiny_genes = ["G-A", "G-B", "G-C", "G-D"]
tiny_B = Int[
    1  1  0  0
    1  0  0  0
    1  0  0  0
    0  1  0  0
    0  0  1  1
    0  0  1  0
    0  0  0  1
]

@assert size(tiny_B) == (length(tiny_patients), length(tiny_genes))
@assert all(value in (0, 1) for value in tiny_B)
println("Tiny fixture: $(size(tiny_B, 1)) patients × $(size(tiny_B, 2)) genes")

Tiny fixture: 7 patients × 4 genes

For each gene, the coverage value is the number of patients containing that gene. The diagonal matrix $D$ stores those counts. For each pair $i<j$ , $A_{ij}$ is the number of patients containing both genes. We assign both $A_{ij}$ and $A_{ji}$ , and leave the diagonal at zero. Thus

B^{\mathsf{T}}B = D + A.

(1)

The explicit pair loop below prevents the stale-index defect found in the companion notebooks.

"""
    coverage_and_comutation(B)

Return the per-gene coverage vector and zero-diagonal symmetric co-mutation
matrix for a binary patient-by-gene incidence matrix.
"""
function coverage_and_comutation(B::AbstractMatrix{<:Integer})
    @assert all(value in (0, 1) for value in B)
    n_genes = size(B, 2)
    coverage = vec(sum(B; dims = 1))
    A = zeros(Int, n_genes, n_genes)

    for i in 1:(n_genes - 1)
        for j in (i + 1):n_genes
            overlap = sum(B[:, i] .* B[:, j])
            A[i, j] = overlap
            A[j, i] = overlap
        end
    end

    return coverage, A
end

coverage_and_comutation

tiny_coverage, tiny_A = coverage_and_comutation(tiny_B)
tiny_D = Diagonal(tiny_coverage)
tiny_gram = tiny_B' * tiny_B
expected_tiny_A = Int[
    0  1  0  0
    1  0  0  0
    0  0  0  1
    0  0  1  0
]

@assert tiny_coverage == [3, 2, 2, 2]
@assert tiny_A == expected_tiny_A
@assert tiny_A == tiny_A'
@assert iszero(diag(tiny_A))
@assert tiny_gram == tiny_D + tiny_A

println("D diagonal = $tiny_coverage")
println("A =")
display(tiny_A)

D diagonal = [3, 2, 2, 2]
A =

4×4 Matrix{Int64}:
 0  1  0  0
 1  0  0  0
 0  0  0  1
 0  0  1  0

Build and interpret the QUBO¶

For binary $x_g$ , where one selects gene $g$ , the tutorial model is

Q(x)=x^{\mathsf{T}}(A-\alpha D)x = 2\sum_{i<j}A_{ij}x_ix_j-\alpha\sum_iD_{ii}x_i.

(2)

The negative diagonal term rewards gene coverage. The positive off-diagonal term penalizes pairwise co-mutation. Because $A$ is symmetric, $x^{\mathsf{T}}Ax$ intentionally counts each undirected gene pair twice; code, analytic energy, and reported metrics all use that convention.

The named parameter $\alpha$ controls the modeling trade-off. It is not a confidence score and does not make the sampled genes clinically meaningful.

function all_binary_states(n::Integer)
    return [
        collect(state)
        for state in Iterators.product(ntuple(_ -> (0, 1), n)...)
    ]
end

function pathway_energy(coverage, A, bits; alpha)
    @assert alpha >= 0
    @assert length(coverage) == size(A, 1) == size(A, 2) == length(bits)
    return dot(bits, A * bits) - alpha * dot(coverage, bits)
end

function distinct_patient_coverage_bounds(coverage, A, selected, patient_count)
    isempty(selected) && return (lower = 0, upper = 0)

    coverage_sum = sum(coverage[selected])
    pairwise_overlap = sum(A[i, j] for i in selected for j in selected if i < j)
    lower = max(maximum(coverage[selected]), coverage_sum - pairwise_overlap)
    upper = min(patient_count, coverage_sum)

    @assert 0 <= lower <= upper <= patient_count
    return (lower = lower, upper = upper)
end

"""
    decoded_pathway_metrics(gene_names, coverage, A, bits; incidence=nothing, patient_count)

Decode a pathway without dividing by the selected size or overlap count. Exact
distinct-patient coverage is returned only when an incidence matrix is present;
otherwise the aggregate supports rigorous Bonferroni lower and upper bounds.
"""
function decoded_pathway_metrics(
    gene_names,
    coverage,
    A,
    bits;
    incidence = nothing,
    patient_count,
)
    selected = findall(==(1), bits)
    pathway_size = length(selected)
    coverage_sum = sum(coverage[selected])
    pairwise_comutation = sum(
        (A[i, j] for i in selected for j in selected if i < j);
        init = 0,
    )
    double_counted_comutation = dot(bits, A * bits)
    pair_count = pathway_size * (pathway_size - 1) ÷ 2
    zero_comutation_pairs = count(
        A[i, j] == 0 for i in selected for j in selected if i < j
    )
    zero_comutation_pair_fraction =
        pair_count == 0 ? missing : zero_comutation_pairs / pair_count

    if isnothing(incidence)
        exact_unique_coverage = missing
        unique_coverage_bounds = distinct_patient_coverage_bounds(
            coverage,
            A,
            selected,
            patient_count,
        )
    else
        @assert size(incidence, 2) == length(gene_names)
        exact_unique_coverage = isempty(selected) ? 0 : count(
            >(0),
            vec(sum(incidence[:, selected]; dims = 2)),
        )
        unique_coverage_bounds = (
            lower = exact_unique_coverage,
            upper = exact_unique_coverage,
        )
    end

    @assert double_counted_comutation == 2 * pairwise_comutation

    return (
        selected_indices = selected,
        selected_genes = gene_names[selected],
        pathway_size = pathway_size,
        coverage_sum = coverage_sum,
        exact_unique_coverage = exact_unique_coverage,
        unique_coverage_bounds = unique_coverage_bounds,
        pairwise_comutation = pairwise_comutation,
        double_counted_comutation = double_counted_comutation,
        zero_comutation_pair_fraction = zero_comutation_pair_fraction,
    )
end

function validate_decoded_metrics(metrics, coverage, A, bits; patient_count)
    @assert metrics.selected_indices == findall(==(1), bits)
    @assert metrics.pathway_size == sum(bits)
    @assert metrics.coverage_sum == dot(coverage, bits)
    @assert metrics.double_counted_comutation == dot(bits, A * bits)
    @assert metrics.double_counted_comutation == 2 * metrics.pairwise_comutation
    @assert 0 <= metrics.unique_coverage_bounds.lower <=
        metrics.unique_coverage_bounds.upper <= patient_count
    if !ismissing(metrics.exact_unique_coverage)
        @assert metrics.unique_coverage_bounds.lower ==
            metrics.exact_unique_coverage ==
            metrics.unique_coverage_bounds.upper
    end
    return true
end

validate_decoded_metrics (generic function with 1 method)

function build_pathway_model(coverage, A; alpha)
    @assert alpha >= 0
    @assert length(coverage) == size(A, 1) == size(A, 2)
    @assert A == A'
    @assert iszero(diag(A))

    model = Model()
    @variable(model, pathway_x[1:length(coverage)], Bin)
    @objective(
        model,
        Min,
        sum(
            A[i, j] * pathway_x[i] * pathway_x[j]
            for i in eachindex(coverage), j in eachindex(coverage)
        ) - alpha * sum(
            coverage[i] * pathway_x[i] for i in eachindex(coverage)
        ),
    )
    return model, pathway_x
end

function evaluate_jump_objective(model::Model, variables, bits)
    assignment = Dict(variables[i] => Float64(bits[i]) for i in eachindex(variables))
    return JuMP.value(variable -> assignment[variable], JuMP.objective_function(model))
end

function exact_sampler_optima!(model::Model, variables; atol = 1e-9)
    set_optimizer(model, QUBO.ExactSampler.Optimizer)
    optimize!(model)
    @assert termination_status(model) in (JuMP.MOI.OPTIMAL, JuMP.MOI.LOCALLY_SOLVED)

    energies = [objective_value(model; result = result) for result in 1:result_count(model)]
    best_energy = minimum(energies)
    optima = [
        round.(Int, value.(variables; result = result))
        for result in eachindex(energies)
        if isapprox(energies[result], best_energy; atol = atol, rtol = 0)
    ]
    return best_energy, optima
end

same_states(left, right) = Set(Tuple.(left)) == Set(Tuple.(right))

same_states (generic function with 1 method)

tiny_alpha = 0.75
tiny_model, tiny_x = build_pathway_model(tiny_coverage, tiny_A; alpha = tiny_alpha)
nothing

Exhaustive validation on the tiny fixture¶

Four genes give $2^4=16$ assignments. For every state, we compare the JuMP objective against the independent matrix calculation and against the explicit pairwise expansion. This proves the factor-of-two convention before any sampler is trusted.

tiny_states = all_binary_states(length(tiny_genes))
tiny_jump_energies = [
    evaluate_jump_objective(tiny_model, tiny_x, bits) for bits in tiny_states
]
tiny_direct_energies = [
    pathway_energy(tiny_coverage, tiny_A, bits; alpha = tiny_alpha)
    for bits in tiny_states
]
tiny_pair_expansion_energies = [
    2 * sum(
        tiny_A[i, j] * bits[i] * bits[j]
        for i in eachindex(bits), j in eachindex(bits) if i < j
    ) - tiny_alpha * sum(tiny_coverage .* bits)
    for bits in tiny_states
]

@assert length(tiny_states) == 16
@assert all(isapprox.(tiny_jump_energies, tiny_direct_energies; atol = 1e-9, rtol = 0))
@assert tiny_pair_expansion_energies == tiny_direct_energies

println("All 16 assignments match the analytic matrix and explicit pairwise energies.")

All 16 assignments match the analytic matrix and explicit pairwise energies.

tiny_enumerated_best_energy = minimum(tiny_direct_energies)
tiny_enumerated_optima = tiny_states[findall(==(tiny_enumerated_best_energy), tiny_direct_energies)]

tiny_exact_best_energy, tiny_exact_optima = exact_sampler_optima!(tiny_model, tiny_x)
@assert tiny_enumerated_best_energy == -3.75
@assert Set(Tuple.(tiny_enumerated_optima)) == Set(((1, 0, 1, 0), (1, 0, 0, 1)))
@assert isapprox(tiny_exact_best_energy, tiny_enumerated_best_energy; atol = 1e-9, rtol = 0)
@assert same_states(tiny_exact_optima, tiny_enumerated_optima)

println("Enumerated optimum energy: $tiny_enumerated_best_energy")
println("ExactSampler returns the same $(length(tiny_exact_optima)) global minima: $tiny_exact_optima")

Enumerated optimum energy: -3.75
ExactSampler returns the same 2 global minima: [[1, 0, 0, 1], [1, 0, 1, 0]]

Decode pathway metrics safely¶

The coverage reward is the sum of per-gene patient counts and can count a patient more than once. Distinct-patient coverage is the union of affected patients. The tiny incidence matrix lets us compute that union exactly. Pairwise co-mutation is reported once per undirected gene pair, while the QUBO uses twice that value.

A zero- or one-gene pathway has no gene pairs, so the pair fraction is reported as missing instead of dividing by zero.

tiny_representative = first(tiny_enumerated_optima)
tiny_metrics = decoded_pathway_metrics(
    tiny_genes,
    tiny_coverage,
    tiny_A,
    tiny_representative;
    incidence = tiny_B,
    patient_count = length(tiny_patients),
)
@assert validate_decoded_metrics(
    tiny_metrics,
    tiny_coverage,
    tiny_A,
    tiny_representative;
    patient_count = length(tiny_patients),
)
@assert tiny_metrics.selected_genes == ["G-A", "G-C"]
@assert tiny_metrics.pathway_size == 2
@assert tiny_metrics.coverage_sum == 5
@assert tiny_metrics.exact_unique_coverage == 5
@assert tiny_metrics.pairwise_comutation == 0
@assert tiny_metrics.zero_comutation_pair_fraction == 1.0

println("Selected genes: $(join(tiny_metrics.selected_genes, ", "))")
println("Coverage sum: $(tiny_metrics.coverage_sum); distinct-patient coverage: $(tiny_metrics.exact_unique_coverage)")
println("Pairwise co-mutation: $(tiny_metrics.pairwise_comutation); raw energy: $(pathway_energy(tiny_coverage, tiny_A, tiny_representative; alpha = tiny_alpha))")

Selected genes: G-A, G-C
Coverage sum: 5; distinct-patient coverage: 5
Pairwise co-mutation: 0; raw energy: -3.75

tiny_empty = zeros(Int, length(tiny_genes))
tiny_empty_metrics = decoded_pathway_metrics(
    tiny_genes,
    tiny_coverage,
    tiny_A,
    tiny_empty;
    incidence = tiny_B,
    patient_count = length(tiny_patients),
)
@assert validate_decoded_metrics(
    tiny_empty_metrics,
    tiny_coverage,
    tiny_A,
    tiny_empty;
    patient_count = length(tiny_patients),
)
@assert tiny_empty_metrics.pathway_size == 0
@assert tiny_empty_metrics.coverage_sum == 0
@assert tiny_empty_metrics.exact_unique_coverage == 0
@assert tiny_empty_metrics.unique_coverage_bounds == (lower = 0, upper = 0)
@assert ismissing(tiny_empty_metrics.zero_comutation_pair_fraction)
@assert pathway_energy(tiny_coverage, tiny_A, tiny_empty; alpha = tiny_alpha) == 0

println("Empty pathway decoded safely: size=0, coverage=0, co-mutation=0, energy=0.")

Empty pathway decoded safely: size=0, coverage=0, co-mutation=0, energy=0.

Provenance-checked TCGA AML aggregate¶

The full educational example loads three committed CSV files and their JSON provenance record. The data refresh was performed separately by the reviewed workflow in issue #91; this notebook never contacts cBioPortal. The files retain only gene labels, coverage counts, pairwise co-mutation counts, and aggregate provenance. They contain no patient or sample identifiers.

The loader cross-checks gene ordering across all files, integer bounds, symmetry, zero diagonal, and the selected-gene and patient counts recorded in provenance.

function split_csv_lines(path)
    return [split(line, ',') for line in readlines(path)]
end

"""
    load_tcga_aml_aggregate(repo_root)

Load and cross-check the committed cancer-genomics aggregate without network
access. The repository artifact uses a simple comma-delimited schema with no
quoted fields.
"""
function load_tcga_aml_aggregate(repo_root)
    data_dir = joinpath(repo_root, "notebooks_data")
    gene_rows = split_csv_lines(joinpath(data_dir, "9-CancerGenomics_genes.csv"))
    coverage_rows = split_csv_lines(joinpath(data_dir, "9-CancerGenomics_coverage.csv"))
    matrix_rows = split_csv_lines(joinpath(data_dir, "9-CancerGenomics_comutation.csv"))
    provenance = JSON.parsefile(joinpath(data_dir, "9-CancerGenomics_provenance.json"))

    @assert gene_rows[1] == ["rank", "gene", "mutation_record_count", "patient_coverage"]
    @assert coverage_rows[1] == ["gene", "patient_count"]

    genes = String[row[2] for row in gene_rows[2:end]]
    ranks = parse.(Int, [row[1] for row in gene_rows[2:end]])
    gene_file_coverage = parse.(Int, [row[4] for row in gene_rows[2:end]])
    coverage_genes = String[row[1] for row in coverage_rows[2:end]]
    coverage = parse.(Int, [row[2] for row in coverage_rows[2:end]])
    matrix_genes = String.(matrix_rows[1][2:end])

    @assert ranks == collect(eachindex(genes))
    @assert length(genes) == length(unique(genes))
    @assert genes == coverage_genes == matrix_genes
    @assert coverage == gene_file_coverage
    @assert all(>=(0), coverage)

    A = zeros(Int, length(genes), length(genes))
    for (i, row) in enumerate(matrix_rows[2:end])
        @assert row[1] == genes[i]
        @assert length(row) == length(genes) + 1
        A[i, :] = parse.(Int, row[2:end])
    end

    @assert A == A'
    @assert iszero(diag(A))
    @assert all(A .>= 0)
    @assert all(A[i, j] <= min(coverage[i], coverage[j]) for i in eachindex(genes), j in eachindex(genes))

    aggregate = provenance["aggregate"]
    patient_count = Int(aggregate["patient_count"])
    @assert length(genes) == Int(aggregate["selected_gene_count"])
    @assert patient_count > 0

    return (
        genes = genes,
        coverage = coverage,
        D = Diagonal(coverage),
        A = A,
        patient_count = patient_count,
        provenance = provenance,
    )
end

load_tcga_aml_aggregate

aml = load_tcga_aml_aggregate(QUBONOTEBOOKS_REPO_DIR)

@assert length(aml.genes) == 33
@assert aml.patient_count == 196
@assert aml.A == aml.A'
@assert iszero(diag(aml.A))
@assert all(
    aml.A[i, j] <= min(aml.coverage[i], aml.coverage[j])
    for i in eachindex(aml.genes), j in eachindex(aml.genes)
)

println("Loaded $(length(aml.genes)) genes for $(aml.patient_count) aggregate patients.")
println("Source retrieval recorded at $(aml.provenance["retrieved_at"]); default execution made no network request.")

Loaded 33 genes for 196 aggregate patients.
Source retrieval recorded at 2026-07-20T13:05:32Z; default execution made no network request.

What the aggregate can and cannot decode¶

For a selected set $S$ , the files determine the coverage reward $\sum_{g\in S}D_g$ , every pairwise co-mutation $A_{ij}$ , and the QUBO energy exactly. They do not determine the exact union of affected patients for arbitrary $|S|>2$ : triple and higher intersections were intentionally not committed.

We therefore report a rigorous interval. The largest selected-gene coverage and the second-order Bonferroni expression $\sum D_g-\sum_{i<j}A_{ij}$ give lower bounds; the patient count and $\sum D_g$ give upper bounds. This preserves the privacy/minimal-data design of issue #91 and avoids fabricating a distinct-patient count.

function sample_aggregate(
    gene_names,
    coverage,
    A;
    alpha,
    patient_count,
    seed,
    num_reads,
    num_sweeps,
)
    model, variables = build_pathway_model(coverage, A; alpha = alpha)
    set_optimizer(model, DWave.Neal.Optimizer)
    set_optimizer_attribute(model, "num_reads", num_reads)
    set_optimizer_attribute(model, "num_sweeps", num_sweeps)
    set_optimizer_attribute(model, "seed", seed)
    optimize!(model)
    @assert termination_status(model) == JuMP.MOI.LOCALLY_SOLVED

    samples = [
        begin
            bits = round.(Int, value.(variables; result = result))
            reported_energy = objective_value(model; result = result)
            recomputed_energy = pathway_energy(coverage, A, bits; alpha = alpha)
            metrics = decoded_pathway_metrics(
                gene_names,
                coverage,
                A,
                bits;
                patient_count = patient_count,
            )

            @assert isapprox(reported_energy, recomputed_energy; atol = 1e-8, rtol = 0)
            @assert validate_decoded_metrics(
                metrics,
                coverage,
                A,
                bits;
                patient_count = patient_count,
            )

            (
                bits = bits,
                reported_energy = reported_energy,
                recomputed_energy = recomputed_energy,
                metrics = metrics,
            )
        end
        for result in 1:result_count(model)
    ]
    sort!(samples; by = sample -> (sample.reported_energy, Tuple(sample.bits)))

    return samples
end

sample_aggregate (generic function with 1 method)

alpha_cases = [0.25, 0.45, 1.0]
sampler_seed = 94094
sampler_reads = 500
sampler_sweeps = 2_000

aggregate_results = [
    begin
        samples = sample_aggregate(
            aml.genes,
            aml.coverage,
            aml.A;
            alpha = alpha,
            patient_count = aml.patient_count,
            seed = sampler_seed,
            num_reads = sampler_reads,
            num_sweeps = sampler_sweeps,
        )
        best = first(samples)
        (
            alpha = alpha,
            returned_states = length(samples),
            validated_states = length(samples),
            bits = best.bits,
            genes = best.metrics.selected_genes,
            pathway_size = best.metrics.pathway_size,
            coverage_sum = best.metrics.coverage_sum,
            unique_coverage_bounds = best.metrics.unique_coverage_bounds,
            pairwise_comutation = best.metrics.pairwise_comutation,
            raw_energy = best.recomputed_energy,
        )
    end
    for alpha in alpha_cases
]

@assert all(result.returned_states == result.validated_states for result in aggregate_results)
@assert length(unique(result.pathway_size for result in aggregate_results)) > 1
@assert length(unique(Tuple(result.bits) for result in aggregate_results)) > 1

println("Seeded local Neal sampler: seed=$sampler_seed, reads=$sampler_reads, sweeps=$sampler_sweeps")
println("alpha  states checked  size  coverage sum  distinct-patient bounds  pair co-mutation  raw energy")
for result in aggregate_results
    bounds = result.unique_coverage_bounds
    @printf(
        "%4.2f %7d %7d %5d %13d %9d–%-3d %17d %11.2f\n",
        result.alpha,
        result.returned_states,
        result.validated_states,
        result.pathway_size,
        result.coverage_sum,
        bounds.lower,
        bounds.upper,
        result.pairwise_comutation,
        result.raw_energy,
    )
    println("      selected genes: $(join(result.genes, ", "))")
end

Seeded local Neal sampler: seed=94094, reads=500, sweeps=2000
alpha  states checked  size  coverage sum  distinct-patient bounds  pair co-mutation  raw energy

0.25      13      13     4            91        91–91                  0      -22.75
      selected genes: NPM1, RUNX1, TP53, FCGBP
0.45       8       8     6           109       106–109                 3      -43.05
      selected genes: FLT3, RUNX1, TP53, KIT, KRAS, FCGBP
1.00      23      23    11           152       133–152                19     -114.00
      selected genes: FLT3, RUNX1, TP53, CEBPA, IDH1, KIT, KRAS, RAD21, CROCC, CSMD1, FCGBP

The sampled pathway size and composition change across the three $\alpha$ values. These are reproducible outputs from a locked, seeded local sampler, not proofs of global optimality. The validity claim is narrower and checked for every returned state: its bit vector, decoded metrics, and independently recomputed raw energy agree.

Increasing $\alpha$ makes coverage more valuable relative to the pairwise co-mutation penalty, so larger selected sets become more attractive. No clinical interpretation is made from any selected name.

Practice checkpoints¶

For the tiny state selecting G-A and G-B, compute the undirected co-mutation count and predict the contribution of $x^{\mathsf{T}}Ax$ . Confirm the factor of two.
Enumerate the tiny fixture at $\alpha=0.5$ and $\alpha=1.5$ . Compare the optimal pathway sizes without relying on stochastic output.
Select G-A, G-B, and G-C. Derive the aggregate-style distinct-patient bounds, then use the incidence matrix to locate the exact union inside them.

# EXERCISE: select G-A and G-B, then compare one undirected overlap
# with the symmetric matrix contribution x' * A * x.
nothing

# SOLUTION (hidden in workshop version):
overlap_bits = [1, 1, 0, 0]
one_counted_overlap = tiny_A[1, 2]
matrix_overlap = dot(overlap_bits, tiny_A * overlap_bits)
@assert one_counted_overlap == 1
@assert matrix_overlap == 2 * one_counted_overlap == 2
println("One undirected co-mutation contributes $matrix_overlap to the symmetric matrix energy.")

One undirected co-mutation contributes 2 to the symmetric matrix energy.

# EXERCISE: enumerate the tiny fixture at alpha=0.5 and alpha=1.5.
# Compare complete sets of minimizers and their pathway sizes.
nothing

# SOLUTION (hidden in workshop version):
tiny_sensitivity = [
    begin
        energies = [
            pathway_energy(tiny_coverage, tiny_A, bits; alpha = alpha)
            for bits in tiny_states
        ]
        best_energy = minimum(energies)
        optima = tiny_states[findall(==(best_energy), energies)]
        (
            alpha = alpha,
            energy = best_energy,
            sizes = sort(unique(sum(bits) for bits in optima)),
            optima = optima,
        )
    end
    for alpha in (0.5, 1.5)
]
@assert tiny_sensitivity[1].sizes == [2]
@assert tiny_sensitivity[2].sizes == [4]
println("Tiny optimum sizes: alpha=0.5 → $(tiny_sensitivity[1].sizes), alpha=1.5 → $(tiny_sensitivity[2].sizes)")

Tiny optimum sizes: alpha=0.5 → [2], alpha=1.5 → [4]

# EXERCISE: select G-A, G-B, and G-C.
# Compute the union bounds from D and A before checking the exact incidence rows.
nothing

# SOLUTION (hidden in workshop version):
bound_bits = [1, 1, 1, 0]
bounded_metrics = decoded_pathway_metrics(
    tiny_genes,
    tiny_coverage,
    tiny_A,
    bound_bits;
    patient_count = length(tiny_patients),
)
exact_metrics = decoded_pathway_metrics(
    tiny_genes,
    tiny_coverage,
    tiny_A,
    bound_bits;
    incidence = tiny_B,
    patient_count = length(tiny_patients),
)
@assert bounded_metrics.unique_coverage_bounds == (lower = 6, upper = 7)
@assert exact_metrics.exact_unique_coverage == 6
@assert bounded_metrics.unique_coverage_bounds.lower <=
    exact_metrics.exact_unique_coverage <=
    bounded_metrics.unique_coverage_bounds.upper
println("Aggregate-style bounds: 6–7 patients; incidence-derived exact union: 6 patients.")

Aggregate-style bounds: 6–7 patients; incidence-derived exact union: 6 patients.

Summary¶

Learning objectives met:

A patient-by-gene incidence matrix yields coverage $D$ and a zero-diagonal symmetric co-mutation matrix $A$ by visiting every $i<j$ pair.
The QUBO $x^{\mathsf{T}}(A-\alpha D)x$ counts each undirected co-mutation twice, consistently in code and mathematics.
All 16 tiny states agree across the analytic matrix formula, explicit pair expansion, JuMP objective, and ExactSampler optimum set.
Empty-pathway metrics are safe, and exact distinct-patient coverage is computed only when incidence rows are available.
The committed 33-gene aggregate loads without a service call; every state returned by the seeded local sampler passes independent energy and metric checks.
Aggregate distinct-patient coverage is reported as a rigorous interval because higher-order patient intersections are not present.

Next steps: Notebook 11 will compare solver interfaces and optional D-Wave hardware. This modeling notebook intentionally stops at a local, credential-free sampler and makes no performance claim.

Further reading:

The Five Starter Problems tutorial gives the displayed QUBO formulation and its coverage/exclusivity motivation.
Alghassi et al. develop quantum and quantum-inspired formulations for altered-pathway discovery.
Ley et al. describe the underlying TCGA AML study; the committed provenance explains the exact cBioPortal aggregate used here.

References¶

A. R. Mazumder and S. Tayur, Five Starter Problems: Solving Quadratic Unconstrained Binary Optimization Models on Quantum Computers, TutORials in Operations Research (2025), pp. 145–183, Mazumder & Tayur (2025).
H. Alghassi, R. Dridi, A. G. Robertson, and S. Tayur, Quantum and Quantum-inspired Methods for de novo Discovery of Altered Cancer Pathways, bioRxiv 845719 (2019), Alghassi et al. (2019).
T. J. Ley et al., Genomic and Epigenomic Landscapes of Adult De Novo Acute Myeloid Leukemia, New England Journal of Medicine 368 (2013), pp. 2059–2074, Massachusetts Medical Society (2013).
E. Cerami et al., The cBio Cancer Genomics Portal: An Open Platform for Exploring Multidimensional Cancer Genomics Data, Cancer Discovery 2 (2012), pp. 401–404, Cerami et al. (2012); current portal: https://www.cbioportal.org/.
Companion materials: https://github.com/arulrhikm/Solving-QUBOs-on-Quantum-Computers.
D-Wave Neal simulated annealing documentation: https://docs.dwavequantum.com/en/latest/ocean/api_ref_samplers/index.html.

The committed artifact’s exact query, retrieval time, aggregation rules, license, and checksums are recorded in notebooks_data/9-CancerGenomics_provenance.json. This notebook contains original Julia code and prose. The unlicensed companion repository is cited as context; no cells, prose, outputs, figures, or assets were copied.

This educational reproduction is not a clinically validated pathway analysis, biomarker study, diagnostic tool, or medical guidance.

References¶

Mazumder, A. R., & Tayur, S. (2025). Five Starter Problems: Solving Quadratic Unconstrained Binary Optimization Models on Quantum Computers. In Tutorials in Operations Research: Advances in Analytics and Operations Research: Improving Decisions to Secure the Future (pp. 145–183). INFORMS. 10.1287/educ.2025.0288
Alghassi, H., Dridi, R., Robertson, A. G., & Tayur, S. (2019). Quantum and Quantum-inspired Methods for de novo Discovery of Altered Cancer Pathways. 10.1101/845719
(2013). New England Journal of Medicine, 368(22), 2059–2074. 10.1056/nejmoa1301689
Cerami, E., Gao, J., Dogrusoz, U., Gross, B. E., Sumer, S. O., Aksoy, B. A., Jacobsen, A., Byrne, C. J., Heuer, M. L., Larsson, E., Antipin, Y., Reva, B., Goldberg, A. P., Sander, C., & Schultz, N. (2012). The cBio Cancer Genomics Portal: An Open Platform for Exploring Multidimensional Cancer Genomics Data. Cancer Discovery, 2(5), 401–404. 10.1158/2159-8290.cd-12-0095