Package 'mSigSpectra'

Description

Extracts seq_context_width bases upstream and downstream of each variant position from the reference genome and attaches them as a new column named ⁠seq.<N>bases⁠ where N = 2 * seq_context_width + 1.

Usage

add_seq_context(df, ref_genome, seq_context_width = 10, name_of_vcf = NULL)

Arguments

Value

df with an added character column ⁠seq.<N>bases⁠.

Description

For each variant, finds overlapping transcripts in trans_ranges via GenomicRanges::findOverlaps() with type = "within", and appends columns trans.start.pos, trans.end.pos, trans.strand, trans.Ensembl.gene.ID, trans.gene.symbol, plus bothstrand (TRUE if the variant falls on transcripts from both strands) and count (number of overlapping transcripts).

Usage

add_transcript_strand(df, ref_genome, trans_ranges = NULL, name_of_vcf = NULL)

Arguments

Details

If trans_ranges is NULL and no shipped table is available for the given ref_genome, returns df as a data.table unchanged (no strand columns added). If df has zero rows, it is returned unchanged.

Value

A data.table with the annotation columns added (left-join semantics: variants outside any transcript have NA values).

Description

Nested list of k-mer counts keyed by BSgenome.Hsapiens.1000genomes.hs37d5 / BSgenome.Hsapiens.UCSC.hg38 / BSgenome.Mmusculus.UCSC.mm10, then by exome / transcript / genome, then by catalog size as a string ("78", "96", "136", "144", "192", "1536"). Each leaf value is a named integer vector: k-mer → count.

Usage

all.abundance

Format

An object of class list of length 3.

Details

Used to convert counts to density (mutations per megabase of context). See the planned transform_catalog() function (not yet ported).

Examples

all.abundance$BSgenome.Hsapiens.UCSC.hg38$transcript$`144`[1:5]

Description

Take an annotated indel VCF data frame (with columns Koh_476 and R as produced by indel classification functions) and produce a single-column data frame of mutation counts in the 476-category Koh classification scheme.

Usage

annot_vcf_to_476_catalog(
  annot_vcf,
  sample_id = "no_sample_id_provided",
  FILTER_PASS = TRUE,
  do_message = FALSE,
  clip_le_9 = TRUE
)

Arguments

Details

The function:

1. Optionally filters to PASS variants.

2. Removes duplicate positions (warns if ALT alleles differ).

3. Collapses single-base indels with repeat count $\ge 9$ into an "R(9,)" bin.

4. Tallies counts per Koh 476 category and returns a data frame with one row per category (using catalog_row_order()$ID476).

Value

A single-column data frame with 476 rows (one per Koh category) and integer mutation counts. Row names are the Koh 476 category strings; the column name is sample_id.

Description

Take an annotated indel VCF data frame (with column COSMIC_83 as produced by indel classification functions) and produce a single-column data frame of mutation counts in the 83-category COSMIC ID classification scheme.

Usage

annot_vcf_to_83_catalog(
  annot_vcf,
  sample_id = "no_sample_id_provided",
  FILTER_PASS = TRUE,
  do_message = FALSE,
  clip_le_9 = TRUE
)

Arguments

Details

The function:

1. Optionally filters to PASS variants.

2. Removes duplicate positions (warns if ALT alleles differ).

3. Tallies counts per COSMIC 83 category and returns a data frame with one row per category (using catalog_row_order()$ID).

Value

A single-column data frame with 83 rows (one per COSMIC ID category) and integer mutation counts. Row names are the COSMIC 83 category strings; the column name is sample_id.

Description

Take an annotated indel VCF data frame (with column Koh_89 as produced by indel classification functions) and produce a single-column data frame of mutation counts in the 89-category Koh classification scheme.

Usage

annot_vcf_to_89_catalog(
  annot_vcf,
  sample_id = "no_sample_id_provided",
  FILTER_PASS = TRUE,
  do_message = FALSE,
  clip_le_9 = TRUE
)

Arguments

Details

The function:

1. Optionally filters to PASS variants.

2. Removes duplicate positions (warns if ALT alleles differ).

3. Tallies counts per Koh 89 category and returns a data frame with one row per category (using catalog_row_order()$ID89).

Value

A single-column data frame with 89 rows (one per Koh category) and integer mutation counts. Row names are the Koh 89 category strings; the column name is sample_id.

Description

Port of ICAMS's AnnotateIDVCF(). Justifies each indel against the reference genome, optionally annotates transcript strand, and categorizes each justified indel into the COSMIC 83, Koh 89, and Koh 476 schemes.

