Adapt SigMA2 to your own data

Know the pipeline

A run is a fixed sequence of steps on one 96-channel matrix. Samples reads the mutations into the matrix and keeps the indels beside it; the decomposition and the likelihood step compute, per tumor, the exposures, the cosine similarities and the likelihoods of the tumor type's reference clusters; the classifier step scores those features with the cell's model and applies its cutoffs. HRD and MMRD/POLE add their own feature steps in front of their own models. How a classifier scores a panel describes the features, and what a call means what each model was trained on.

Everything a run used is on the object it returns:

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

import sigma2 as sg2

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

(GBMModel(breast/SBS3/msk), 'cosmic_v3p2_inhouse', 5, 'msk', SBS3)

Build a matrix in your own pipeline

When your mutations do not come as a MAF or VCF, pass the 96-channel matrix as a DataFrame: one row per tumor, one column per channel, named by trinucleotide context in COSMIC spelling (A[C>A]A) or the package's (a[c>a]a), in any order, as counts. A matrix carries no indels, so MMRD/POLE and HRD need the indel counts as columns beside it; sg2.indel_features computes them from a table of indel records with the columns Chromosome, Position, REF, ALT and Sample, in VCF style (the anchor base included) or MAF style (- for the empty allele, a deletion placed at its first deleted base and an insertion at the base before it). On the example file, the matrix route gives the same indel counts and the same MMRD/POLE calls as passing the MAF itself:

import pandas as pd

maf = pd.read_csv(
    "test_mutations_50sample.maf",
    sep="\t",
    comment="#",
    dtype=str,
    keep_default_na=False,
)
indels = maf.loc[maf["Variant_Type"].isin(["INS", "DEL"])].rename(columns={
    "Tumor_Sample_Barcode": "Sample",
    "Start_Position": "Position",
    "Reference_Allele": "REF",
    "Tumor_Seq_Allele2": "ALT",
})[["Chromosome", "Position", "REF", "ALT", "Sample"]]
features = sg2.indel_features(
    indels=indels,
    ref_genome="hg19",
    samples=maf["Tumor_Sample_Barcode"].unique().tolist(),
)
spectra = sg2.Samples("test_mutations_50sample.maf", ref_genome="hg19").df.iloc[:, :96]
mm = sg2.predict_mmrd_ppd(
    input=spectra.join(features),
    tumor_type="breast",
    seq_type="msk",
)
mm.df["MMRD_category"].value_counts()

MMRD_category
MMRP    37
MMRD     1
Name: count, dtype: int64

The indel table's Sample names must be the matrix's row names. An indel whose sample matches no row is counted for no tumor, without a message: that tumor is scored as if it had no indels, and its MMRD/POLE and HRD calls change (renaming the samples in the indel table above turns the MMRD tumor into an MMRP call). An empty set means every name matched:

set(indels["Sample"]) - set(spectra.index)

set()

The same table serves HRD: with its mmej and nhej columns, run(soi="HRD") on the matrix gives the same calls as on the MAF. The repeat loci SigMA2 ships are hg19; on another build nmsi_ins and nmsi_del are left empty unless repeat-locus BED files of that build are passed with beds=, and the MMRD/POLE classifier then scores without them, printing nothing at the default verbosity. Split multi-allelic records before counting: indel_features skips an indel whose ALT lists two alleles, without a message.

Check what the reader accepts

The reader refuses what it cannot score and, in most cases, names the reason. Observed on the current release:

Filter a VCF to PASS before scoring. 1. Prepare your input files lists the columns of each format and the builds each input accepts.

Run it in a pipeline

The Python call returns the table and writes nothing; a batch step writes what it keeps:

out.df.to_csv("sbs3_full.tsv", sep="\t")
out.to_lite().to_csv("sbs3_lite.tsv", sep="\t")
out.df.shape, out.to_lite().shape

((38, 186), (38, 14))

An HRD run is also a command, sigma2 hrd score, for a workflow manager that calls programs (Command line). Pass mxw="max" on a cohort to use every core.

Record what a result depends on

A result depends on the package release, the classifier and its catalog, the floor, the reference genome build and, for indel features, the repeat loci. Keep sg2.__version__, out.model.name, out.catalog_name and out.snv_floor with every table; the model's card (out.model.config) holds its cutoffs and its training record. Pin the release when you install it, with pip install "sigma2==0.0.16a20".

Know what is public

• Public: the functions on the Python API pages (scoring, models and the store, HRD, training, inputs, the SigMA v1 tools) and the columns of 5. Read the output.
• Internal: sigma2.api, sigma2.steps, sigma2.core.Pipeline, sigma2.tools, sigma2.utils and the private members of the model classes. sigma2.__all__ lists api, steps and tools, so they import like the public names; the reference's Advanced branch lists them as the package has them; nothing promises their names or behavior from one release to the next, and a pipeline that imports them can break on an upgrade.

Extend what the package ships

Reference: errors and their causes