5. Read the output
Read every column of a result table, what a call means at each operating point, and what NA means.
Before you begin, you need a run's output.
Get a result table
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",
)
| 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 |
Read the columns
| Column | Meaning |
|---|---|
total_snvs | SNVs counted for the sample |
<soi>_mva | the classifier's score, in [0, 1]. SBS3_mva, APOBEC_mva, HRD_mva, … |
pass_mva | the call at the released operating point, a 10% false-positive rate |
pass_mva_strict | the call at the strict operating point, 1%. NA when the classifier has no strict cutoff |
confidence | high (strict call), low (loose call only), negative (no call). NA for a loose call when there is no strict cutoff, and for a NaN score |
call_evidence | informative or weak: what a positive call from this classifier is worth, from its card |
exp_<soi>, rat_<soi> | the signature's exposure (mutations) and exposure fraction in the best decomposition |
sigs_all, exps_all | the best decomposition: signatures joined with ., exposures joined with _ |
cutoff_readjusted | the loose cutoff a cohort was called at; present with readjust=True only |
motif_count | the count of the required motif; present when the card requires one (lung APOBEC on msk) |
A sample with fewer SNVs than the floor is absent from out.df of a signature run, for every input format: it was not scored. An HRD run keeps it, with NA features and a call that is not an assessment (Detect HRD). verbose=2 lists these samples with their counts; out.input_matrix.df keeps every sample that was read.
cols = [
"total_snvs",
"SBS3_mva",
"pass_mva",
"pass_mva_strict",
"confidence",
"exp_SBS3",
"rat_SBS3",
"sigs_all",
]
out.df[cols].head(4)
total_snvs SBS3_mva pass_mva pass_mva_strict confidence exp_SBS3 rat_SBS3 sigs_all
P-0008709-T01-IM5 18 0.916549 True True high 6.518799 0.362155 SBS3.SBS8.SBS18.SBS30.SBS85
P-0010938-T01-IM5 41 0.991418 True True high 21.509469 0.524621 SBS2.SBS3.SBS5.SBS13.SBS17b
P-0010560-T01-IM5 5 0.082124 False False negative 0.000000 0.000000 SBS1.SBS5.SBS8.SBS18
P-0008385-T01-IM5 5 0.098018 False False negative 0.000000 0.000000 SBS5.SBS34
Read a call
pass_mvaTrue,pass_mva_strictTrue,confidencehigh: the score clears both cutoffs.pass_mvaTrue,pass_mva_strictFalse,confidencelow: the score clears the loose cutoff only. The loose cutoff lets one in ten of the classifier's held-out simulated negatives through; how many of your own calls are false depends on how many of your tumors lack the signature, and when the signature is rare most calls can be false.- both False,
confidencenegative: no call. pass_mva_strictNA: this classifier has a loose operating point only, andpass_mvais its call. Its strict tier made no call in cross-validation.confidenceis thennegativebelow the loose cutoff and NA above it.- every call column NA: a score-only classifier run with
force=True. The score ranks; nothing is called.
call_evidence comes from the classifier's card, not from your sample: weak means a positive call from this cell changes the odds of the signature little, at every score. Release status explains the grade.
A missing call is NA, never 0 or False: run() neither drops nor zero-fills the call columns.
Optional: the required motif
One card requires a motif: lung APOBEC on msk asks for at least one TCW>T or TCW>G mutation. Below that count pass_mva and pass_mva_strict are False and confidence is negative, the score is unchanged, and motif_count gives the count. A lung panel with no TCW mutation is not called APOBEC-positive, however high its score.
Get the short table with to_lite()
lite = out.to_lite()
list(lite.columns)
['total_snvs', 'SBS3_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']
categ | When |
|---|---|
MMRD, POLE | the MMRD/POLE classifier called it (HRD runs and predict_mmrd_ppd; msi=True does not run that classifier) |
| an MSI or POLE signature | its likelihood is decisive (SigMA v1's rule). to_lite() leaves this category out of a signature run (Migrate from SigMA v1); on predict_mmrd_ppd output it gives msi or Signature_msi to tumors the classifier calls MMRP, four of the fifty example panels, so read that call from MMRD_category |
<soi>_hc | a strict call |
<soi>_lc | a loose call |
| another signature | the most likely other signature, or the most exposed one in the decomposition |
Reference: the working columns
The remaining columns are the features the classifier read. They can be inspected; they are not calls.
| Pattern | Meaning |
|---|---|
<sig>_c | cosine similarity of the spectrum to the signature |
<sig>_l_rat, <sig>_c_diff, <sig>_wl | likelihood ratio, cosine difference, weighted likelihood of the signature |
<soi>_high_<k>_ml, <soi>_low_<k>_ml | likelihood of the tumor type's clusters with and without the signature |
max_ml, sig_max_ml, max_c, sig_max_c, max_wl, sig_max_wl | the best cluster or signature by each measure |
comb_all_c, comb_all_l | the decomposition's cosine and likelihood |
pass_ml | SigMA v1's likelihood-only call; internal, not a call to report |
Reference: HRD and MMRD/POLE columns
run(soi="HRD") adds sig3_score, sig3_ml, sig3_exp, the copy-number features (gis, hrd_loh, tai, lst, or fga, lst_lr), the indel features (mmej, nhej), hrd_group, hrd_variant, hrd_score_pan, hrd_score_type, hrd_score, hrd_call and MMRD_category. Its pass_mva, pass_mva_strict and confidence refer to HRD_mva, and pass_mva is True for HRD-low as well as HRD, so read hrd_call; the Signature 3 step's own calls are SBS3_pass_mva, SBS3_pass_mva_strict and SBS3_confidence. Detect HRD describes them.
predict_mmrd_ppd() returns MMRD_category, MMRD_prob, POLE_prob, MMRP_prob, the indel counts and SigMA v1's subtype columns. Detect MMRD or POLE describes them.
Save the table
out.df.to_csv("scores.tsv", sep="\t")
out.to_lite().to_csv("scores_lite.tsv", sep="\t")
The run's settings are on the object: out.tumor_type, out.seq_type, out.soi, out.catalog_name, out.model, out.snv_floor (the floor applied) and out.snv_ctf (a floor you passed, else None). Record them with the table.