Usage

annotate_id_vcf(
  vcf,
  ref_genome,
  trans_ranges = NULL,
  name_of_vcf = NULL,
  suppress_discarded_variants_warnings = TRUE,
  explain_indels = 1,
  context_width_multiplier = 20L,
  add_transcript_ranges = TRUE
)

Arguments

Value

A list with:

• annotated.vcf: the input VCF with new columns seq.context, seq.context.width, pos_shift, COSMIC_83, Koh_89, Koh_476, and (when transcript ranges are available) trans.strand / bothstrand.

• discarded.variants: rows that could not be justified, or NULL if none were discarded.

Description

Adds flanking sequence context via add_seq_context() and, when a transcript-ranges table is available, transcript strand via add_transcript_strand(). The same pipeline is appropriate for SBS and DBS; DBS-specific validation (e.g. no N in the tetranucleotide context) happens downstream in vcf_to_dbs_catalog().

Usage

annotate_sbs_or_dbs_vcf(
  vcf,
  ref_genome,
  trans_ranges = NULL,
  seq_context_width = 10L,
  name_of_vcf = NULL
)

Arguments

Value

A list with:

• annotated.vcf: the input VCF with new columns ⁠seq.<N>bases⁠ and (when transcript ranges are available) trans.strand / bothstrand.

• discarded.variants: currently always NULL for this path (included for symmetry with annotate_id_vcf()).

Description

Attaches the standard catalog attributes (type, counts_or_density, ref_genome, region, abundance) to x and returns it. No S3 class is set — mSigSpectra catalogs are plain matrices with attributes; functions key off attr(x, "type") via explicit checks.

Usage

as_catalog(
  x,
  type = NULL,
  ref_genome = NULL,
  region = "unknown",
  abundance = NULL,
  counts_or_density = "counts",
  infer_rownames = FALSE
)

Arguments

Value

x as a numeric matrix with attributes set.

Examples

m <- matrix(
  1, nrow = 96, ncol = 1,
  dimnames = list(catalog_row_order()$SBS96, "sample1")
)
cat96 <- as_catalog(m)
attr(cat96, "type")   # "SBS96"
attr(cat96, "region") # "unknown"

Description

Report the attributes of an mSigSpectra catalog

Usage

catalog_attrs(x)

Arguments

Value

A named list with elements type, counts_or_density, ref_genome, region, abundance.

Description

Returns a named list containing the canonical row ordering for each catalog type. These are used for validation and ordering of mutation catalogs.

Usage

catalog_row_order()

Details

Row names use a compact 4-letter format for SBS types: e.g. ACAA encodes the trinucleotide context as ⁠<5' base><ref><3' base><alt>⁠. SBS288 row names add a strand prefix: T:ACAA (transcribed), U:ACAA (untranscribed), N:ACAA (non-transcribed).

Value

A named list with elements: SBS96, SBS192, SBS288, SBS1536, DBS78, DBS136, DBS144, ID (83-category COSMIC indels), ID166, ID89, ID476.

Examples

cro <- catalog_row_order()
head(cro$SBS96)
length(cro$DBS78)

Description

This function is primarily for internal use, but we export it to document the underlying logic.

Usage

categorize_1_justified_indel(
  context,
  ins_or_del,
  ins_or_del_seq,
  pos,
  chrom = NULL,
  genomic_pos = NULL
)

Arguments

Details

See https://github.com/steverozen/ICAMS/blob/v3.0.9-branch/data-raw/PCAWG7_indel_classification_2021_09_03.xlsx for additional information on deletion mutation classification.

This function first handles deletions in homopolymers, then handles deletions in simple repeats with longer repeat units (e.g. CACACACA), and if the deletion is not in a simple repeat, looks for microhomology.

Value

A string that is the canonical representation of the given deletion type. Return NA and raise a warning if there is an un-normalized representation of the deletion of a repeat unit. See FindDelMH for details. (This seems to be very rare.)

Examples

simplify = function(ll) unlist(ll[c("COSMIC_83", "Koh_89", "Koh_476")])
categorize_1_justified_indel("GGAAAGG", "d", ins_or_del_seq = "A", pos = 3) # "DEL:T:1:2"
categorize_1_justified_indel("GGAAAGG", "d", ins_or_del_seq = "A", pos = 4) # "DEL:T:1:2"
simplify(
  categorize_1_justified_indel("TTATT", "d", ins_or_del_seq = "A", pos = 3))
simplify(
  categorize_1_justified_indel("TTATATAT", "d", ins_or_del_seq = "TATA", pos = 2))
