Migrate from SigMA v1
Run a SigMA (R) analysis in SigMA2, translate its arguments, columns and functions, and see what behaves differently.
Before you begin, you need the package installed, the hg19 genome; about 5 minutes.
Run your first SigMA2 analysis
SigMA v1's run() took a file, a platform and a tumor type and wrote a CSV beside the input; SigMA2's takes the same and returns the table. With model="legacy", SigMA2 runs SigMA v1's own classifier on SigMA v1's catalog. On SigMA v1's example file:
curl -LO https://raw.githubusercontent.com/parklab/SigMA/master/inst/extdata/examples/test_mutations_50sample.maf
import sigma2 as sg2
v1 = sg2.run(
input="test_mutations_50sample.maf",
ref_genome="hg19",
tumor_type="breast",
seq_type="msk",
soi="SBS3",
model="legacy",
)
v1.df[["total_snvs", "SBS3_mva", "pass_mva", "pass_mva_strict"]].head(3)
| 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 |
| modelstr or GBMModel | Which classifier runs: "auto" for the released cell, "legacy" for SigMA v1's, or a model loaded with sg2.model(). | 'auto' |
total_snvs SBS3_mva pass_mva pass_mva_strict
P-0008709-T01-IM5 18 0.892601 True True
P-0010938-T01-IM5 41 0.937323 True True
P-0010560-T01-IM5 5 0.086188 False False
The package's tests run this call on SigMA v1's MSK-IMPACT breast example and compare it with the output R wrote: the short table and every Signature 3 column agree to within 1e-13, and every call agrees. The likelihood ratios and cosine differences of the other signatures do not: they differ from R's by up to 0.19, and by up to 0.41 in the MSI columns of a run with msi=True. Leave model out to run SigMA2's own classifier for the cell instead, a different model on another catalog whose calls can differ (Reproduce SigMA v1).
Caution The agreement holds for the classifiers that read only features SigMA2's port computes: 30 of the 86 SigMA v1 classifiers that run, MSK-IMPACT breast among them. The other 56 read a feature that the port passes to the classifier as missing while R computes it: SigMA v1's Signature 3 ratio
rat_sig3, or an APOBEC, clock or oxidative cluster likelihood. Their scores and calls can differ from R's. On R's output for 1,018 TCGA breast exomes, scored by thetcga_mc3breast classifier with one APOBEC feature missing, 8 loose and 3 strict calls differ; with the feature supplied, the scores agree with R's to within 1e-14.verbose=1prints the features a run passed as missing, and Reproduce SigMA v1 lists them for every classifier.
Translate the arguments
SigMA v1 run() | SigMA2 run() | What changes |
|---|---|---|
genome_file, input_df | input= | a MAF or VCF path, a folder of VCFs, a matrix CSV or a DataFrame |
data="msk" | seq_type="msk" | SigMA v1's classifiers keep their platform names (msk, seqcap, seqcap_probe, tcga_mc3, wgs, wgs_pancan); SigMA2's own are on fo, msk, op, seqcap, tcga_mc3 and wgs. SigMA v1's other settings (msk_on_target, ampliseq, the Illumina and Invitae panels) run only with mva=False. R's default was "msk"; SigMA2's is "auto", chosen by mutation burden, so name the platform |
tumor_type="breast" | tumor_type="breast" | SigMA v1's names except other, which has no classifier in SigMA2; R's default was "breast", SigMA2 requires the argument |
catalog_name= (required) | catalog= (optional) | left out, the classifier's own: cosmic_v2_inhouse for SigMA v1's, cosmic_v3p2_inhouse for SigMA2's |
do_mva=T | mva=True | mva=False returns the decomposition and the likelihoods without a score. Its only call column is then pass_ml, the Signature 3 likelihood above 0.5: what SigMA v1 called with do_mva=F for a tumor type without a classifier, and what SigMA2 does not report as a call |
do_assign=T | assign=True | the calls at the cutoffs |
check_msi=T | msi=True | the MSI and POLE likelihoods on a signature run; the MMRD/POLE call itself is predict_mmrd_ppd() |
weight_cf=T | wcf=True | SigMA v1's cluster-size prior, with model="legacy" on msk |
lite_format=T | out.to_lite() | the short table, from the returned object |
snv_cutoff=5 | snv_ctf=5 | the default differs: SigMA2 applies the tumor type's floor, 4 for prostate and 3 for osteosarcoma and pancreatic neuroendocrine tumors |
readjust=T | readjust=True | SigMA2's own classifiers readjust from their cutoff tables; with model="legacy" SigMA2 has no port of SigMA v1's readjusted cutoffs, warns, and applies the published ones. SigMA v1 also counted the tumors below the floor in the cohort's median and could switch the platform first; SigMA2's median counts the tumors above the floor only |
return_df=T, output_file= | the returned object | out.df.to_csv(...) writes the table |
custom=T | model= a model of your own store | train it with sg2.Add (Train a classifier) |
add_sig3=T | none | a signature is scored only against the tumor type's own signature set and reference spectra |
norm96= | none | no public weighting for a custom panel's trinucleotide content; score a custom panel as the closest bundled platform. run() still accepts norm96 without a message, and the scores change with it: leave it out |
The signature of interest is new: soi="SBS3" scores Signature 3, and any other signature with a classifier can be named (2. Detect a signature).
Translate the output
SigMA v1 names a column by Signature_3 or sig3; SigMA2 by the COSMIC name. The short table of SigMA v1's lite_format and of SigMA2's to_lite():
| SigMA v1 | SigMA2 |
|---|---|
tumor | the row index |
total_snvs | total_snvs |
Signature_3_ml, Signature_8_ml, Signature_17_ml, Signature_clock_ml, Signature_APOBEC_ml | SBS3_ml, SBS8_ml, SBS17_ml, clock_ml, APOBEC_ml |
Signature_3_c | SBS3_c |
exp_sig3 | exp_SBS3 |
Signature_3_l_rat | SBS3_l_rat |
Signature_3_mva | SBS3_mva |
pass_mva, pass_mva_strict, pass_ml, sigs_all, exps_all | the same names; sigs_all lists SigMA2's signature names (SBS3.SBS8.SBS18 for SigMA v1's Signature_3.Signature_8.Signature_18) |
categ (Signature_3_hc, Signature_3_lc, Signature_clock, …) | categ (SBS3_hc, SBS3_lc, clock, …) |
| none | rat_SBS3 (SigMA v1's full table has it as rat_sig3) and confidence (high, low, negative) |
Signature_msi_ml, Signature_pole_ml, categ values Signature_msi, Signature_pole | none in to_lite(); see below |
list(v1.to_lite().columns)
['total_snvs', 'APOBEC_ml', 'SBS17_ml', 'SBS3_ml', 'SBS8_ml', 'clock_ml', 'SBS3_c', 'exp_SBS3', 'rat_SBS3', 'SBS3_l_rat', 'SBS3_mva', 'pass_mva', 'pass_mva_strict', 'confidence', 'sigs_all', 'exps_all', 'pass_ml', 'categ']
SigMA v1's get_sig_exps unpacked sigs_all and exps_all into one exposure column per signature. SigMA2 has no function for it; pandas does it:
import pandas as pd
rows = [
dict(zip(sigs.split("."), map(float, exps.split("_"))))
for sigs, exps in zip(v1.df["sigs_all"], v1.df["exps_all"])
]
exposures = pd.DataFrame(rows, index=v1.df.index).fillna(0.0).add_prefix("exp_")
exposures.iloc[:2, :4]
exp_SBS3 exp_SBS8 exp_SBS18 exp_SBS30
P-0008709-T01-IM5 6.141477 0.160019 2.994815 5.134828
P-0010938-T01-IM5 21.900193 0.000000 0.000000 5.115141
to_lite() does not carry SigMA v1's MSI and POLE likelihood sums, and its categ has no MSI or POLE category; a tumor SigMA v1 categorizes Signature_msi reads by its best other group. The likelihoods are in the full table of a run with msi=True, and their sum is SigMA v1's Signature_msi_ml:
m = sg2.run(
input="test_mutations_50sample.maf",
ref_genome="hg19",
tumor_type="breast",
seq_type="msk",
soi="SBS3",
model="legacy",
msi=True,
)
import pandas as pd
msi_ml = m.df.filter(regex=r"^msi_.*_ml_msi$").sum(axis=1).rename("msi_ml")
pd.concat([msi_ml, m.to_lite()["categ"]], axis=1).sort_values("msi_ml", ascending=False).head(2)
| 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 |
| modelstr or GBMModel | Which classifier runs: "auto" for the released cell, "legacy" for SigMA v1's, or a model loaded with sg2.model(). | 'auto' |
| msibool | SigMA v1's MSI features on a signature run; for MMRD/POLE calls use predict_mmrd_ppd. | False |
msi_ml categ
P-0007087-T01-IM5 1.000000 clock
P-0007288-T01-IM5 0.995339 clock
Find each SigMA v1 function
Every function the R package exports, and what carries it in SigMA2. Reproduced means the package's tests compare the SigMA2 function with R on the same input and they agree; ported means the function was carried over without such a comparison in this release's tests.
| SigMA v1 | SigMA2 | Status |
|---|---|---|
run() | sg2.run(model="legacy") | reproduced, the short table and the Signature 3 columns, for the 30 classifiers that read every feature; for the other 56, ported, and their scores can differ from R's (the caution above) |
lite_df() | Output.to_lite() | ported, under SigMA2's names, without the MSI and POLE columns above |
predict_mmrd() | sg2.predict_mmrd_ppd(model="legacy") | ported for msk and tcga_mc3; the tests match SigMA v1's published MMRD/POLE calls on TCGA exomes |
get_threshold(), sen_fpr() | sg2.get_threshold(), sg2.sen_fpr() | reproduced |
llh_max_characteristics_wgs_data() | sg2.match_training_genomes() | reproduced |
llh_max_characteristics() | sg2.cluster_characteristics() | reproduced; the comparison needs R with its gbm and nnls packages |
quick_simulation() | sg2.quick_simulation() | reproduced in its deterministic steps, a comparison that needs R with gbm and nnls; the random placement of mutations differs (seed=), and maf_percent is refused |
match_to_catalog(), calc_llh(), decompose() | sg2.run(mva=False) | replaced: one call returns the likelihoods, cosine similarities and exposures |
get_sig_exps() | none | the exposures stay packed in sigs_all and exps_all; Translate the output unpacks them |
make_matrix(), conv_snv_matrix_to_df() | sg2.Samples() | replaced: the 96-channel matrix from a MAF, VCFs or a matrix file |
mmej_nhej_indels_from_vcfs(), calculate_mh(), mh_from_vcf() | sg2.indel_features() | replaced; the microhomology rule follows calculate_mh and the tests check it |
overlap_repeats(), run_overlap_repeats() | sg2.indel_features() | replaced: the microsatellite counts against the bundled repeat loci |
find_data_setting() | run(seq_type="auto") | replaced: the platform is chosen from the cohort's mutation burden, as find_data_setting chose the setting whose SNV counts lie closest |
list_data_options(), list_tumor_types(), has_model() | sg2.describe(), sg2.models(), sigma2.classifiers.cards() | replaced: the classifiers as data |
tune_new_gbm(), tune_gbm_model(), add_gbm_model(), remove_gbm_model() | sg2.Add, Add.tune(), Add.export(), sg2.delete_model() | replaced by another method: whole genomes simulated as panels, nested cross-validation and a release gate (Train a classifier) |
simulate_from_wgs(), subsample(), subsample_id(), get_base_matrix(), get_inter_tt_suppl() | inside sg2.quick_simulation() and Add.simulate() | internal |
plot_summary(), plot_detailed(), plot_tribase_dist() | sigma2.plot.exps(), sigma2.plot.sbs() | replaced: the stacked exposures and the spectrum; no summary page (6. Plot spectra and exposures) |
get_catalog(), get_trinuc_norm(), get_trinuc_counts(), get_tumor_types_for_signature(), get_wgs_data_for_llh(), get_Sig3_llh_from_wgs_data(), get_exps_mmrd(), info_stat(), predict_prob(), predict_mva() | none public | the reference data and the classifier internals are not public functions; sg2.get() reads the reference tables |
Choose a cutoff on panels like yours
SigMA v1's tools around a shipped classifier run on SigMA2's classifiers too. A quick simulation brings the package's whole genomes of known Signature 3 status to your cohort's depth; scored by the Signature 3 classifier for your platform, it gives the cutoff at a false-positive rate of your choice:
sim = sg2.quick_simulation(
input="test_mutations_50sample.maf",
ref_genome="hg19",
tumor_type="breast",
seq_type="msk",
seed=1,
)
scored = sg2.run(
input=sim.set_index("tumor").iloc[:, :96],
tumor_type="breast",
seq_type="msk",
soi="SBS3",
).df.join(sim.set_index("tumor")["is_sig3"])
cut = sg2.get_threshold(
scored,
limits=[0.1, 0.01],
var="SBS3_mva",
signal="is_sig3",
cut_var="fpr",
)
cut[["limit", "cutoff"]]
| Parameter | Description | Default |
|---|---|---|
| inputpath, list of paths or DataFrame | The cohort to imitate: a MAF or VCF path, a list of VCF paths, or a 96-channel matrix. | Required |
| ref_genomestr | The build of the coordinates (hg19, hg38) or a FASTA path; unused for a matrix input. | 'hg19' |
| tumor_typestr | The tumor type whose built-in whole genomes are brought to the cohort's depth. | Required |
| seq_typestr | The platform the cohort was sequenced with; the whole genomes are brought to its depth. | Required |
| seedint | Makes the random subsampling reproducible. | None |
| Parameter | Description | Default |
|---|---|---|
| limitslist of float | The rates to read the cutoff at, one row of the result per limit. | Required |
| varstr | The score column. | 'prob' |
| signalstr | The column of known labels: 1 or True for a positive. | 'is_sig3' |
| cut_varstr | "fpr" or "fdr" for an upper bound on that rate, "sen" for a lower bound on sensitivity. | 'fpr' |
limit cutoff
0 0.10 0.482495
1 0.01 0.950058
Of the 471 simulated panels above the floor, 296 are negatives, and the 1% cutoff falls between the second- and third-highest of their scores.
The run above takes input, tumor_type, seq_type and soi as in Run your first SigMA2 analysis. Only Signature 3 can be simulated, as in SigMA v1.
Find the genomes and clusters a tumor resembles
SigMA v1's llh_max_characteristics_wgs_data matched each tumor to the single training genome it resembles most, and llh_max_characteristics read the signatures of its most likely cluster:
spectra = sg2.Samples("test_mutations_50sample.maf", ref_genome="hg19").df.iloc[:, :96]
match = sg2.match_training_genomes(spectra, tumor_type="breast")
cc = sg2.cluster_characteristics(v1.df, tumor_type="breast")
match[["max_ml_sample", "max_ml_sample_is_sig3", "sig3_genomes_ml"]].join(cc["max_clust"]).head(2)
| Parameter | Description | Default |
|---|---|---|
| ref_genomepath or str | The build of the coordinates (hg19, hg38) or a FASTA path. | 'hg19' |
| Parameter | Description | Default |
|---|---|---|
| tumor_typestr | The tumor type whose training genomes are matched against. | Required |
| Parameter | Description | Default |
|---|---|---|
| tumor_typestr | The tumor type whose clusters and signatures are used. | Required |
max_ml_sample max_ml_sample_is_sig3 sig3_genomes_ml max_clust
P-0008709-T01-IM5 PD14457a False 0.255751 SBS3_2
P-0010938-T01-IM5 PD6684a True 0.975855 SBS3_1
Know what changed underneath
- The classifiers. SigMA2's own classifiers were trained on
cosmic_v3p2_inhouseand carry a release status, an evidence grade and a recorded floor per cell (release status); SigMA v1's are Signature 3 classifiers oncosmic_v2_inhousewithout either. The two can call the same tumor differently. - The floor. SigMA v1 applied five SNVs to every tumor type; SigMA2 applies the tumor type's floor. Pass
snv_ctf=5to match R. - MMRD/POLE. SigMA2 ships its own classifiers for five platforms;
model="legacy"keeps SigMA v1's formskandtcga_mc3. The microsatellite counts need the repeat-locus files, the same six files SigMA v1 ships; without them the classifier runs on zero microsatellite indels and prints nothing (Install SigMA2). - HRD. SigMA v1 had no HRD classifier. 3. Detect HRD scores it, and what a call means states how an HRD call differs from a Signature 3 call.
- Removed.
add_sig3, the Shiny interface, and the CSV written beside the input;norm96is still accepted, as the table above says. - What a result records. A SigMA2 run returns its classifier, catalog, floor and cutoffs on the object (Adapt SigMA2 to your own data).