Worked example: fifty breast panels
Run every classifier that applies to one public file of fifty breast panels, step by step.
Before you begin, you need the package installed, the hg19 genome, the repeat-locus files; about 10 minutes.
One public file, every classifier that applies to it: Signature 3 and APOBEC, the new SBS3 classifier beside SigMA v1's, HRD with and without copy-number segments, and mismatch-repair deficiency. Steps 5 and 6 need the repeat-locus files (sigma2 data fetch --all, Install).
1. Get the data
Fifty MSK-IMPACT breast panels, SigMA (R)'s example MAF, and their copy-number segments from this site:
curl -LO https://raw.githubusercontent.com/parklab/SigMA/master/inst/extdata/examples/test_mutations_50sample.maf
curl -LO https://sigscape.org/sigma2/examples/msk_breast_50_segments.seg
import sigma2 as sg2
MAF = "test_mutations_50sample.maf"
MSK-IMPACT is seq_type="msk", the tumors are tumor_type="breast", and the coordinates are hg19.
Note Thirty-nine of these fifty panels were training panels of the
mskHRD and MMRD/POLE models, so the HRD and MMRD outputs below show how to run and read the models, not how they perform on tumors they have not seen.
2. Score Signature 3
sbs3 = sg2.run(
input=MAF,
ref_genome="hg19",
tumor_type="breast",
seq_type="msk",
soi="SBS3",
verbose=1,
)
| Parameter | Description | Default |
|---|---|---|
| inputpath, list of paths, folder or DataFrame | The mutations: a MAF path, a VCF path, a list of VCF paths, a folder of VCFs, a CSV, or a DataFrame of SBS-96 counts. | Required |
| ref_genomestr or path | The build of the coordinates (hg19, hg38, GRCh37, GRCh38) or a FASTA path; unused for a matrix input. | 'hg19' |
| tumor_typestr | A tumor type of the vocabulary, in any case (the package lists 39 names; msi and pole among them are signature sets, not tumor types); it sets the signature set, the reference spectra and the classifier. | Required |
| seq_typestr | The platform (msk, fo, op, seqcap, tcga_mc3, wgs); "auto" chooses by SNV burden and reports the choice. | 'auto' |
| soistr or CompositeSignature | The signature of interest: a COSMIC name, a composite preset (APOBEC, SBS7, SBS10, SBS17, CLOCK), an alias (UV, POLE), or HRD. | Required |
| verboseint | 0 prints errors and some warnings (a catalog other than the classifier's, a readjust without a table); 1 also prints the classifier, the sample counts and the other warnings, such as a legacy classifier's features passed as missing; 2 also lists the samples below the floor. | 0 |
Classifier: tumor_type="breast", soi="SBS3", seq_type="msk" (verified)
Input: 50 samples → 38 after filtering (12 dropped)
Using classifier: breast/SBS3/msk
Twelve panels have fewer than five SNVs and are not scored; verbose=2 lists them. The result for the other 38:
sbs3.df[[
"total_snvs",
"SBS3_mva",
"pass_mva",
"pass_mva_strict",
"confidence",
"call_evidence",
]].sort_values("SBS3_mva", ascending=False).head(3)
total_snvs SBS3_mva pass_mva pass_mva_strict confidence call_evidence
P-0006809-T01-IM5 33 0.991699 True True high informative
P-0010938-T01-IM5 41 0.991418 True True high informative
P-0008488-T02-IM5 26 0.962242 True True high informative
Three high-confidence calls at the top; the loose calls come lower in the list, P-0007087 at rank six with a score of 0.78. call_evidence is informative on every row: this is a released cell. The breast SBS3 classifier on msk was trained on strong SBS3 only, so a sample with weak SBS3 is expected to be negative here; its card records this (sg2.model("breast/SBS3/msk").config["label_curation"]).
Look at the calls
The spectra behind the three highest scores, drawn by the package, each titled with its score and confidence. Each spreads its 26 to 41 SNVs over 22 to 30 of the 96 channels and five or six of the substitution classes, the way Signature 3 spreads its own:
from sigma2 import plot as sp
top = sbs3.df.sort_values("SBS3_mva", ascending=False).head(3)
for sid, row in top.iterrows():
sp.sbs(
sbs3.df.loc[[sid]].iloc[:, :96],
title=f"{sid}: Signature 3 score {row['SBS3_mva']:.2f}, {row['confidence']}",
path=f"{sid}.svg",
)
3. Run the same panels through SigMA v1
v1 = sg2.run(
input=MAF,
ref_genome="hg19",
tumor_type="breast",
seq_type="msk",
soi="SBS3",
model="legacy",
)
both = sbs3.df[[
"total_snvs",
"SBS3_mva",
"pass_mva",
]].join(v1.df[["SBS3_mva", "pass_mva"]], rsuffix="_v1")
both.sort_values("SBS3_mva", ascending=False).head(3)
total_snvs SBS3_mva pass_mva SBS3_mva_v1 pass_mva_v1
P-0006809-T01-IM5 33 0.991699 True 0.956447 True
P-0010938-T01-IM5 41 0.991418 True 0.937323 True
P-0008488-T02-IM5 26 0.962242 True 0.935358 True
The legacy run used SigMA v1's catalog (cosmic_v2_inhouse) by default; Reproduce SigMA v1 lists the settings that reproduce R's result in full. The two classifiers agree at the top. Lower in the list they part: P-0007087 is a loose call from the new classifier and none from the old, and step 6 calls it MMRD.
4. Score APOBEC
apobec = sg2.run(
input=MAF,
ref_genome="hg19",
tumor_type="breast",
seq_type="msk",
soi="APOBEC",
)
apobec.df[[
"total_snvs",
"APOBEC_mva",
"pass_mva",
"pass_mva_strict",
"confidence",
]].sort_values("APOBEC_mva", ascending=False).head(3)
total_snvs APOBEC_mva pass_mva pass_mva_strict confidence
P-0001939-T02-IM5 13 0.990404 True True high
P-0009841-T01-IM5 13 0.990404 True True high
P-0004246-T01-IM5 57 0.990264 True True high
APOBEC is the SBS2 and SBS13 composite; the columns are named by it. Two panels of 13 SNVs score as high as the panel with 57. The panel with 57 SNVs was negative for Signature 3; P-0007087, the loose Signature 3 call, is negative for APOBEC.
Look at the APOBEC calls
The cohort as stacked exposures sorted by SBS13, one of the composite's two members, and the spectrum of the highest APOBEC score; APOBEC's signatures put C>T and C>G at a TpC context, and the top panel has 7 of its 10 C>T and C>G there:
sp.exps(
apobec.df,
soi=sg2.resolve_signature("SBS13"),
path="apobec_exposures.svg",
show_sample_names=False,
)
sid = apobec.df["APOBEC_mva"].idxmax()
row = apobec.df.loc[sid]
sp.sbs(
apobec.df.loc[[sid]].iloc[:, :96],
title=f"{sid}: APOBEC score {row['APOBEC_mva']:.2f}, {row['confidence']}",
path="apobec_top.svg",
)
5. Score HRD, without and with copy number
Without segments, the HRD model runs its nogis variant on Signature 3 and the indel features. The indel features and the MMRD/POLE pre-step need the repeat-locus files:
hrd0 = sg2.run(
input=MAF,
ref_genome="hg19",
tumor_type="breast",
seq_type="msk",
soi="HRD",
verbose=1,
)
hrd0.df[[
"total_snvs",
"sig3_score",
"mmej",
"nhej",
"hrd_variant",
"hrd_score",
"hrd_call",
"MMRD_category",
]].head(3)
HRD: no copy-number input (segments= or gis=); scoring with the model trained without copy number.
An HRD run keeps every sample: the twelve below the floor have NA Signature 3 features and NA MMRD_category, and their HRP is not an assessment (Detect HRD). The pre-step ran, because a MAF has indels: MMRD_category is MMRP for the rest, and an MMRD tumor carries that word as its hrd_call.
With the panels' copy-number segments, an IMPACT-style .seg file (ID, chrom, loc.start, loc.end, num.mark, seg.mean), the model's tcn variant reads the fraction of genome altered and a log-ratio LST on its gene grid:
hrd = sg2.run(
input=MAF,
ref_genome="hg19",
tumor_type="breast",
seq_type="msk",
soi="HRD",
segments="msk_breast_50_segments.seg",
)
print(hrd.df[[
"total_snvs",
"sig3_score",
"fga",
"lst_lr",
"mmej",
"nhej",
"hrd_variant",
"hrd_score",
"hrd_call",
]].sort_values("hrd_score", ascending=False).head(1))
hrd.df["hrd_call"].value_counts()
total_snvs sig3_score fga lst_lr mmej nhej hrd_variant hrd_score hrd_call
P-0006809-T01-IM5 33 0.956447 0.469237 14.0 0.0 1.0 tcn 0.988345 HRD
hrd_call
HRP 41
HRD-low 4
HRD 4
MMRD 1
Name: count, dtype: int64
HRD-low means the pan-cancer and the tumor-type model disagree, one above 0.5 and the other not; an orthogonal assay settles it. Allele-specific segments (ASCAT, FACETS) give the gis variant instead; Detect HRD describes the three variants and HRD on exomes runs the gis one.
6. Classify mismatch repair and POLE
mm = sg2.predict_mmrd_ppd(
input=MAF,
ref_genome="hg19",
tumor_type="breast",
seq_type="msk",
)
print(mm.df["MMRD_category"].value_counts())
mm.df.loc[
mm.df["MMRD_category"] == "MMRD",
["total_snvs", "nmsi_ins", "nmsi_del", "MMRD_prob", "POLE_prob", "MMRP_prob"],
]
| Parameter | Description | Default |
|---|---|---|
| inputpath, list of paths, folder or DataFrame | The mutations with their indels: a MAF or VCF path, or a matrix with nins, ndel, nmsi_ins and nmsi_del columns. | Required |
| ref_genomestr or path | The build of the coordinates (hg19, hg38, GRCh37, GRCh38) or a FASTA path. | 'hg19' |
| tumor_typestr | A tumor type of the vocabulary (msi and pole are signature sets, not tumor types); the point-mutation features need its signature set. | Required |
| seq_typestr | The platform: a model exists for msk, tcga_mc3, seqcap, fo and op. | Required |
MMRD_category
MMRP 37
MMRD 1
Name: count, dtype: int64
total_snvs nmsi_ins nmsi_del MMRD_prob POLE_prob MMRP_prob
P-0007087-T01-IM5 43 1.0 8.0 0.999076 0.00005 0.000874
The panel with the loose Signature 3 call is called MMRD, and the HRD run set it aside for that reason (3. Detect HRD). model="legacy" runs SigMA v1's classifier instead.
Look at the MMRD panel
The spectrum of the one MMRD call, drawn by the package: 23 of its 43 SNVs are T>C, and the MMRD run gives signatures 21 and 26, two of the signatures associated with defective mismatch repair, its largest exposures (4. Detect MMRD or POLE shows them). The panel's loose Signature 3 call came from this spectrum:
sid = mm.df.index[mm.df["MMRD_category"] == "MMRD"][0]
sp.sbs(
mm.df.loc[[sid]].iloc[:, :96],
title=f"{sid}: MMRD, probability {mm.df.loc[sid, 'MMRD_prob']:.2f}",
path="mmrd.svg",
)
7. Write one table
import pandas as pd
table = pd.concat({
"SBS3": sbs3.df["pass_mva"], "SBS3_v1": v1.df["pass_mva"], "APOBEC": apobec.df["pass_mva"],
"HRD": hrd.df["hrd_call"], "MMRD": mm.df["MMRD_category"],
}, axis=1)
table.to_csv("breast_panels_calls.tsv", sep="\t")
Record with it the release (sg2.__version__), the platform, the tumor type, and that the floors and catalogs were the classifiers' own. Every call here is at the released operating point of a released cell, and breast APOBEC on MSK-IMPACT was released before the release gate, so its calls carry no evidence grade. Before reporting, set aside the Signature 3 and APOBEC calls on the panel called MMRD (P-0007087), whose spectrum another process shaped, and the HRD calls of the twelve panels below the floor; what it means for a patient is the study's question.
Take it to your own data
- Your own panels: Prepare your input files for the MAF and VCF rules, then the same calls.
- A cohort sequenced deeper or shallower than MSK-IMPACT 2017: read readjust before switching it on.
- The same tumors in sigscape studio, with the model card beside the table and every run kept: In sigscape studio.