simplify(
  categorize_1_justified_indel("TTATATAT", "d", ins_or_del_seq = "TATAT", pos = 2))

Description

Checks that all input catalogs share the same type, counts_or_density, ref_genome, and region, then cbinds their matrices and re-applies the shared attributes.

Usage

cbind_catalogs(catalogs)

Arguments

Value

A single catalog with ncol equal to the sum of the input ncols.

Description

Replaces bounded repeat-count suffixes like ⁠:R(X,9)⁠ with right-open equivalents like ⁠:R(X,)⁠, making the classification agnostic as to whether repeats longer than 9 were discarded from the input data.

Usage

change_476_type_ids_to_open_intervals(type_476_indel_type_identifiers)

Arguments

Value

Character vector the same length as the input, with ⁠:R(X,9)⁠ at the end of each string replaced by ⁠:R(X,)⁠.

Examples

change_476_type_ids_to_open_intervals(
  c("Del(C):Ins(C):R(5,9)", "Del(T):R(3,5)", "Ins(C):R(5,9)")
)

Description

Replaces bounded repeat-count suffixes like R(5,9) with right-open equivalents like R(5,), making the classification agnostic as to whether repeats longer than 9 were discarded from the input data.

Usage

change_89_type_ids_to_open_intervals(type_89_indel_type_identifiers)

Arguments

Value

Character vector the same length as the input, with bounded repeat suffixes replaced by right-open equivalents.

Examples

change_89_type_ids_to_open_intervals(
  c("Ins(2,):R(5,9)", "Ins(C):R(7,9)", "[Del(T):R(8,9)]",
    "[Ins(T):R(8,9)]", "Del(2,):U(1,2):R(5,9)", "Del(3,):U(3,):R(3,9)")
)

Description

Removes:

• Rows with identical REF and ALT.

• Stray ⁠#CHROM⁠ header repeats.

• Duplicated (CHROM, POS, REF, ALT) rows (keeping one copy).

• Multiple-ALT rows (comma-separated ALT field).

• Non-standard chromosome names (or those outside chr_names_to_process when supplied).

• Substitutions of length > 2 (e.g. ACT>TGA).

• Complex indels (REF[1] != ALT[1]).

• Wrong DBS rows where REF and ALT share a base at the same position.

• Variants whose REF base is not in ⁠{A, C, G, T}⁠.

Usage

check_and_remove_discarded_variants(
  vcf,
  name_of_vcf = NULL,
  chr_names_to_process = NULL
)

Arguments

Value

A list with element df, the retained rows of vcf (same class as the input), and, only when at least one row was removed, element discarded.variants, a data.table of the removed rows with an added character column discarded.reason explaining why each row was discarded.

Description

Supports:

• SBS1536 -> SBS96: collapse pentanucleotide → trinucleotide contexts.

• SBS192 -> SBS96: sum across the two transcript strands.

• DBS144 -> DBS78: sum across the two transcript strands.

Usage

collapse_catalog(catalog, to = c("SBS96", "DBS78"))

Arguments

Value

A new catalog with the collapsed rows.

Description

If trans_ranges is supplied, it is returned unchanged. Otherwise, returns the shipped transcript-ranges table for the given reference genome (GRCh37 / GRCh38 / GRCm38).

Usage

infer_trans_ranges(ref_genome, trans_ranges = NULL)

Arguments

Value

A data.table of transcript ranges, or NULL if no shipped table is available and none was supplied.

Description

Returns TRUE if x is a numeric matrix with the five catalog attributes (type, counts_or_density, ref_genome, region, abundance) and canonical rownames for its type.

Usage

is_catalog(x)

Arguments

Value

A single logical value: TRUE if x is a numeric matrix carrying the catalog attributes with the canonical row names for its type, otherwise FALSE.

Description

Add sequence context and transcript information to an in-memory ID (insertion/deletion) VCF, and confirm that they match the given reference genome

Usage

justify_id_vcf(
  ID.vcf,
  ref.genome,
  name.of.VCF = NULL,
  suppress.discarded.variants.warnings = TRUE,
  explain_indels = 1,
  context_width_multiplier = 20L
)

Arguments

Value

A list of elements:

• annotated.vcf: The original VCF data frame in which POS, REF and ALT with new columns added to the input data frame:

  • seq.context: The sequence embedding the variant.

  • seq.context.width: The width of seq.context to the left of the indel

  • pos_shift

• discarded.variants: Non-NULL only if there are variants that were excluded from the analysis. See the added extra column discarded.reason for more details.

