Metadata-Version: 2.4
Name: sigma2
Version: 0.0.16a16
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: License :: Other/Proprietary License
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.9
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Programming Language :: Python :: 3.14
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Requires-Dist: pandas>=1.3.0
Requires-Dist: numpy>=1.20.0
Requires-Dist: scipy>=1.7.0
Requires-Dist: pyranges>=0.0.68
Requires-Dist: tqdm>=4.60.0
Requires-Dist: pyyaml>=5.4.0
Requires-Dist: appdirs>=1.4.0
Requires-Dist: h5py>=3.0
Requires-Dist: joblib>=1.0
Requires-Dist: scikit-learn>=1.4,<1.8
Requires-Dist: tabulate>=0.8
Requires-Dist: pysam>=0.19
Requires-Dist: matplotlib>=3.5
Requires-Dist: seaborn>=0.12
Requires-Dist: kneed>=0.8
Requires-Dist: psutil>=5.8
Requires-Dist: shap>=0.41
Requires-Dist: diptest>=0.10
Requires-Dist: ipython>=8.0 ; extra == 'notebook'
Requires-Dist: pybedtools>=0.9 ; extra == 'panel'
Provides-Extra: notebook
Provides-Extra: panel
License-File: LICENSE
Summary: Mutational signature analysis for low statistics SNV data
Keywords: bioinformatics,genomics,mutational-signatures,cancer
Author-email: Jan P Hummel <jan_hummel@hms.harvard.edu>, Doga C Gülhan <dgulhan@mgh.harvard.edu>
Requires-Python: >=3.9
Description-Content-Type: text/markdown; charset=UTF-8; variant=GFM
Project-URL: Documentation, https://www.sigscape.org/docs#sigma2
Project-URL: Homepage, https://www.sigscape.org/packages/sigma2
Project-URL: Repository, https://github.com/sigscape/SigMA2

SigMA2 - Mutational Signature Analysis in Python
=================================================

SigMA2 detects mutational signatures and homologous-recombination deficiency (HRD) in tumours with few
single nucleotide variants (SNVs): hybrid-capture panels, whole exomes and whole genomes. Trained classifiers
ship with the package: install it, point it at your mutations, and say which question you are asking.

