2. Detect a signature (SBS3, APOBEC, any other)

sg2.run() scores one signature in a set of tumors. You give it the mutations, the tumor type, the signature and the platform; it finds the trained classifier for that combination, checks that the classifier may be used, and returns a table with one row per tumor. The signature is any of the catalog: Signature 3, APOBEC, SBS8, SBS18, a composite. HRD and mismatch repair have calls of their own:

Get the example file

curl -LO https://raw.githubusercontent.com/parklab/SigMA/master/inst/extdata/examples/test_mutations_50sample.maf

Run the classifier

import sigma2 as sg2

out = sg2.run(
    input="test_mutations_50sample.maf",
    ref_genome="hg19",
    tumor_type="breast",
    seq_type="msk",
    soi="SBS3",
)
out.df

Look at the result

The package draws the cohort as stacked exposures sorted by the signature's share, and any one panel as its spectrum:

from sigma2 import plot as sp

sp.exps(out.df, soi=sg2.resolve_signature("SBS3"), path="exposures.svg", show_sample_names=False)
top = out.df["SBS3_mva"].idxmax()
row = out.df.loc[top]
sp.sbs(
    out.df.loc[[top]].iloc[:, :96],
    title=f"{top}: Signature 3 score {row['SBS3_mva']:.2f}, {row['confidence']}",
    path="top.svg",
)

6. Plot spectra and exposures has both functions in full.

Check what the panel delivered

A call needs the signature's own mutations on the panel, and most panels hold few. exp_SBS3 in the result is the count the decomposition attributed to Signature 3 in each panel, the package's assigned burden; a panel at zero gave the classifier nothing of the signature:

cols = ["total_snvs", "exp_SBS3", "SBS3_mva", "confidence"]
burden = out.df[cols].sort_values("exp_SBS3", ascending=False)
print(burden.describe().loc[["50%", "mean", "max"]])
none = int((out.df["exp_SBS3"] == 0).sum())
print("panels with no Signature 3 mutation attributed:", none, "of", len(out.df))

      total_snvs   exp_SBS3  SBS3_mva
50%         11.0   0.000000  0.141468
mean   18.289474   2.120242  0.321848
max        104.0  21.509469  0.991699
panels with no Signature 3 mutation attributed: 28 of 38

On this cohort, 28 of 38 scored panels gave the classifier no Signature 3 mutation to read, and those are where a true positive is missed. Most calls sit at the top of this list; four loose calls come from panels with no Signature 3 attributed, on the likelihood of the breast Signature 3 reference spectra alone, and deserve a look at the spectrum. A negative call on a panel is weak evidence of absence, so a count of calls is not the signature's prevalence in the cohort. The expected count behind that, the signature's share times the panel's SNVs, is the quantity the package's detection law is fitted on; When a call is possible says what it means for a cohort, and on your own training the package draws the law and your positives' burdens against it:

from sigma2 import plot as sp

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

Check the calls before you report them

• Hypermutated tumors. A signature run does not look for mismatch-repair deficiency or a POLE defect. Run predict_mmrd_ppd() on the same file: a Signature 3 call on a tumor it calls MMRD or POLE, here the loose call on P-0007087, is read off a spectrum another process shaped. A Signature 3 call on a panel dominated by APOBEC mutations deserves the same look.
• How the cell was released. 39 breast, lung and skin cells were released before the release gate, among them breast APOBEC, SBS8 and SBS18 on MSK-IMPACT: their calls carry no evidence grade and no call_evidence column (release status).

Choose the arguments

input, tumor_type and ref_genome are what the table says. Beyond it: every SNV in a VCF counts whatever its FILTER, so filter to PASS before scoring (Prepare your input files); the tumor type sets the signature set, the reference spectra and the classifier, so a wrong one gives a wrong answer without an error, and a tumor type without a signature set in the catalog fails early with an error naming it (the vocabulary lists all 39); a MAF's NCBI_Build is checked against ref_genome, and a matrix input needs none (genome builds). The platform and the signature take more:

Choose the platform with seq_type

The platform the tumors were sequenced with. Pass it whenever you know it:

Names match in any case. seq_type="auto" picks, among the platforms with a classifier for this tumor type and signature, the one whose training cohort's median SNV count is closest to your input's, and reports the choice:

out = sg2.run(
    input="test_mutations_50sample.maf",
    ref_genome="hg19",
    tumor_type="breast",
    seq_type="auto",
    soi="SBS3",
    verbose=1,
)

seq_type="auto" chose "msk": the input's median is 8 SNVs, the msk classifier's training median 9
(others: seqcap 36, tcga_mc3 189, wgs 4310). Burden cannot tell panels of similar size apart;
pass seq_type= when the platform is known.

