1. Prepare your input files

run() and predict_mmrd_ppd() read mutations from a MAF, from VCFs, or from a matrix of SBS-96 counts. HRD adds copy-number segments. Every input is read by column name, and a missing column is named in the error. Two things pass without a message: mutations outside chromosomes 1 to 22, X and Y are left out, and a reference allele is not compared with the genome (check the reference alleles).

Prepare a MAF

A tab-separated MAF, plain or gzipped (.maf, .maf.gz), one row per mutation, many tumors in one file.

A TRUE allele (Excel's or R's reading of T) is mapped back to T. A single-tumor MAF works.

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",
)

Prepare VCFs

One file (.vcf, or .vcf.gz and .vcf.bgz compressed with bgzip), a list of paths, or a folder (one sample per file). A file compressed with gzip fails with "seek not implemented in files compressed by method 1".

• Every SNV record counts, whatever its FILTER. Filter to PASS before scoring. Indels are read only when PASS or ..
• Multi-allelic records are split into one row per ALT allele.
• A VCF with variants in several sample columns gives one sample per column that carries a non-reference genotype (FORMAT GT), named by the column. A tumor/normal pair, with the normal 0/0, stays one sample. Two sample columns without GT (Strelka-style) are read as one sample, without a message at the default verbosity.
• A single-sample VCF is named by its file name without the extension (T1.vcf.gz is T1). Files with the same name in different folders keep their paths. An exposure table or a segments table matches a folder of VCFs by those names.
• Sample IDs are strings.

sg2.run(
    input="calls/",
    ref_genome="hg38",
    tumor_type="lung",
    seq_type="tcga_mc3",
    soi="SBS4",
)
sg2.run(
    input=["a.vcf.gz", "b.vcf.gz"],
    ref_genome="hg19",
    tumor_type="breast",
    seq_type="msk",
    soi="SBS3",
)

Prepare an SBS-96 matrix

A CSV path or a pandas.DataFrame, samples as rows and the 96 channels as columns. No reference genome is needed.

• Channels are found by name, in any spelling and case: COSMIC A[C>A]A, the package's a[c>a]a, or SigMA v1's caaa. A missing or repeated channel is refused with the list.
• A transposed matrix (channels as rows) is detected and flipped.
• Other columns (total_snvs, indel counts, …) are kept wherever they sit.
• Sample IDs come from the index, from the first text column of a CSV, or from SigMA's tumor column.

spectra = sg2.Samples("test_mutations_50sample.maf", ref_genome="hg19").df.iloc[:, :96]
spectra.to_csv("spectra.csv")
out = sg2.run(input="spectra.csv", tumor_type="breast", seq_type="msk", soi="SBS3")

MMRD/POLE on a matrix needs indel counts in it: nins, ndel, nmsi_ins, nmsi_del. Without them the error names the missing columns.

Prepare copy-number segments for HRD

A segments table beside the mutations, passed as segments= (path or DataFrame). A sample column is required; add one to single-sample files such as PURPLE output. Column names are matched without regard to case:

Allele-specific segments give the genomic instability scores (HRD-LOH, TAI, LST). Recognized column pairs:

Log-ratio segments (cBioPortal and IMPACT .seg files, seg.mean and its aliases) give the tcn variant of the HRD model. Allele-specific segments with purity and ploidy columns are converted to a log ratio when the model needs one.

sample        chr  startpos  endpos    nMajor  nMinor
TCGA-A2-A04P  1    61735     6656960   4       0
TCGA-A2-A04P  1    6657355   8027397   2       0

A table with no recognized copy-number column is an error listing the accepted names; copy number that matches no tumor is a warning naming the samples; a tumor without segments is scored by the nogis variant. Detect HRD describes the variants.

Check the genome build

ref_genome takes hg19, hg38, GRCh37, GRCh38, or a FASTA path. On Windows the FASTA must be plain text; the genome sigma2 data genome installs is.

Check the reference alleles

The reader takes the trinucleotide context from the genome at each position and does not compare the file's reference allele with the base there, so a file with 0-based positions, or in another build without an NCBI_Build column, is scored without a message. Shifting the example file's positions by one base changed 16 of its 38 Signature 3 calls. Compare the alleles with the genome before scoring:

import pandas as pd
import pysam

maf = pd.read_csv(
    "test_mutations_50sample.maf",
    sep="\t",
    comment="#",
    dtype=str,
    keep_default_na=False,
)
snv = maf[maf["Variant_Type"] == "SNP"]
genome = pysam.FastaFile("/path/to/data_dir/reference/hg19_genome.fa")


def base(chrom, pos):
    name = chrom if chrom.startswith("chr") else "chr" + chrom
    return genome.fetch(name, pos - 1, pos).upper() if name in genome.references else None


found = [base(c, int(p)) for c, p in zip(snv["Chromosome"], snv["Start_Position"])]
agree = sum(g == r.upper() for g, r in zip(found, snv["Reference_Allele"]) if g is not None)
print(f"{agree} of {sum(g is not None for g in found)} reference alleles agree with the genome")

On the example file it prints 729 of 729 reference alleles agree with the genome; the same file with every position lowered by one base prints 129 of 729.

What is not accepted

• A custom panel name. seq_type accepts the vocabulary only; score a custom panel as the closest bundled platform.
• A VCF without positions. Matrices skip the genome; every other input needs coordinates and a build.
• Germline calls. SigMA2 expects somatic mutations.