SigMA2 is the successor to [SigMA](https://www.sigscape.org/packages/sigma); the original SigMA (v1)
classifiers are included too, so earlier analyses stay reproducible.

## Choose a classifier

Signatures and HRD go through `sg2.run()`, where the `soi` ("signature of interest") argument picks the
classifier; mismatch-repair deficiency and POLE through `sg2.predict_mmrd_ppd()`:

| You want to know | `soi=` | Classifier | Main output columns |
|---|---|---|---|
| Does the tumour carry Signature 3 (the HRD mutational signature)? | `"SBS3"` | SigMA2 SBS3 | `SBS3_mva`, `pass_mva`, `confidence` |
| Is the tumour homologous-recombination deficient? | `"HRD"` | CSI-HRD: Signature 3 + deletions + copy number | `HRD_mva`, `hrd_call`, `confidence` |
| Does it carry another signature (APOBEC, SBS4 tobacco, SBS7 UV, SBS8, SBS18)? | that signature | SigMA2 signature classifier | `<soi>_mva`, `pass_mva`, `confidence` |
| Is it mismatch-repair (MMR) deficient, or POLE (DNA polymerase epsilon) mutated? | – | `sg2.predict_mmrd_ppd()` | `MMRD_category` |

Ask for HRD rather than SBS3 when the decision is about HR deficiency. On real MSK-IMPACT panels CSI-HRD finds
63 % of biallelic BRCA1/2 breast tumours at a 10 % false-positive rate, SigMA2 SBS3 33 % (see
[HRD](#homologous-recombination-deficiency-hrd)).

## Installation

SigMA2 is an **alpha pre-release**, so `--pre` is required:

```bash
pip install --pre sigma2
```

Python 3.9 or newer. One dependency, `sorted-nearest`, compiles during install, so a C toolchain is needed
(`build-essential` on Debian/Ubuntu, `base-devel` on Arch, the Xcode command line tools on macOS).

## Quickstart

```python
import sigma2 as sg2

# Signature 3 on MSK-IMPACT breast panels
sbs3 = sg2.run(input="panels.maf", ref_genome="hg19", tumor_type="breast", seq_type="msk", soi="SBS3")

# HRD on the same panels, with their copy-number segments
hrd = sg2.run(input="panels.maf", ref_genome="hg19", tumor_type="breast", seq_type="msk", soi="HRD",
              segments="panels.seg")

hrd.df[["HRD_mva", "hrd_call", "confidence"]]
```

`input` is a MAF, a VCF (or a directory of them), or a 96-channel SNV matrix (CSV or DataFrame).
`seq_type` is the sequencing platform:

| `seq_type` | Platform |
|---|---|
| `msk` | MSK-IMPACT |
| `fo` | FoundationOne |
| `op` | OncoPanel |
| `seqcap` | SeqCap EZ exome |
| `tcga_mc3` | TCGA whole exomes (MC3 calls) |
| `wgs` | whole genome (signature reference models) |

## See what is installed

```python
sg2.describe()                          # every classifier per tumour type and platform, with its status
sg2.models(tumor_type="breast")         # "breast/SBS3/msk", "breast/HRD/msk", "breast/MMRD/msk", ...
sg2.model("breast/HRD/msk")             # load one
sg2.verify_model("breast/HRD/msk")      # integrity check (SHA-256 fingerprint)
```

`describe()` shows one table per tumour type, headed by its SNV floor (samples with fewer SNVs are not
scored by any classifier). Rows are classifiers (signatures, HRD and MMRD/POLE alike), columns are platforms:

```
breast  ·  SNV floor 5
        fo    msk    op    seqcap    tcga_mc3    wgs
------  ----  -----  ----  --------  ----------  -----
APOBEC  ●     ●      ●     ●         ●           ●
SBS3    ●     ●      ●     ●         ●           ●
HRD     ●     ●      ●     ●         ●           —
MMRD    —     ●      —     —         ●           —
```

● released · ◐ released, weak evidence: a positive call raises the odds little · ▲ experimental: only its ranking
holds, and `run()` refuses it without `force=True` · `—` none. Weak and experimental cells add the positive
likelihood ratio (LR+) of a call. Tumour types whose tables would be identical share one. The original SigMA
(v1) classifiers, which `model="legacy"` loads, are listed separately.

`models()`, `model()` and `verify_model()` take every kind by the same name, `tumour/classifier/platform`:
`"breast/SBS3/msk"`, `"breast/HRD/msk"`, `"breast/MMRD/msk"`.

Bundled classifiers (the GEL-trained signature classifiers also have `wgs` reference models; the PCAWG-trained ones list `wgs` where it passed the release gate):

| Tumour type | Classifier | Platforms |
|---|---|---|
| breast | SBS3, SBS8, SBS18, APOBEC, HRD | `fo`, `msk`, `op`, `seqcap`, `tcga_mc3` |
| lung | SBS4, APOBEC, HRD | `fo`, `msk`, `op`, `seqcap`, `tcga_mc3` |
| lung | SBS3, SBS18 | `fo`, `op`, `seqcap`, `tcga_mc3` |
| prostate (`prost`) | HRD | `fo`, `msk`, `op`, `seqcap`, `tcga_mc3` |
| prostate (`prost`) | SBS3 | `seqcap`, `tcga_mc3` |
| prostate (`prost`) | SBS18 | `seqcap` |
| skin | SBS7, HRD | `fo`, `msk`, `op`, `seqcap`, `tcga_mc3` |
| colorectal (`crc`) | SBS18, APOBEC, HRD | `fo`, `msk`, `op`, `seqcap`, `tcga_mc3` |
| colorectal (`crc`) | SBS8 | `seqcap`, `tcga_mc3` |
| ovary, bladder, cervix, head and neck (`head`), stomach, uterus | HRD | `fo`, `msk`, `op`, `seqcap`, `tcga_mc3` |
| every other tumour type | HRD, from a model pooled across them | `fo`, `msk`, `op`, `seqcap`, `tcga_mc3`[^hrd] |
| oesophagus (`eso`) | SBS28 | `msk`, `seqcap`, `tcga_mc3`, `wgs` |
| oesophagus (`eso`) | SBS93 | `fo`, `op`, `tcga_mc3` |
| oesophagus (`eso`) | SBS36 | `msk`, `seqcap` |
| kidney | SBS1 | `msk`, `seqcap`, `tcga_mc3`, `wgs` |
| kidney | APOBEC | `msk`, `seqcap`, `tcga_mc3` |
| kidney | SBS40 | `msk`, `wgs` |
| kidney | SBS18 | `tcga_mc3` |
| low-grade glioma (`lgg`) | SBS18 | `wgs` |
| low-grade glioma (`lgg`) | SBS30, SBS37 | `tcga_mc3` |
| liver | SBS4, SBS12 | `fo`, `msk`, `op`, `seqcap`, `tcga_mc3`, `wgs` |
| liver | SBS16 | `fo`, `op`, `seqcap`, `tcga_mc3`, `wgs` |
| liver | SBS1 | `msk`, `seqcap`, `tcga_mc3`, `wgs` |
| liver | SBS22, SBS36 | `msk`, `seqcap`, `tcga_mc3` |
| liver | SBS8 | `fo`, `msk`, `op` |
| liver | SBS85 | `seqcap`, `tcga_mc3`, `wgs` |
| liver | SBS18 | `tcga_mc3`, `wgs` |
| liver | SBS30 | `tcga_mc3` |
| lymphoma (`lymph`) | SBS36 | `msk`, `seqcap`, `tcga_mc3`, `wgs` |
| lymphoma (`lymph`) | SBS17, SBS84 | `msk`, `seqcap`, `tcga_mc3` |
| lymphoma (`lymph`) | SBS9, SBS30 | `seqcap`, `tcga_mc3`, `wgs` |
| lymphoma (`lymph`) | APOBEC | `seqcap`, `tcga_mc3` |
| lymphoma (`lymph`) | SBS34 | `msk`, `seqcap` |
| lymphoma (`lymph`) | SBS85 | `msk`, `tcga_mc3` |
| lymphoma (`lymph`) | SBS89 | `wgs` |
| medulloblastoma (`medullo`) | SBS8 | `seqcap`, `tcga_mc3`, `wgs` |
| medulloblastoma (`medullo`) | SBS3 | `seqcap`, `wgs` |
| pancreatic adenocarcinoma (`panc_ad`) | SBS3 | `fo`, `op`, `seqcap`, `wgs` |
| pancreatic adenocarcinoma (`panc_ad`) | SBS9, SBS17 | `msk`, `seqcap`, `tcga_mc3` |
| pancreatic adenocarcinoma (`panc_ad`) | SBS30 | `tcga_mc3`, `wgs` |
| pancreatic adenocarcinoma (`panc_ad`) | SBS36 | `seqcap`, `tcga_mc3` |
| pancreatic adenocarcinoma (`panc_ad`) | SBS18 | `tcga_mc3` |
| pancreatic neuroendocrine (`panc_en`) | SBS36 | `msk`, `seqcap`, `tcga_mc3`, `wgs` |
| pancreatic neuroendocrine (`panc_en`) | APOBEC | `seqcap`, `tcga_mc3` |
| pancreatic neuroendocrine (`panc_en`) | SBS18 | `tcga_mc3`, `wgs` |
| pancreatic neuroendocrine (`panc_en`) | SBS30 | `seqcap`, `wgs` |
| stomach | SBS28 | `msk`, `seqcap`, `tcga_mc3` |
| stomach | SBS41 | `msk`, `op` |
| every tumour type | MMRD/POLE | `fo`, `msk`, `op`, `seqcap`, `tcga_mc3` |

[^hrd]: Except tumour types a platform's training cohort did not contain, which `run()` refuses; for example
uveal melanoma on `msk`.

## Mutational signatures

`model="auto"` (the default) picks the classifier for your tumour type, signature and platform.

**Output.** `<soi>_mva` is the score. `pass_mva` is a call at the cutoff for a 10 % false-positive rate and
`pass_mva_strict` at 1 %. `confidence` is high (past 1 %), low (past 10 % only) or negative.

**Experimental classifiers.** Some combinations ship as experimental: lung SBS3 and SBS18 on `msk`, most
prostate panel cells, and colorectal SBS8 on `fo`, `msk` and `op`. Their ranking is usable but their calls are
not, so `run()` refuses them unless you pass `force=True`.

**What a classifier calls** is set by its training labels, recorded on its model card (`label_curation`).
- The breast SBS3 classifiers on `msk` and `tcga_mc3` call strong SBS3 only. They are **not intended to call
  weak SBS3**.
- The colorectal SBS18 classifiers call the summed SBS18 + SBS36 process.

`readjust=True` re-derives a model's cutoffs for the SNV depth of the cohort you score, as SigMA did.

**SigMA v1.** `model="legacy"` runs the original SigMA classifiers (Signature 3 only), for reproducing
earlier analyses:

```python
sg2.run(input="...", model="legacy", soi="SBS3", tumor_type="breast", seq_type="msk")
```

## Homologous-recombination deficiency (HRD)

`soi="HRD"` runs CSI-HRD (Gulhan et al., bioRxiv 10.1101/2025.10.13.681986). It combines three kinds of
evidence:

- SigMA's Signature 3 score, likelihood and exposure;
- deletions of at least 5 bp with microhomology (`mmej`) and without (`nhej`);
- copy-number instability.

CSI stands for copy number, single-base substitutions and indels. Each tumour gets a pan-cancer score and a
tumour-type score. The call combines them:

| `hrd_call` | Rule | `pass_mva` | `pass_mva_strict` | `confidence` |
|---|---|---|---|---|
| HRD | both scores > 0.5 | True | True | high |
| HRD-low | one score > 0.5 | True | False | low |
| HRP (HR-proficient) | neither | False | False | negative |
| MMRD / POLE | hypermutated (MMR-deficient or POLE-mutated), set aside first | – | – | – |

`HRD_mva` is the mean of the two scores; `hrd_score_pan` and `hrd_score_type` are the parts. The Signature 3
run behind it stays in the output as `SBS3_*`.

**What you pass.**

- **Mutations.** A VCF or MAF, so SigMA2 can count the deletions. With a matrix input, add `mmej` and `nhej`
  columns, or pass `indels=`.
- **Copy number,** optional but strongly recommended. Pass `segments=`; SigMA2 picks the model variant from
  what the file holds:
  - **Allele-specific** segments (ASCAT, FACETS): the three parts of the genomic instability score (GIS) --
    loss-of-heterozygosity regions (HRD-LOH), telomeric allelic imbalance (TAI) and large-scale state
    transitions (LST). On `msk` they are read on the IMPACT gene grid, as MSK's gene-level calls are.
  - **Log-ratio** segments (a panel `.seg`): `fga` (fraction of the genome with |log ratio| ≥ 0.2) and
    `lst_lr` (large-scale transitions).
  - **None:** a model without copy number, which is weaker.
- `hrd_loh` / `tai` / `lst` or `fga` / `lst_lr` columns in a matrix input are used as given. A GIS sum alone
  is not enough for the allele-specific variant; SigMA2 warns and falls back.

**What `run()` refuses, before computing anything:**
- `seq_type="auto"`, or a platform without an HRD model;
- a tumour type that platform's model never saw;
- `hrd=False`, `mva=False` and `readjust=True`, which do not apply to HRD.

**Platforms and tumour types:** see [See what is installed](#see-what-is-installed).

**How it performs.** Real MSK-IMPACT 50K panels, biallelic BRCA1/2 truth, scored out of fold; AUC:

| Tumours (biallelic) | SigMA SBS3 | Best single input | CSI-HRD |
|---|---|---|---|
| All 12 types (288) | 0.58 | 0.87 | 0.90 |
| Breast (72) | 0.63 | 0.87 | 0.88 |
| Ovary (85) | 0.60 | 0.82 | 0.83 |
| Prostate (76) | 0.57 | 0.78 | 0.81 |

On TCGA exomes, out of fold, CSI-HRD reaches 0.95, against 0.82 for SigMA's Sig3 score.

There is no public real-panel truth for `seqcap`, `fo` and `op`. Their models are trained and tested on TCGA
exomes reduced to each footprint and reach 0.95–0.96. `describe()` shows each tumour type's cross-validated
AUC, so weak ones are visible. For example, colorectal on `msk` has an AUC of 0.48: only 5 biallelic tumours.

`hrd=True` adds the HRD columns to a run for another `soi`.

## Mismatch-repair deficiency and POLE

```python
mmrd = sg2.predict_mmrd_ppd(input="cohort.maf", ref_genome="hg19", tumor_type="uterus",
                            seq_type="msk", soi="SBS3", catalog="cosmic_v2_inhouse", mva=False)
mmrd.df[["MMRD_category", "MMRD_prob", "POLE_prob", "MMRP_prob"]]
```

`MMRD_category` is MMRD, POLE (a POLE exonuclease-domain mutator) or MMRP, one call per tumour at a 1%
false-positive rate; tumours below the SNV floor are "Unclassified". SigMA2 has a classifier for `msk`, `fo`,
`op`, `seqcap` and `tcga_mc3`; `model="legacy"` runs SigMA v1's (`msk`, `tcga_mc3`). Both read the same
features; SigMA2's are retrained on larger truth sets:

| Platform | Trained on | MMRD / POLE / MMRP tumours |
|---|---|---|
| `msk` | real MSK-IMPACT 50K panels, MSIsensor calls | 1,177 / 138 / 32,955 |
| `tcga_mc3` | TCGA exomes, MANTIS and MSIsensor agreeing | 238 / 42 / 7,337 |
| `seqcap`, `fo`, `op` | the same exomes reduced to each footprint | 238 / 42 / about 7,000 |

POLE truth: a pathogenic exonuclease-domain hotspot (P286R, S297F, F367S, V411L, L424I, P436R, M444K, A456P,
S459F) in a microsatellite-stable tumour. Out of fold, against SigMA v1 on the same tumours, sensitivity at
1% FPR on real MSK-IMPACT panels: MMRD 0.997 vs 0.994, POLE 0.993 vs 0.855. Tumours called MSI-high whose input
holds no microsatellite indel are left out of these comparisons for both classifiers.

It reads the indel counts SigMA2 takes from a VCF or MAF (`sigma2.indel_features`):

| Count | Meaning |
|---|---|
| `nins`, `ndel` | all insertions and deletions |
| `nmsi_ins`, `nmsi_del` | insertions and deletions in a microsatellite locus |
| `mmej`, `nhej` | deletions ≥ 5 bp with and without ≥ 2 bp microhomology |

VCF-style (anchor base) and MAF-style (`-`) records give the same counts. The repeat loci in the data
directory are hg19; on hg38 pass `beds=` in that build.

## Train your own

- **Signature classifier:** the `Add` workflow annotates, clusters and simulates your cohort, then tunes a
  classifier you pass back as `model=`. See the example notebooks.

  ```python
  a = sg2.Add(input=mutations, exposures=exposures, tumor_type="breast", seq_type=["msk"], soi=["SBS3"])
  a.annotate(); a.cluster(); a.simulate(); models = a.tune(output_dir="...")
  ```

- **HRD classifier:** for another assay or cohort, build the feature table and give biallelic BRCA1/2 status
  per tumour. `Add` does not train HRD; HRD is a classifier, not a signature.

  ```python
  F = sg2.hrd_features(input="cohort.maf", tumor_type="breast", seq_type="my_assay",
                       segments="segments.tsv", ref_genome="hg19")
  F["brca"] = brca_biallelic                            # bool per sample
  sg2.tune_hrd(F, labels="brca", seq_type="my_assay")   # then run(..., seq_type="my_assay", soi="HRD")
  ```

From the shell: `sigma2 hrd score | features | tune | gis | models`.

## Reference data

The package ships everything a run needs in the way of models: the data store,
the classifiers, and the panel definitions. Bulk reference data (~835 MB) is
downloaded on demand instead:

```bash
sigma2 data list          # what exists and what is installed
sigma2 data fetch --all   # download everything
sigma2 data verify        # re-check checksums
sigma2 info               # versions, data directory, disk usage
```

Files land in a per-user data directory; set `$SIGMA2_DATA_DIR` to override it,
which is the usual choice on a cluster where one pre-seeded copy is shared.

Reference genomes are separate and larger (~3 GB each): `sigma2 data genome hg19`.

Two things are deliberately not distributed, because they are not ours to
redistribute: reference genomes, fetched from UCSC by the command above, and
gnomAD germline VCFs, which you supply yourself. Obtain the latter from
[gnomAD](https://gnomad.broadinstitute.org/downloads) and point SigMA2 at your
own copy.

> **Access note.** Until the accompanying paper is published, the bulk data
> lives in a private repository. `sigma2 data fetch` therefore needs GitHub
> credentials with access to it — either run `gh auth login`, or set
> `SIGMA2_GITHUB_TOKEN`. Without access the fetch fails with a 404; the code
> itself is fully functional, and any data you already hold can be used by
> pointing `$SIGMA2_DATA_DIR` at it. Please contact the authors (below) to
> request access.

## Example notebooks

Work through the [example notebooks](https://github.com/sigscape/SigMA2/tree/master/examples)
in numerical order to get started with SigMA2 and its built-in functions.

They are not included in the installed package; clone the repository to run them
locally:

```bash
git clone https://github.com/sigscape/SigMA2.git
jupyter notebook SigMA2/examples
```

## Documentation & Further Information
Find more [here](https://www.sigscape.org/packages/sigma2)

## Contact

Requests for use of the Software for or on behalf of for-profit entities or for any commercial purposes, please contact:

Office of Technology Development
Harvard University
Smith Campus Center, Suite 727E
1350 Massachusetts Avenue
Cambridge, MA 02138 USA
Telephone: (617) 495-3067
E-mail: otd@harvard.edu

## Citation

**Citation pending.** The manuscript describing SigMA2 is in preparation; this
section will be updated with the reference once the preprint is available.

If you use SigMA2 in your work in the meantime, please contact the authors for
the citation to use.

---

Copyright (c) 2019-2026 - President and Fellows of Harvard College. All rights reserved.

