Worked example: HRD on exomes

Call homologous-recombination deficiency on TCGA breast exomes with allele-specific copy-number segments, and compare the HRD call with the Signature 3 call on the same tumors. The exome is tcga_mc3; the segments are ASCAT's.

You need the hg19 genome and the repeat-locus files for the MMRD/POLE step an HRD run starts with (sigma2 data fetch --all for a GitHub account the maintainers have given access, or from the public SigMA (R) repository). The mutations are SigMA (R)'s public TCGA breast MAF; the segments, for six of its tumors, come from this site.

Note These six tumors are in the training cohort of the tcga_mc3 HRD model, and three of them, TCGA-A2-A04P, TCGA-AO-A0J2 and TCGA-C8-A134, were labeled positives through BRCA1 silencing. The example shows how to run and read the model, not how it performs on tumors it has not seen.

1. Get the two inputs

curl -LO https://raw.githubusercontent.com/parklab/SigMA/master/inst/extdata/examples/tcga_mc3_brca.maf
curl -LO https://sigscape.org/sigma2/examples/tcga_brca_ascat_segments.tsv

The MAF holds about a thousand TCGA breast exomes (NCBI_Build GRCh37, Variant_Type SNP, INS and DEL, so indels are read from it). The segments table:

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

This example restricts the MAF to the tumors that have segments, so the tables stay short and the run takes seconds. segments= itself restricts nothing: a tumor without segments is scored without copy number. On a whole cohort pass mxw="max"; an HRD run costs about a second per tumor on one worker.

import pandas as pd, sigma2 as sg2

seg = pd.read_csv("tcga_brca_ascat_segments.tsv", sep="\t")
maf = pd.read_csv("tcga_mc3_brca.maf", sep="\t", comment="#", low_memory=False)
maf[maf["Tumor_Sample_Barcode"].isin(seg["sample"])].to_csv(
    "brca_subset.maf",
    sep="\t",
    index=False,
)

2. Compute genomic instability on its own

The copy-number features before the classifier:

gis = sg2.genomic_instability("tcga_brca_ascat_segments.tsv", genome="hg19")
gis

              hrd_loh  tai  lst  gis
TCGA-A2-A04P       20   24   24   68
TCGA-AO-A0J2       22   29   19   70
TCGA-C8-A134       21   22   21   64
TCGA-AO-A0J8        1    2    0    3
TCGA-AN-A0FN        0    0    0    0
TCGA-LD-A9QF        0    0    0    0

sigma2 hrd gis --segments tcga_brca_ascat_segments.tsv -o gis.tsv is the same from the shell.

3. Score HRD

hrd = sg2.run(input="brca_subset.maf", ref_genome="hg19", tumor_type="breast", seq_type="tcga_mc3",
              soi="HRD", segments="tcga_brca_ascat_segments.tsv", verbose=1)
cols = [
    "total_snvs",
    "sig3_score",
    "gis",
    "mmej",
    "nhej",
    "hrd_variant",
    "hrd_score",
    "hrd_call",
    "MMRD_category",
]
hrd.df[cols]

              total_snvs  sig3_score   gis  mmej  nhej hrd_variant  hrd_score hrd_call MMRD_category
TCGA-A2-A04P         169    0.976025  68.0   4.0   3.0         gis   0.998613      HRD          MMRP
TCGA-AN-A0FN          53    0.024308   0.0   0.0   0.0         gis   0.000092      HRP          MMRP
TCGA-AO-A0J2          82    0.889912  70.0   2.0   1.0         gis   0.999091      HRD          MMRP
TCGA-AO-A0J8          33    0.015610   3.0   0.0   1.0         gis   0.000208      HRP          MMRP
TCGA-C8-A134         117    0.975525  64.0   2.0   1.0         gis   0.997923      HRD          MMRP
TCGA-LD-A9QF          40    0.026845   0.0   0.0   0.0         gis   0.000083      HRP          MMRP

• hrd_variant is gis: allele-specific segments were found, so the full model ran. Without segments it would say nogis; with a .seg file, tcn.
• sig3_score is the Signature 3 step's score (SigMA v1's SBS3 classifier on cosmic_v2_inhouse, as recorded on the HRD model); mmej and nhej count the microhomology-mediated and non-homologous end-joining deletions.
• MMRD_category is MMRP for all six: the pre-step found no mismatch-repair deficient or POLE tumor, so every tumor was scored for HRD.
• hrd_call is HRD for the three tumors with scars and a high Signature 3 score, HRP for the three without. Breast has its own HRD model on tcga_mc3 (hrd_group is breast, not pooled).