Examples

# See tests/testthat/test_indel_classification.R for an end-to-end example
# against a real Strelka ID VCF.

Description

short_string should be generated by a deletion in long-str at (1-based) position pos. If supplied, the expected_delta should be the deleted sequence. If expected_delta is not NULL, the function checks this and aborts if this is not the case.

Dually, pos is the 1-based position in short string in front o which one can make an insertion to get long_str.

This function moves pos as far to left as possible so that a deletion at that position still results in an edit of long_str to generate short_str.

Usage

justify_indel(long_str, short_str, pos, expected_delta = NULL)

Arguments

Value

A list with elements

• leftmost_pos

• del_string

• del_string_plus del_string plus one base to the left the position of that base to the left is leftmost_pos - 1

• edge_warning See the third example. We are not sure how far to the left we could have moved the indel if we had more context on the left.

• error if non-NULL we did not see the expected value for expected_delta

Examples

justify_indel("CAAAG", "CAAG", pos = 2, expected_delta = "A")
justify_indel("CAAAG", "CAAG", pos = 3, expected_delta = "A")
justify_indel("CACAG", "CAG", pos = 3, expected_delta = "CA")

Description

Auto-detects the file format from its first row. Recognized formats:

• ICAMS native CSV — first N columns are header columns (mutation type, context, etc), remaining columns are per-sample counts / densities.

• SigProfiler TSV / CSV — single header column with bracketed mutation labels like ⁠A[C>A]A⁠, ⁠AA[C>A]AA⁠, or PCAWG indel codes like 1:Del:C:0.

• COSMIC CSV — SBS96 uses the ICAMS-external row-header format; SBS192 uses a stranded variant of the same; DBS78 / ID83 also share their SigProfiler-style layout.

Usage

read_catalog(
  file,
  ref_genome = NULL,
  region = "unknown",
  counts_or_density = "counts",
  format = c("auto", "ICAMS", "SigProfiler", "COSMIC")
)

Arguments

Details

The matrix is reordered to the canonical rownames for its type and wrapped in a catalog via as_catalog().

Value

A catalog matrix with attributes (see as_catalog()).

Description

