Package {mSigSpectra}

mSigSpectra: Build mutational-spectrum catalogs from VCF files

Description

Reads variant call files (VCFs) caller-agnostically, annotates variants with flanking sequence context and transcriptional strand, and builds mutational- spectrum catalogs (SBS, DBS, indel) at multiple resolutions and in both counts and density representations.

Details

Catalogs are plain numeric matrices (rows = mutation categories, columns = samples) carrying attributes type, ref_genome, region, counts_or_density, and optionally abundance. mSigSpectra deliberately does not use S3 classes or method dispatch on catalogs; behavior is selected by explicit attr(x, "type") checks inside regular functions.

Plotting is intentionally out of scope and is handled by the separate mSigPlot package.

Author(s)

Maintainer: Steve Rozen steverozen@pm.me

Authors:

• Steve Rozen steverozen@pm.me

• Nanhai Jiang

• Arnoud Boot

• Mo Liu

• Yang Wu

• Mi Ni Huang

• Jia Geng Chang

See Also

Useful links:

• https://github.com/steverozen/mSigSpectra

• https://steverozen.github.io/mSigSpectra/

• Report bugs at https://github.com/steverozen/mSigSpectra/issues

Add flanking sequence context to a VCF data frame

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⁠.

Annotate a VCF data frame with transcript strand information

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).

K-mer abundances for density calculations

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]

Convert an annotated indel VCF to a Koh 476-category catalog

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.

Convert an annotated indel VCF to a COSMIC 83-category catalog

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.

Convert an annotated indel VCF to a Koh 89-category catalog

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.

Annotate an in-memory ID (indel) VCF with sequence context, transcript strand, and COSMIC / Koh indel categories

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.

Annotate an SBS or DBS VCF with flanking sequence context and transcript strand

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()).

Turn a numeric matrix into a mutational-spectrum catalog

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"

Report the attributes of an mSigSpectra catalog

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.

Return catalog row orders for all supported catalog types

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)

Given a indel and its sequence context, categorize it

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

Combine catalogs across samples (column-bind)

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.

Change 476-type indel category identifiers to use right-open repeat intervals

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

Change 89-type indel category identifiers to use right-open repeat intervals

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

Check and, if possible, correct the chromosome names in a VCF data.frame

Description

Harmonizes the CHROM column of vcf.df with the chromosome names used by ref.genome (a BSgenome object). Adds or strips the chr prefix as needed; for human and mouse, translates 23/24 or 20/21 to X/Y.

Usage

check_and_fix_chrom_names(vcf.df, ref.genome, name.of.VCF = NULL)

Arguments

Value

A character vector of chromosome names that can be used as a replacement for vcf.df$CHROM. Errors if reconciliation is impossible.

Check and, if possible, correct the chromosome names in a trans.ranges table

Description

Harmonizes trans.ranges$chrom with the chromosome naming used by vcf.df$CHROM, adding or stripping chr as needed and translating organism-specific numeric X/Y encodings.

Usage

check_and_fix_chrom_names_for_trans_ranges(
  trans.ranges,
  vcf.df,
  ref.genome,
  name.of.VCF = NULL
)

Arguments

Value

A character vector of chromosome names that can be used as a replacement for trans.ranges$chrom. Errors if reconciliation is impossible.

Check a VCF for common variant-level problems and remove the offenders

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.

Reorder catalog rows to the canonical order for its type

Description

Reorder catalog rows to the canonical order for its type

Usage

check_and_reorder_rownames(x, type)

Arguments

Value

x with rows reordered to match catalog_row_order()[[type]]. Errors if any canonical rowname is missing.

Collapse a higher-resolution catalog to a lower-resolution one

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.

Infer k-mer abundance from attributes

Description

Looks up the appropriate entry of all.abundance keyed on ref_genome and region. Returns NULL if no entry is found, if the catalog is a signature / density catalog with a flat abundance, or if the catalog is COMPOSITE (which has no meaningful abundance).

Usage

infer_abundance(x, ref_genome, region, counts_or_density)

Value

