Reproduce SigMA v1
Reproduce a SigMA (R) result in Python, and compare it with the SigMA2 classifier on the same tumors.
Before you begin, you need the package installed, the hg19 genome.
Get the example file
curl -LO https://raw.githubusercontent.com/parklab/SigMA/master/inst/extdata/examples/test_mutations_50sample.maf
Run SigMA v1's classifier
import sigma2 as sg2
out = sg2.run(input="test_mutations_50sample.maf", ref_genome="hg19", tumor_type="breast",
seq_type="msk", soi="SBS3", model="legacy", 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 |
| modelstr or GBMModel | Which classifier runs: "auto" for the released cell, "legacy" for SigMA v1's, or a model loaded with sg2.model(). | 'auto' |
| 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 |
[SigMA2|v0.0.16-alpha.20|run] Using SBS3 classifier from SigMA v1 (legacy, native Python)
Check the legacy grid
All SigMA v1 classifiers are Signature 3 classifiers. By platform, with the features SigMA2's port passes to each as missing:
| platform | reads every feature | a feature passed as missing | does not run |
|---|---|---|---|
| msk | bladder, breast, eso, ovary, prost, stomach | crc (APOBEC_c1, rat_sig3), osteo (clock_c1, clock_c2, clock_c3), panc_ad (APOBEC_c1), panc_en (APOBEC_c1), uterus (APOBEC_c1, APOBEC_c2) | none |
| msk_cf | bladder, breast, eso, ovary, panc_ad, panc_en, prost, stomach | osteo (clock_c1, clock_c2, clock_c3), uterus (APOBEC_c1) | none |
| seqcap | bladder, breast, eso, ewing, ovary, panc_ad, stomach | medullo (rat_sig3), osteo (APOBEC_c1), panc_en (APOBEC_c1), prost (rat_sig3), uterus (APOBEC_c1) | none |
| seqcap_probe | bladder, eso, ewing, ovary, panc_ad, prost, stomach | breast (APOBEC_c1, APOBEC_c2), medullo (rat_sig3), osteo (APOBEC_c1), panc_en (APOBEC_c1), uterus (APOBEC_c1) | none |
| tcga_mc3 | ovary | adrenal (rat_sig3), biliary (APOBEC_c1, rat_sig3), bladder (APOBEC_c2, rat_sig3), breast (APOBEC_c1), cervix (APOBEC_c3, rat_sig3), crc (APOBEC_c1, rat_sig3), eso (APOBEC_c1, APOBEC_c2, rat_sig3), glioma (rat_sig3), head (APOBEC_c1, rat_sig3), kidney (rat_sig3), lgg (rat_sig3), liver (rat_sig3), lung (rat_sig3), lymph (rat_sig3), meso (rat_sig3), panc_ad (rat_sig3), prost (rat_sig3), skin (APOBEC_c1, oxi_c1, rat_sig3), stomach (rat_sig3), uterus (APOBEC_c1, APOBEC_c2, rat_sig3) | none |
| wgs | none | biliary (rat_sig3), bladder (rat_sig3), bone_other (rat_sig3), breast (rat_sig3), cervix (rat_sig3), crc (rat_sig3), eso (APOBEC_c2, rat_sig3), ewing (rat_sig3), lms (rat_sig3), lung (APOBEC_c3, rat_sig3), medullo (rat_sig3), osteo (rat_sig3), ovary (rat_sig3), panc_ad (rat_sig3), panc_en (rat_sig3), prost (rat_sig3), skin (rat_sig3), stomach (rat_sig3), uterus (APOBEC_c1, rat_sig3) | kidney |
| wgs_pancan | breast | none | generic |
A classifier in the middle column scores without that feature, through the missing branch of each of its trees, while R computes the feature, so its scores and calls can differ from R's; Migrate from SigMA v1 gives the size of the difference on TCGA breast exomes. With verbose=1 a run prints the features it passed as missing and adds that SigMA v1 did the same; R's output for those exomes carries the feature, and supplying it reproduces R's scores. The classifiers in the last column stop with an IndexError before they score.
Plus the msk_cf variants on msk (below) and SigMA v1's MMRD/POLE multinomials for msk and tcga_mc3. sg2.describe() prints the grid under "SigMA v1 Classifiers".
sg2.legacy_store.available()
Set the four settings that reproduce R
1. Catalog. SigMA v1 was trained on cosmic_v2_inhouse, and model="legacy" uses it by default. Another catalog produces a warning.
2. The cluster-size prior (weight_cf). R weights the cluster likelihoods by cluster size on some platforms. SigMA2 does the same: the prior is always on for seqcap and seqcap_probe, and your choice on msk, where R offered the msk_cf classifiers:
cf = sg2.run(input="test_mutations_50sample.maf", ref_genome="hg19", tumor_type="breast",
seq_type="msk", soi="SBS3", model="legacy", wcf=True)
3. The SNV floor. R applies a flat floor of five SNVs to every tumor type. A legacy run in SigMA2 applies the tumor type's floor, which is 4 for prostate and 3 for osteosarcoma and pancreatic neuroendocrine tumors. To match R on those three, pass the flat floor:
sg2.run(
input="prostate.maf",
ref_genome="hg19",
tumor_type="prost",
seq_type="msk",
soi="SBS3",
model="legacy",
snv_ctf=5,
)
4. MMRD/POLE. R's multinomial classifiers run through predict_mmrd_ppd(model="legacy"), for msk and tcga_mc3. The microsatellite counts need the repeat-locus files (Install SigMA2); without them the run proceeds on zero microsatellite indels and prints nothing. A matrix input with an msisensor column (MSIsensor's score per sample) selects the R model trained with that score, as combine_SBS_indel(msisensor_file=) did:
mm = sg2.predict_mmrd_ppd(
input="test_mutations_50sample.maf",
ref_genome="hg19",
tumor_type="breast",
seq_type="msk",
model="legacy",
)
| 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 |
| modelstr or GBMModel | "auto" for SigMA2's model of the platform, "legacy" for SigMA v1's (msk and tcga_mc3 only). | 'auto' |
Read the output
The same columns as a SigMA2 run, in SigMA v1's vocabulary: SBS3_mva, pass_mva, pass_mva_strict, confidence, exp_SBS3, rat_SBS3, sigs_all, exps_all. The cluster likelihoods carry v1's cluster names (SBS3_1_ml, APOBEC_1_ml, clock_1_ml, …), and to_lite() sums them per signature:
list(out.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']
call_evidence is absent: a legacy card carries no evidence grade.
Compare with the SigMA2 classifier
The same tumors through both classifiers:
kw = dict(
input="test_mutations_50sample.maf",
ref_genome="hg19",
tumor_type="breast",
seq_type="msk",
soi="SBS3",
verbose=0,
)
new = sg2.run(**kw).df
old = sg2.run(**kw, model="legacy").df
both = new[[
"total_snvs",
"SBS3_mva",
"pass_mva",
]].join(old[["SBS3_mva", "pass_mva"]].add_suffix("_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
For a classifier that reads every feature, the legacy path reproduces R's short table and Signature 3 columns on the same input; for the others, the table under Check the legacy grid names what differs. The example notebook examples/7_sigma1_concordance.ipynb in the SigMA2 repository compares the two on TCGA breast exomes.
Simulate panels like your cohort
SigMA v1's quick_simulation brings the package's whole genomes of the tumor type, of known Signature 3 status, to your cohort's depth, so a cutoff can be read at a false-positive rate on panels that look like yours. Only Signature 3 is simulated:
sim = sg2.quick_simulation(
input="test_mutations_50sample.maf",
ref_genome="hg19",
tumor_type="breast",
seq_type="msk",
seed=1,
)
print(sim.shape)
print(sim["is_sig3"].value_counts())
sim.iloc[:, :96].sum(axis=1).describe().round(1)
| 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 |
(671, 100)
is_sig3
False 471
True 200
Name: count, dtype: int64
count 671.0
mean 10.9
std 13.1
min 0.0
25% 4.0
50% 8.0
75% 14.0
max 168.0
dtype: float64
Score the simulated panels with a classifier and keep the truth beside the score; a cutoff is then read off the negatives at the false-positive rate you accept:
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"])
scored[["total_snvs", "SBS3_mva", "pass_mva", "is_sig3"]].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 |
| 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 |
total_snvs SBS3_mva pass_mva is_sig3
SP7692 50 0.044079 False False
SP5052 32 0.984625 True True
SP7785 20 0.671830 True True
Choose between legacy and SigMA2
- Reproducing a SigMA (R) result:
model="legacy", with the four settings above, on a classifier that reads every feature (Check the legacy grid); on the others the legacy run is SigMA v1's classifier without the features it lacks. - A tumor type or signature SigMA2 has no released cell for, but SigMA v1 has:
model="legacy"is a validated classifier where there is no other. Therun()refusal for a missing cell says when this applies. - Everything else:
model="auto". The SigMA2 cells cover more signatures and carry a status, an evidence grade and a recorded floor; the legacy cells carry none of these.