The run also carries HRD_mva with pass_mva, pass_mva_strict and confidence, and the Signature 3 step's own calls as SBS3_pass_mva, SBS3_pass_mva_strict, SBS3_confidence.

4. Score Signature 3 alone

Note This cell's card records a known limit: its exome panels were simulated at hg19 coordinates from genomes aligned to hg38, and it is to be retrained (known limits).

sbs3 = sg2.run(
    input="brca_subset.maf",
    ref_genome="hg19",
    tumor_type="breast",
    seq_type="tcga_mc3",
    soi="SBS3",
)
sbs3.df[["total_snvs", "SBS3_mva", "pass_mva", "pass_mva_strict", "confidence"]]

              total_snvs  SBS3_mva  pass_mva  pass_mva_strict confidence
TCGA-A2-A04P         169  0.841614      True             True       high
TCGA-AN-A0FN          53  0.000825     False            False   negative
TCGA-AO-A0J2          82  0.177295      True            False        low
TCGA-AO-A0J8          33  0.000133     False            False   negative
TCGA-C8-A134         117  0.554682      True             True       high
TCGA-LD-A9QF          40  0.000062     False            False   negative

The SigMA2 SBS3 classifier for tcga_mc3 agrees on five tumors and gives the sixth a loose call. The HRD run calls that tumor HRD, and its own Signature 3 input, SigMA v1's classifier, called it strictly (SBS3_confidence is high in the HRD run of step 3); without segments the HRD call is the same, so it does not rest on copy number. The two Signature 3 classifiers disagree on this tumor, and the SigMA2 cell's card records a known limit in its training.

Report hrd_call with a Signature 3 call beside it, and name the classifier it came from. HRD-low, had it appeared, would mean the pan-cancer and the tumor-type model disagree, and the result is indeterminate.

Look at the six tumors

The six exomes as stacked exposures sorted by the Signature 3 share, and the spectra of two HRD calls: TCGA-A2-A04P, with a high Signature 3 score, and TCGA-AO-A0J2, which SigMA2's Signature 3 classifier gave only a loose call. The spectra show the point mutations alone; the HRD model read their deletions and copy number as well:

from sigma2 import plot as sp

sp.exps(sbs3.df, soi=sg2.resolve_signature("SBS3"), path="exome_exposures.svg")
for sid in ["TCGA-A2-A04P", "TCGA-AO-A0J2"]:
    r, h = sbs3.df.loc[sid], hrd.df.loc[sid]
    sp.sbs(
        sbs3.df.loc[[sid]].iloc[:, :96],
        title=f"{sid}: Signature 3 score {r['SBS3_mva']:.2f}, {r['confidence']}; "
        f"HRD call {h['hrd_call']}",
        path=f"{sid}.svg",
    )

5. Run the same thing in two steps

feats = sg2.hrd_features("brca_subset.maf", tumor_type="breast", seq_type="tcga_mc3",
                         segments="tcga_brca_ascat_segments.tsv", ref_genome="hg19")
scores = sg2.score_hrd(feats, seq_type="tcga_mc3", tumor_type="breast")
list(feats.columns), list(scores.columns)

(['sig3_score', 'sig3_ml', 'sig3_exp', 'gis', 'hrd_loh', 'tai', 'lst', 'fga', 'lst_lr', 'mmej', 'nhej', 'tumor_type'],
 ['hrd_group', 'hrd_variant', 'hrd_score_pan', 'hrd_score_type', 'hrd_score', 'hrd_call'])

Keep the feature table when you will score more than once, inspect a feature, or train your own model with sg2.tune_hrd from tumors with biallelic BRCA1/2 truth.

Optional: from the shell

sigma2 hrd score --input brca_subset.maf --tumor-type breast --seq-type tcga_mc3 \
                 --segments tcga_brca_ascat_segments.tsv -o hrd.tsv

Take it to your own data

• hg38 segments work on tcga_mc3, seqcap, fo and op; the msk model's gene grid is hg19 and refuses them.
• Segments from a panel on fo, op or seqcap are coarser than the SNP-array segments those models learned from; report the call without copy number beside the call with it.
• A single-sample file (PURPLE) needs a sample column added.
• Pooled tumor types rest on other tumor types' training data; hrd_group says pooled when that is the case.
• Tumor types the model cannot score are refused with the reason (Detect HRD).