A named integer vector of k-mer counts (the abundance for the catalog's context size), or NULL if no abundance applies.

Infer catalog type from the number of rows of a matrix

Description

Infer catalog type from the number of rows of a matrix

Usage

infer_catalog_type(n_rows)

Arguments

Value

A type identifier (e.g. "SBS96", "DBS78", "ID83"). Errors if n_rows does not correspond to a supported catalog type.

Map a character ref_genome argument to its canonical BSgenome package name

Description

Map a character ref_genome argument to its canonical BSgenome package name

Usage

infer_ref_genome_name(ref_genome)

Value

A single character string giving the canonical 'BSgenome' package name, e.g. "BSgenome.Hsapiens.UCSC.hg38". Errors if ref_genome is not recognized.

Infer transcript ranges for a reference genome

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.

Check whether an object looks like an mSigSpectra catalog

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.

Is this reference genome GRCh37 (1000 Genomes hs37d5)?

Description

Is this reference genome GRCh37 (1000 Genomes hs37d5)?

Usage

is_grch37(x)

Arguments

Value

A single logical value, TRUE if x identifies that genome, otherwise FALSE.

Is this reference genome GRCh38 (UCSC hg38)?

Description

Is this reference genome GRCh38 (UCSC hg38)?

Usage

is_grch38(x)

Arguments

Value

A single logical value, TRUE if x identifies that genome, otherwise FALSE.

Is this reference genome GRCm38 (UCSC mm10)?

Description

Is this reference genome GRCm38 (UCSC mm10)?

Usage

is_grcm38(x)

Arguments

Value

A single logical value, TRUE if x identifies that genome, otherwise FALSE.

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

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.

Move the notional position of a deletion or insertion as far left as possible.

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

Justify indels in a VCF and update positions accordingly

Description

For each indel in the input VCF, this function:

1. Justifies the indel position (moves it as far left as possible)

2. Calculates how much the position was shifted

3. Updates the POS column by decrementing it by the shift amount

4. Updates the seq.context.width column accordingly

5. Adds a new column 'pos_shift' showing how much each position was moved

6. Categorizes the justified indel

Usage

justify_indels_in_id_vcf_with_contexts(vcf, explain_indels = 1)

Arguments

Value

A data.frame with:

• All original columns from the input vcf

• Updated POS values (decremented by pos_shift)

• Updated REF and ALT values if the position of the indel has been moved

• Updated seq.context.width values (decremented by pos_shift)

• New column 'pos_shift': Amount the position was moved (usually 0)

Important invariant

The seq.context string is never modified. Both POS and seq.context.width are decremented by the same amount (pos_shift), which maintains the relationship: genomic position POS corresponds to position (seq.context.width + 1) in seq.context.

For example, if seq.context was extracted from genomic positions (POS - seq.context.width, POS + var.width + seq.context.width), then after justification by shift S:

• seq.context remains unchanged

• POS becomes POS - S

• seq.context.width becomes seq.context.width - S

• The variant is now at position (seq.context.width - S) + 1 in seq.context

Compute length of longest common prefix of two strings

Description

Fast alternative to Biostrings::lcprefix that avoids S4 method dispatch overhead.

Usage

lcprefix_fast(a, b)

Arguments

Value

An integer: the number of leading characters that match.

Normalize a reference-genome argument to a BSgenome object

Description

Accepts either a BSgenome object directly, or one of the string identifiers "GRCh37" / "hg19" / "BSgenome.Hsapiens.1000genomes.hs37d5", "GRCh38" / "hg38" / "BSgenome.Hsapiens.UCSC.hg38", or "GRCm38" / "mm10" / "BSgenome.Mmusculus.UCSC.mm10". The relevant BSgenome package must be installed; an informative error is raised if not.

Usage

normalize_genome_arg(ref_genome)

Arguments

Value

A BSgenome object.

Normalize SBS1536 pentanucleotide + ALT strings to pyrimidine form

Description

Input strings are 6 characters: a 5-base pentanucleotide context followed by a 1-base ALT (e.g. "ATGCTT" = ATGCT>T). If the center of the pentanucleotide (position 3) is A or G, the pentanucleotide and the ALT are reverse-complemented so that the center becomes C or T — the pyrimidine form used in canonical SBS96/1536 row names.

Usage

pyr_penta(mutstring)

Arguments

Value

A character vector the same length as mutstring, with each element in pyrimidine-centered form.

Filter, deduplicate, and check an annotated indel VCF

Description

Helper used by annot_vcf_to_83_catalog, annot_vcf_to_89_catalog, and annot_vcf_to_476_catalog to perform shared preprocessing: rename #CHROM to CHROM, optionally filter to PASS variants, warn about positions with differing ALT alleles, and deduplicate by position.

Usage

quick_check_vcf(annot_vcf, FILTER_PASS = FALSE, do_message = FALSE)

Arguments

Value

A data frame deduplicated by position, with a pos_id column added.

Read a mutational-spectrum catalog from a file

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()).

Read transcript ranges from a GENCODE-derived CSV

Description

Reads a CSV with columns chrom, start, end, strand, Ensembl.gene.ID, gene.symbol (1-based coordinates), orders the chromosome factor canonically (1..22/19, X, Y), and returns a keyed data.table.

Usage

read_transcript_ranges(file)

Arguments