Reads the body of a VCF file (lines after ⁠#CHROM⁠) into a data.table. The resulting table has whatever columns the VCF has (CHROM, POS, ID, REF, ALT, QUAL, FILTER, INFO, and optionally FORMAT plus one or more sample columns), with ⁠#CHROM⁠ renamed to CHROM.

Usage

read_vcf(file, filter = TRUE, name_of_vcf = NULL)

Arguments

Details

Caller-agnostic. read_vcf() does not know or care which variant caller produced the VCF. It does not parse FORMAT/sample columns and does not extract VAF or read depth. The only caller-dependent semantics is the default value of the filter argument (see below).

Uses data.table::fread() with check.names=FALSE, ⁠na.strings = ""``, ⁠fill = TRUE', to parse the VCF body. Handles uncompressed and gzipped files; does not handle bgzipped/tabix.

Value

A data.table with one row per variant. The name of the first column is 'CHROM', not '#CHROM'. Other column names

Description

Thin batch wrapper around read_vcf().

Usage

read_vcfs(files, ..., names_of_vcfs = NULL)

Arguments

Value

A named list of data.tables, one per input file.

Description

Segment a single indel sequence using Rcpp interface

Usage

seg_simple(ins_or_del, string, context)

Arguments

Details

This function segments an indel sequence by finding the optimal repeat unit that best explains the sequence structure. The algorithm tries all possible repeat unit sizes and selects the best segmentation based on:

• Highest 3' flanking repeat count (most important)

• Highest internal repeat count

• Lowest spacer length (prefer clean repeats)

• Lowest unit length (prefer simpler units)

Value

A list with the following elements:

Examples

# Simple AT repeat
result <- seg_simple("d", "ATATAT", "ATATGG")
print(result)

# CG repeat
result <- seg_simple("d", "CGCGCG", "CGCGAA")
print(result)

Description

This function provides direct C++ interface for indel segmentation without calling an external binary via system2.

Usage

segment_simple_cpp(ins_or_del, string, context)

Arguments

Value

A list with the following elements:

Description

Classifies each row by REF/ALT length alone:

• SBS (single base substitution): nchar(REF) == 1 && nchar(ALT) == 1

• DBS (double base substitution): nchar(REF) == 2 && nchar(ALT) == 2

• ID (indel): nchar(REF) != nchar(ALT)

• Other (e.g. 3+bp substitutions): discarded with a recorded reason.

Usage

split_vcf(vcf, name_of_vcf = NULL)

Arguments

Details

This is caller-agnostic. In particular, mSigSpectra does not merge adjacent SBSs into DBSs based on VAF similarity (which ICAMS did via SplitOneVCF for Strelka-style VCFs). Users who want that behavior should apply it as a post-processing step with their own VAF column.

Value

A list with elements

• SBS: data.table of SBS rows.

• DBS: data.table of DBS rows.

• ID: data.table of indel rows.

• discarded: data.table of rows that did not fit any of the above, with a discarded.reason column — NULL when no rows were discarded.

Description

Base [ on a matrix drops attributes other than dim / dimnames. Use subset_catalog() when you want to keep the catalog's type / ref_genome / region / abundance metadata through a subset operation.

Usage

subset_catalog(x, rows = NULL, cols = NULL)

Arguments

Value

A catalog (numeric matrix) containing the selected rows and columns, with the type, counts_or_density, ref_genome, and region attributes of x preserved. The abundance attribute is preserved only when all rows are kept, because it is not meaningful for a subset of mutation types.

Description

Precomputed transcript ranges (one row per gene) used by add_transcript_strand() to determine the coding strand for each variant. Sources: GENCODE v30 (human) and vM21 (mouse). Only genes with CCDS IDs are retained.

Usage

trans.ranges.GRCh37

trans.ranges.GRCh38

trans.ranges.GRCm38

Format

A data.table::data.table with columns chrom, start, end, strand, Ensembl.gene.ID, gene.symbol. One-based coordinates.

An object of class data.table (inherits from data.frame) with 19083 rows and 6 columns.

An object of class data.table (inherits from data.frame) with 19096 rows and 6 columns.

An object of class data.table (inherits from data.frame) with 20325 rows and 6 columns.

Source

https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_30/GRCh37_mapping/gencode.v30lift37.annotation.gff3.gz

https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_30/gencode.v30.annotation.gff3.gz

https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M21/gencode.vM21.annotation.gff3.gz

Description

Converts between counts (raw mutation counts per category) and density (mutations per megabase of context) representations, and between catalog <-> signature (column-normalized) forms. The transformation is expressed by multiplying each category's count by target_abundance[source.n.mer] / source_abundance[source.n.mer], where source.n.mer is the reference context encoded in the row name (first 3 characters for SBS96/192, 5 for SBS1536, 2 for DBS78/144, 4 for DBS136).

Usage

transform_catalog(
  catalog,
  target_ref_genome = NULL,
  target_region = NULL,
  target_counts_or_density = NULL,
  target_abundance = NULL
)

Arguments

Value

A new catalog (plain matrix with attributes).

Description

Build a DBS mutational-spectrum catalog from an annotated DBS VCF

Usage

vcf_to_dbs_catalog(
  annotated_vcf,
  type = c("DBS78", "DBS136", "DBS144"),
  ref_genome = NULL,
  region = "unknown",
  sample_name = "count"
)

Arguments

Value

A single-column numeric matrix with catalog attributes (see as_catalog()).

Description

Turns an indel-annotated VCF (with COSMIC_83 / Koh_89 / Koh_476 columns as produced by annotate_id_vcf()) into a count matrix for the requested ID classification scheme.

Usage

vcf_to_id_catalog(
  annotated_vcf,
  type = c("ID83", "ID89", "ID476"),
  ref_genome = NULL,
  region = "unknown",
  sample_name = "count",
  FILTER_PASS = TRUE,
  clip_le_9 = TRUE
)

Arguments

Value

A single-column numeric matrix with catalog attributes (see as_catalog()).

Description

Returns a single catalog matrix of the requested type. Intermediate matrices for the other SBS resolutions are still computed (cheap) but not returned, keeping the public API focused on "one call, one catalog type".

Usage

vcf_to_sbs_catalog(
  annotated_vcf,
  type = c("SBS96", "SBS192", "SBS1536"),
  ref_genome = NULL,
  region = "unknown",
  sample_name = "count"
)

Arguments

Value

A single-column numeric matrix with catalog attributes (see as_catalog()).

Description

Writes catalog in ICAMS-native CSV format (the only format currently supported for writing). The row-header columns that precede the sample columns are taken from the shipped catalog.row.headers object for the catalog's type.

Usage

write_catalog(
  catalog,
  file,
  format = c("ICAMS", "SigProfiler", "COSMIC"),
  sep = ","
)

Arguments

Value

file invisibly.