Train a classifier

SigMA2's signature classifiers cover 14 tumor types on 6 platforms. For a tumor type, a signature or a panel it does not cover, you train your own: from whole genomes with known signature exposures, each genome is labeled from its exposure, reduced to what the panel would have seen of it, and the simulated panels train one classifier per signature and platform. The Add object holds the pipeline, and five method calls run it: annotate, cluster, simulate, tune, then publish or export.

The blocks on this page are templates: they run on a cohort of your own in the layout below. The worked example runs them on forty genomes with their printed output; the exposure table's columns are in the template file.

Check the requirements

Step 1: build the object

import pandas as pd, sigma2 as sg2

genomes = sg2.Samples(
    "genomes/",
    ref_genome="hg19",
)          # a folder of VCFs, or a MAF
exposures = pd.read_csv(
    "exposures.csv",
    index_col=0,
)         # samples × catalog signatures
a = sg2.Add(
    genomes,
    exposures=exposures,
    tumor_type="breast",
    seq_type=["msk", "seqcap"],
    soi=["SBS3", "APOBEC"],
)

Add(...) prints an overview on construction: samples, median SNVs, each signature's zero fraction and modality, and the cells that will be trained.

Step 2: annotate the labels

# a Gaussian mixture on the exposures, BIC picks 2 to 5 components
a.annotate("gauss")
a.annotate("gauss", n_components=3)   # or fix the count

Add.annotate is called with its defaults. All parameters of Add.annotate

The exposure distribution of a signature is cut into bands, and each sample gets a bin: 0 is negative, the others are band percentiles (P60, P100, …) and are positive. The output table says how many positives, negatives and bins each signature got.

min_exposure (a relative floor) and below="exclude" are hand-set alternatives; the curation axes (Curate training labels) move the boundary by rule. a.label_report() tabulates positives, negatives, exact zeros and prevalence per signature, and warns when a negative class is too small to certify the operating point.

Step 3: cluster the reference spectra

a.cluster()                                       # elbow heuristic per band
a.cluster(method="distance", distance_threshold=0.3)

Add.cluster is called with its defaults. All parameters of Add.cluster

Within each signature and band, samples are clustered on their exposure fractions, and the cluster centroids become the spectra the likelihood features are computed against. "elbow" picks the count per band and does not grow with cohort size; "distance" holds cluster tightness fixed and grows with the exposure space the cohort covers. Recluster after any change to the labels.

Step 4: simulate the panels

a.simulate(scaling_factor=2.0, include_wgs=True)
a.simulate(do=["msk"])                             # one platform

Each genome's mutations inside a platform's capture regions, topped up by resampling the genome's own spectrum (addition_source to draw from something else) by the scaling_factor, give one simulated panel per genome per platform. include_wgs=True also keeps the genome itself as a wgs cell. The output table gives the median and range of simulated SNVs per platform and the share below the floor; a.sparse_report() breaks that down by bin.

The five bundled platforms (msk, fo, op, seqcap, tcga_mc3) are hg19. Another design is a SeqType with your BED: sg2.SeqType("my_panel.bed", name="msk") reuses a vocabulary name, because custom names are not accepted. Exomes need no top-up. The default factor is SigMA's breast value; sigma2.tools.simulator.calibrate_scaling_factor(base, reference) fits it against a real cohort's burdens on the panel.

Step 5: tune the classifiers

a.tune(
    estimator="monologit",
    cross_fit=True,
    cv_repeats=20,
    cutoff_rule="median_repeat",
    output_dir="results",
)

One learner per signature and platform, in nested cross-validation (5 outer folds, 3 inner), repeated cv_repeats times with new splits. The released cutoffs are set on the out-of-fold scores; median_repeat takes the median of the per-repeat cutoffs. cross_fit=True fits the released model on folds rather than on everything.

The package defaults (estimator="gbm", one repeat, mean_score, no cross-fit) reproduce older runs; the call above is the recipe most shipped cells were trained with and the one to start from. A model's card records its own: config["cross_fit"], config["cv_repeats"] and config["cutoff_rule"], and type(model.estimator).__name__ names its learner. In this release about three quarters of the released panel and exome cells carry this recipe; the cards of the rest record none of the three fields. monologit is a sign-constrained logistic baseline with monotone boosted trees; hist is faster and routes missing values itself; gbm is SigMA's learner. Twenty repeats of a few candidates take tens of minutes on a 16-core server; five repeats a few minutes; memory about 1 GB per candidate. mxw sets the workers.

The run writes a result tree (Files and formats) and prints a summary table per cell. It also prints the gate's warnings as it goes: a cell whose negatives cannot certify the operating point is reported before the training ends.

Training at another floor: snv_cutoff= trains at a fixed floor; snv_constraint={0.10: 0.75} searches the lowest floor meeting the constraints, at a full training pass per candidate. A model trained at a non-default floor records it, and run() applies it. No shipped model used the search.

Step 6: publish or export

a.publish(result_dir="results")
a.export("my_store")

Add.export is called with its defaults. All parameters of Add.export

After a single cross-validation run, what tune() does by default (cv_repeats=1), the gate withholds a cell whose likelihood ratio at the 10% cutoff is below 5, whose ratio's lower bound is at most 1, or that has fewer than ten positives. Under repeated cross-validation the gate withholds a cell only when its false-positive claim is unsupported (fewer than 30 distinct negatives, an unstable operating point, a realized false-positive rate over the target by more than one negative) or its calls carry no information; otherwise it releases, graded weak when there are fewer than ten positives or the likelihood ratio is below the bar. Release status has the rule; the publish() docstring states it as the code acts. Without result_dir nothing can be judged, so everything publishes to the release namespace with a warning.

The cells trained on PCAWG genomes passed two further bars outside publish(): no single input of the classifier beat it on held-out genomes, and, for SBS3, it matched SigMA v1 where SigMA v1 has a classifier. The package's changelog records both; no test in the package checks them.

Check the result

m = sg2.model("breast/SBS3/msk", root_dir="my_store")
m.config["evidence"], m.config["gate_reason"], m.config["snv_cutoff"], m.config["catalog"]
out = sg2.run(
    input="new_panels.maf",
    ref_genome="hg19",
    tumor_type="breast",
    seq_type="msk",
    soi="SBS3",
    model=m,
)

The package draws a finished run's detectability: the detection law on the pooled positives, one burden density per assay, the share of positives the assay captured nothing of, and each cell's AUC beside the AUC from SNV count alone (the detection law):

from sigma2 import plot as sp

sp.dtc("cohort/results", tumor_type="breast", output="detectability.svg", n_boot=100)

The figure above is from a run on the 211 public PCAWG breast genomes with this page's recipe at five repeats. training_summary.json in the result tree holds the gate's verdict; cv_scores.tsv the out-of-fold scores. a.export_threshold_table("results") re-scores the run at three floors without retraining, for a floor decision.

Save and resume

a.save("breast.pkl")              # the whole state, any stage
a = sg2.Add.load("breast.pkl")
a.save("breast.json")             # the editor export, for sigscape studio's curation editor

Add.save is called with its defaults. All parameters of Add.save

Add.load is called with its defaults. All parameters of Add.load

A pickle resumes at the stage it was saved. The JSON is for curation in the studio (Curate training labels) and for inject_simulations() plus tune() when the simulations were made elsewhere; it cannot simulate(), because the per-mutation table is not in it.

Optional: run it in sigscape studio

sigscape studio runs Build (annotate and cluster), Simulate and Tune as kept runs on the local runner, with the curation editor between Build and Simulate and the result tree exportable as a zip. In sigscape studio.