Value

A data.table keyed on chrom, start, end.

Read a VCF file into a data.table, caller-agnostically

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

Read multiple VCF files

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.

Remove variants with duplicated CHROM+POS

Description

Two passes:

1. Rows with identical (CHROM, POS, REF, ALT): keep one copy, discard the rest.

2. Rows sharing (CHROM, POS, REF) but different ALT: discard all (treated as unresolved multiallelic / inconsistent records).

Usage

remove_rows_with_duplicated_chrom_and_pos(df, name_of_vcf = NULL)

Value

A list with element df (the retained rows) and, only when rows were removed, element discarded.variants (the removed rows with an added character column discarded.reason).

Remove rows that look like stray "#CHROM" header repeats

Description

Occasionally a VCF body contains a spurious row whose CHROM column equals "#CHROM" (e.g. when a concatenated VCF keeps the header line of each input). Drop those rows and record them in discarded.variants.

Usage

remove_rows_with_pound_sign(df, name_of_vcf = NULL)

Value

A list with element df (the retained rows) and, only when rows were removed, element discarded.variants (the removed rows with an added character column discarded.reason).

Reverse complement strings that represent stranded DBSs

Description

Input is a 4-character string where characters 1-2 are the REF dinucleotide and characters 3-4 are the ALT (e.g. "AATC" = AA>TC). Returns the reverse complement of the first 2 characters concatenated with the reverse complement of the last 2 characters, e.g. "AATC" returns "TTGA".

Usage

revc_dbs144(mutstring)

Arguments

Value

A character vector the same length as mutstring containing the reverse-complemented strings.

Reverse complement strings that represent stranded SBSs

Description

Input is a 4-character string where characters 1-3 are the trinucleotide context and character 4 is the ALT (e.g. "AATC" = AAT>ACT). Returns the reverse complement of the first 3 characters concatenated with the reverse complement of the last character, e.g. "AATC" returns "ATTG".

Usage

revc_sbs96(mutstring)

Arguments

Value

A character vector the same length as mutstring containing the reverse-complemented strings.

Segment a single indel sequence using Rcpp interface

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)

Segment a single indel using Rcpp interface

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:

Restrict a VCF data frame to a user-specified set of chromosome names

Description

Restrict a VCF data frame to a user-specified set of chromosome names

Usage

select_variants_by_chrom_name(df, chr.names.to.process, name.of.VCF = NULL)

Arguments

Value

A list with elements

• df: data frame of rows whose CHROM is in chr.names.to.process.

• discarded.variants: non-NULL only if any rows were discarded.

Split a mixed-mutation VCF into SBS / DBS / ID sub-tables

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.

Standardize chromosome names in the first column of a data frame

Description

Drops rows whose chromosome name contains any of GL, KI, random, Hs (anchored at start), M, or JH, and strips any leading chr prefix from the remainder.

Usage

standard_chrom_name(df)

Arguments

Value

df restricted to rows with canonical chromosome names (1:22, X, Y), with chr prefixes removed.

Standardize the chromosome names in a VCF data.frame

Description

Splits df into rows with canonical chromosome names and rows with non-standard chromosome names (those containing any of GL, KI, random, Hs, M, JH, fix, alt). Emits a warning when any rows are discarded.

Usage

standard_chrom_name_new(df, name.of.VCF = NULL)

Arguments

Value

A list with elements

• df: data frame of rows with canonical chromosome names.

• discarded.variants: non-NULL only if any rows were discarded; each discarded row gets a discarded.reason column.

Validate a counts_or_density argument

Description

Validate a counts_or_density argument

Usage

stop_if_counts_or_density_illegal(counts_or_density)

Value

NULL, invisibly. Called for its side effect of raising an error when counts_or_density is not a legal value.

Validate a region argument

Description

Validate a region argument

Usage

stop_if_region_illegal(region)

Arguments

Value

NULL invisibly; raises an error if region is not one of "genome", "exome", "transcript", "unknown".

Validate a region argument for catalog types that require transcript strand

Description

SBS192, DBS144 and similar stranded catalogs cannot be built from region = "genome", since variants outside transcripts have no strand.

Usage

stop_if_transcribed_region_illegal(region)

Arguments

Value

NULL, invisibly. Called for its side effect of raising an error when region is not legal for a stranded catalog.

Subset a catalog while preserving attributes

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.

Transcript ranges for transcriptional strand annotation

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

Transform a catalog between counts and density

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).

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

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()).

Build an ID (indel) mutational-spectrum catalog from an annotated ID VCF

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()).

Build an SBS mutational-spectrum catalog from an annotated SBS VCF

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()).

Write a mutational-spectrum catalog to a file

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.