Panels of similar size (msk, fo, op) have similar burdens, so auto cannot tell them apart. A custom panel has no name in the vocabulary: score it as the closest bundled platform and report it as such.

Choose the signature with soi

The signature of interest: a COSMIC name (SBS3, SBS8, SBS18, …), a composite (APOBEC for SBS2 and SBS13, SBS7 for SBS7a to d, SBS10 for SBS10a to d, SBS17 for SBS17a and b, CLOCK for SBS1 and SBS5), an alias (UV for SBS7, POLE for SBS10), or HRD for the HRD classifier (Detect HRD). Any case. One signature per run; for several, loop over the cells:

from sigma2 import classifiers

calls = {}
for card in classifiers.cards(tumor_type="breast", seq_type="msk"):
    if card.kind != "signature" or card.status == "experimental":
        continue
    out = sg2.run(input="test_mutations_50sample.maf", ref_genome="hg19", tumor_type="breast",
                  seq_type="msk", soi=card.classifier)
    calls[card.classifier] = out.df[[
        f"{card.classifier}_mva",
        "pass_mva",
        "pass_mva_strict",
    ]]

Choose which classifier runs

model="auto" (the default) takes the released classifier for the tumor type, signature and platform. The run names it:

Classifier: tumor_type="breast", soi="SBS3", seq_type="msk" (verified)
Using classifier: breast/SBS3/msk

verified means the stored weights match their recorded checksum. The alternatives:

m = sg2.model("breast/SBS3/msk")
out = sg2.run(
    input="test_mutations_50sample.maf",
    ref_genome="hg19",
    tumor_type="breast",
    seq_type="msk",
    soi="SBS3",
    model=m,
)

model is called with its defaults. All parameters of model

Handle a refusal

No classifier for the cell. The tumor type has no classifier for that signature on that platform. The error lists the platforms that have one and says when SigMA v1 covers the cell:

try:
    sg2.run(
        input="test_mutations_50sample.maf",
        ref_genome="hg19",
        tumor_type="breast",
        seq_type="seqcap_probe",
        soi="SBS3",
    )
except ValueError as e:
    print(e)

No classifier for tumor_type="breast", soi="SBS3", seq_type="seqcap_probe".
  breast SBS3 has one on: fo, msk, op, seqcap, tcga_mc3, wgs.
  SigMA v1 has one for this cell: pass model="legacy".
  sg2.describe() lists every classifier.

classifiers.cards(tumor_type=...) returns the same list as data.

An experimental classifier. The only classifier for the cell was withheld by the release gate:

try:
    sg2.run(
        input="test_mutations_50sample.maf",
        ref_genome="hg19",
        tumor_type="prost",
        seq_type="op",
        soi="SBS3",
    )
except ValueError as e:
    print(e)

the only classifier for tumor_type="prost", soi="SBS3", seq_type="op" is EXPERIMENTAL: withheld by the
release gate (...). Pass force=True to use it anyway; its scores rank samples, but pass_mva and the
confidence labels from it are not supported.

With force=True the run returns <soi>_mva and leaves pass_mva, pass_mva_strict and confidence as NA. The score ranks the samples of your cohort; there is no validated cutoff and no call.

A catalog the classifier was not trained on. Without catalog=, the classifier's own catalog is used. Another catalog produces a warning that the scores are not the validated ones. A catalog that lacks the signature of interest, or that has no signature set for the tumor type (cosmic_v2), is refused. Floor, cutoffs and catalog has the rule.

A tumor type the HRD model cannot score. run(soi="HRD") refuses, with the reason, the tumor types its model saw in training but has no signature set for. Detect HRD lists them.

Optional: run without the classifier

mva=False skips the classifier and returns the features: the likelihoods under the tumor type's reference spectra, the cosine similarities, and the decomposition with the signature's exposure and share. No score, no call. It runs for any signature of the catalog and any tumor type with a signature set, with or without a trained cell, and takes the catalog by name:

feats = sg2.run(
    input="test_mutations_50sample.maf",
    ref_genome="hg19",
    tumor_type="breast",
    seq_type="msk",
    soi="SBS8",
    mva=False,
    catalog="cosmic_v3p2_inhouse",
)
feats.df[["total_snvs", "SBS8_c", "SBS8_l_rat", "exp_SBS8", "rat_SBS8"]].head(3)

The working columns explains each. assign=False leaves the likelihood assignment out as well.

Reference: every argument

Every argument of run(), including the ones the examples above do not pass:

Full signatures and types: Python API.

Read the result object

out.df is the table: the 96 channels, total_snvs, the decomposition features, the score and the calls, one row per scored sample. out.to_lite() returns the columns a reader needs as a DataFrame. out.snv_floor, out.catalog_name, out.seq_type, out.soi and out.model record what ran. Read the output describes every column.