Package 'MutSeqR'

Description

A simple method to test whether your trinucleotide context contains a CpG site. Vectorized version of Biostrings::vcountPattern is used.

Usage

annotate_cpg_sites(mutation_data, motif = "CG", column_query = "context", ...)

Arguments

Value

A data frame with the same number of rows as there were ranges in the input, but with an additional metadata column indicating CpG sites in the target sequence of the mutation.

Description

Calculate the benchmark dose (BMD) of continuous, individual-level data with optional model averaging. This function is intended to model the dose-response of mutation frequency. This function is an extension of the PROAST software (copyright RIVM National Institute for Public Health and the Environment).

Usage

bmd_proast(
  mf_data,
  dose_col = "dose",
  response_col = "mf_min",
  covariate_col = NULL,
  bmr = 0.5,
  adjust_bmr_to_group_sd = FALSE,
  model_averaging = TRUE,
  num_bootstraps = 200,
  plot_results = FALSE,
  output_path = NULL,
  raw_results = FALSE,
  seed = 125
)

Arguments

Details

This function is a modified version of the original interactive PROAST software (https://www.rivm.nl/en/proast that allows for batch processing of data. The function is designed to be used with the output of calculate_mf for the purpose of calculating the Benchmark Dose of mutation frequency data. As such, some functionality of the original PROAST software has been removed.

This function will accept continuous data, with an observation for each individual subject. It is assumed that data are lognormally distributed. The response data is log-transformed, then back-transformed after the statistical analysis. The function will fit model 3 or 5 from various families of models (Exponential, Hill, Inverse Exponential, LogNormal). It will then compare the fits of models 3 and 5 for each model family and select the model with the lowest AIC. The BMD 90% confidence intervals will be calculated based on the selected model (3 or 5) for each model family using the profile likelihood method. The BMD 90% confidence interval may also be calculated using the bootstrap method if model_averaging = TRUE. It is recommended to use 200 bootstraps for model averaging.

To replicate these results in the PROAST interactive software, select the following menu options:

1. f.proast(mf_data)

2. What type of response data do you want to consider? 1: continuous, individual data

3. Do you want to fit a single model or fit various nested families of models? 3: select model 3 or 5 from various families of models

4. Q1: Which variable do you want to consider as the independent variable? # : dose_col

5. Give number(s) of the response(s) you want to analyse. # : response_col

6. Give number of factor serving as potential covariate (e.g.sex) type 0 if none. # : covariate_col

7. Do you want to adjust CES to within group SD? 1: no, 2: yes | adjust_bmr_to_group_sd: FALSE/TRUE

8. Give value for CES (always positive) type 0 to avoid calculation of CIs. bmr

9. Do you want to calculate the BMD confidence interval by model averaging? 1: no 2: yes | model_averaging: FALSE/TRUE

10. give number of bootstrap runs for calculating BMD confidence interval based on MA (e.g. 200) num_bootstraps

11. Which models do you want to be fitted? 4 : previous option with lognormal DR model added

Value

A summmary data frame of final results. If plots or raw results are selected, all data will be returned within a list.

The summary will include the following for each response variable and covariate subgroup (if applicable):

• Model: The m3 or m5 model selected for each model family (Exponential, Hill, Inverse Exponential, LogNormal).

• Response: The response variable.

• Covariate: The covariate subgroup, if applicable.

• bmr: The specified Benchmark Response.

• BMD: The Benchmark Dose, in original dose units, estimated for the given model.

• BMDL: The lower bound of the 90% confidence interval for the BMD, calculated by the profile likelihood method.

• BMDU: The upper bound of the 90% confidence interval for the BMD, calculated by the profile likelihood method.

• AIC: The Akaike Information Criterion for the selected model. Lower values indicate a better fit. It is advised to choose the BMD value from the model with the lowest AIC.

• weights: The weight of the model in the model averaging process, if applicable.

• ⁠Model averaging⁠: The BMDL and BMDU calculated by the bootstrap method if model_averaging = TRUE.

If there is no significant response in the data, the function will return an empty data frame.

If plot_results = TRUE the function will create the following plots for each response variable. The plots will be saved to the output_path. If no output_path is provided, then they will be returned within a list alongside the summary data frame.

• Model Plots. The following plot will be created for each model family (Exponential, Hill, Inverse Exponential, LogNormal): The fitted curve of the selected (3 or 5) model. Data is log-transformed. Individual data points are plotted using small triangles. The geometric mean (median) at each dose is plotted as a large triangle. The BMD is indicated by the dotted line. If applicable, the covariate subgroup is indicated by color. When output_path = NULL, plots are returned alongside results as recordedPlot objects. View plots with replayPlot().

• bootstrap_curves If model_averaging = TRUE, the bootstrap curves based on model averaging. The geometric mean (median) at each dose is plotted as a large triangle. Data is log-transformed. When output_path = NULL, plots are returned alongside results as recordedPlot objects. View plots with replayPlot().

• cleveland plot if model_averaging = TRUE The BMD estimate for each model is plotted as a red point alongside the 90% confidence intervals. The size of the BMD point represents the model weight assigned during model averaging, based on the AIC. When output_path = NULL, plots are returned alongside results as ggplots. View plots by calling object name.

If raw_results = TRUE, the function will return the raw results of the PROAST analysis alongside the summary data frame. PROAST raw_results is a list of variables and data that is continuously modified as it is passed through the proast functions. It can be given to f.proast() to resume analysis.

Examples

# Example data  consists of 24 mouse bone marrow
# samples exposed to three doses of BaP alongside vehicle controls.
# Libraries were sequenced with Duplex Sequencing using
# the TwinStrand Mouse Mutagenesis Panel which consists of 20 2.4kb
# targets = 48kb of sequence. Example data can be retrieved from
# MutSeqRData, an ExperimentHub data package:
## library(ExperimentHub)
## eh <- ExperimentHub()
## query(eh, "MutSeqRData")
# Mutation frequency data was precalculated using
## mf_data_global <- calculate_mf(mutation_data = eh[["EH9861"]],
##   cols_to_group = "sample",
##   retain_metadata_cols = c("dose_group", "dose"))
mf_example <- readRDS(system.file("extdata/Example_files/mf_data_global.rds",
  package = "MutSeqR"
))
# We will calculate the BMD for a 50% increase in mutation frequency from
# control with Model averaging.
# For the purpose of this example, num_bootstraps is set to 3 to reduce
# run time. 200 bootstraps is recommended.
bmd <- bmd_proast(
  mf_data = mf_example,
  dose_col = "dose",
  response_col = "mf_min",
  bmr = 0.5,
  model_averaging = TRUE,
  num_bootstraps = 3
)

Description

A dictionary to map common organism names to their corresponding BSgenome organism names.

Usage

BS_org_map

Format

A list with organism names as keys and corresponding BSgenome organism names as values.

Examples

BS_org_map["Hsapiens"]

Description

Calculates mutation frequencies for arbitrary groupings and creates a new dataframe with the results. Mutation frequency is calculated by dividing the sum of mutations by the sum of the total_depth for a given group (mutations/bp). The operation is run using both the minimum and maximum independent mutation counting methods.

Usage

calculate_mf(
  mutation_data,
  cols_to_group = "sample",
  subtype_resolution = "none",
  variant_types = c("snv", "deletion", "insertion", "complex", "mnv", "sv", "ambiguous",
    "uncategorized"),
  calculate_depth = TRUE,
  correct_depth = TRUE,
  correct_depth_by_indel_priority = FALSE,
  precalc_depth_data = NULL,
  d_sep = "\t",
  summary = TRUE,
  retain_metadata_cols = NULL
)

Arguments

Details

Required columns:

• contig: (or seqnames) The reference sequence name.

• start: 1-based start position of the feature.

• alt_depth: The read depth supporting the alternate allele.

• variation_type: The category to which this variant is assigned.

• subtype_col: The column containing the mutation subtype. This column depends on the subtype_resolution parameter.

• reference context: The column containing the referene base(s) for the mutation. This column depends on the subtype_resolution parameter.

• cols to group: all columns across which you want to calculate the mutation frequency. Ex. c("tissue", "dose"). These columns should be listed in cols_to_group.

It is also required to include the total_depth column if you are calculating depth from the mutation data. If you are using precalculated depth data, the total_depth column is not required.

Subtype Resolutions:

• "none" calculates mutation frequencies across all selected grouping columns.

• "type" calculates mutation frequencies across all selected grouping columns for each variation_type seperately; snv, mnv, deletion, insertion, complex, sv, ambiguous, uncategorized.

• "base_6" calculates mutation frequencies across all selected grouping columns for each variation_type with snv mutations separated by normalized_subtype; C>A, C>G, C>T, T>A, T>C, T>G. The reference context is normalized_ref.

• "base_12" calculates mutation frequencies across all selected grouping columns for each variation_type with snv mutations separated by subtype; A>C, A>G, A>T, C>A, C>G, C>T, G>A, G>C, G>T, T>A, T>C, T>G. The reference context is short_ref.

• "base_96" calculates mutation frequencies across all selected grouping columns for each variation_type with snv mutations separated by normalized_context_with_mutation, i.e. the 96-base trinucleotide context. Ex. A[C>T]A. The reference context is normalized_context.

• "base_192" calculates mutation frequencies across all selected grouping columns for each variation_type with snv mutations separated by context_with_mutation, i.e. the 192-base trinucleotide context. Ex A[G>A]A. The reference context is context.

Subtype depth: For SNV subtypes, the total_depth is summed based on the sequence context in which the SNV subtype occurs. Ex. for base_6, the two possible reference bases are C or T; hence, the total_depth is summed seperately for C:G positions and T:A positions. The MF for C>T mutations is calculated as total # C>T mutations / total_depth for C>G positions (sum / subtype_depth). Non-SNV mutation types will be caluclated as their sum / group_depth, since they can occur in the context of any nucleotide.

retain_metadata_cols at subtype_resolution: The summary table uses a pre-defined list of possible subtypes for each resolution. If a particular subtype within a given group is not recorded in the mutation data, the summary table will have no frame of reference for populating the metadata_cols. Thus, for subtypes that do not occur in the mutation data for a given group, the corresponding metadata_col will be NA.

Variant filtering: Variants flagged as TRUE in the filter_mut column will be excluded from the mutation counts. However, the total_depth of these variants will be included in the group/subtype depths if calculating depth.

Depth correction is important for preventing double-counting of reads in mutation data when summing the total_depth across samples or other groups. Generally, when several mutations have been detected at the same genomic position, within a sample, the total_depth value will be the same for all of them. However, in some datasets, whenever a deletion is detected, the data may contain an additional row with the same genomic position calling a "no_variant". The total_depth will differ between the deletion and the no_variant. In these cases, correct_depth_by_indel_priority == TRUE will ensure that the total_depth value for the deletion is retained, while the total_depth value for the no_variant is removed.

Value

A data frame with the mutation frequency calculated. If summary is set to TRUE, the data frame will be a summary table with the mutation frequency calculated for each group. If summary is set to FALSE, the mutation frequency will be appended to each row of the original mutation_data.

• sum_min: The sum of all mutations within the group, calculated using the "min" method for mutation counting. All identical mutations within a samples are assumed to be the result of clonal expansion and are thus only counted once.

• sum_max: The sum of all mutations within the group, calculated using the "max" method for mutaiton counting. All identical mutations within a sample are assumed to be idenpendant mutational evens and are included in the mutation frequency calculation.

• group_depth: The total_depth summed across groups.

• subtype_depth: The total_depth summed across groups for a given sequence context. Used for calculating subtype frequencies.

• mf_min: The mutation frequency calculated using the "min" method for mutation counting. mf_min = sum_min / depth.

• mf_max: The mutation frequency calculated using the "max" method for mutation counting. mf_max = sum_max / depth.

• proportion_min: The proportion of each mutation subtype within the group, normalized to the depth. Calculated using the "min" method. This is only calculated if subtype_resolution is not "none". If no depth is calculated or provided, proportion is calculated without normalization to the depth.

• proportion_max: The proportion of each mutation subtype within the group, normalized to its read depth. Calculated using the "max" method. This is only calculated if subtype_resolution is not "none". If no depth is calculated or provided, proportion is calculated without normalization to the depth.

Examples

# Mutation data is just for example purposes. It does not reflect real data
mutation_data <- readRDS(system.file("extdata", "Example_files",
                                     "filtered_simple_mutation_data.rds",
                                     package = "MutSeqR"))
# Calculate mutation frequency by sample.
# Calculate depth from the mutation data (default)
# Correct the Depth (default) with indel priority (set)
mf_example <- calculate_mf(
  mutation_data = mutation_data,
  cols_to_group = "sample",
  correct_depth_by_indel_priority = TRUE
)

Description

This function generates additional columns for the mutation data, including a breakdown of the mutation subtypes at various resolutions.

Usage

characterize_variants(mutation_data)

Arguments

Value

A data frame with additional columns for variant characterization.

Description

A utility function that will check that all required columns are present.

Usage

check_required_columns(data, required_columns)

Arguments

Value

an error

Examples

df <- data.frame(
  contig = c("chr1", "chr2", "chr3"),
  start = c(100, 200, 300),
  end = c(100, 200, 300),
  sample = c("S1", "S2", "S3"),
  ref = c("G", "C", "T"),
  alt = c("A", "T", "G")
)
check_required_columns(df, required_columns = op$base_required_mut_cols)

Description

Classify the variation type of a mutation based on its ref and alt values.

Usage

classify_variation(ref, alt)

Arguments

Value

A character indicating the type of variation.

Examples

df <- data.frame(
  ref = c("A", "CAGT", "GCC", "T", "ACG", "C", "G", "T", "A"),
  alt = c("R", "TGA", "G", "TC", "TAC", "C", "<DEL>", "G", "???")
)
df$variation_type <- mapply(classify_variation, df$ref, df$alt)
df

Description

Make a Cleveland plot for the PROAST results.

Usage

cleveland_plot(results, covariate_col = NULL, output_path = NULL)

Arguments

Value

A single ggplot object with facets for each response.

Description

perform hierarchical clustering of samples based on the mutation spectra.

Usage

cluster_spectra(
  mf_data,
  group_col = "sample",
  response_col = "proportion_min",
  subtype_col = "normalized_subtype",
  dist = "cosine",
  cluster_method = "ward.D"
)

Arguments

Details

The cosine distance measure represents the inverted cosine similarity between samples:

$\text{Cosine Dissimilarity} = 1 - \frac{\mathbf{A} \cdot \mathbf{B}}{\| \mathbf{A} \| \cdot \| \mathbf{B} \|}$

This equation calculates the cosine dissimilarity between two vectors A and B.

Leaves are sorted using dendsort, if installed, otherwise leaves are unsorted.

Value

A dendrogram object representing the hierarchical clustering of the samples.

Description

A list of reference contexts at different resolutions

Usage

context_list

Format

A list with corresponding values

Examples

context_list[["base_6"]]
context_list[["base_12"]]

Description

These values are used to cross reference base substitution types to their appropriate denominators for calculations. That is", "for example, the 6 base substitution frequency should be subsetted based on the normalized_ref column which would contain only T or C (i.e., the pyrimidine context for base substitutions).

Usage

denominator_dict

Format

A vector with corresponding values

Examples

denominator_dict["type"]

Description

Computes the robust, ridge-penalized log-likelihood for a Dirichlet-Multinomial distribution. It safely handles massive counts utilizing lgamma components and matrix operations. The ridge penalty shrinks probability vectors toward a uniform distribution to stabilize optimization for sparse or heavily overdispersed data.

Usage

dm_loglik_ridge_stable(x, p, theta, lambda = 1)

Arguments

Value

A numeric scalar representing the computed log-likelihood value, adjusted by the penalty term.

Description

Runs through the plotting functions depending on data type and plot type.

Usage

f.plot.gui(
  ans.all,
  HTML = FALSE,
  model.summ = TRUE,
  display_plots = TRUE,
  .proast_env = NULL,
  output_type = NULL,
  filename = NULL,
  interactive_mode = TRUE,
  knitting = FALSE,
  return_plots = FALSE
)

Arguments

Value

Either the proast object (default) or a list of recorded plots if return_plots = TRUE.

Description

Independently generate the model plots from the raw results.

Usage

f.plot.result(
  proast_results_list,
  output_path = NULL,
  output_type = "svg",
  prefix = NULL,
  model_averaging = FALSE
)

Arguments

Value

Generates plots. Either saves them to an output path or records them and returns them as a list.

Description

Run dose-response modeling using PROAST.

Usage

f.proast(
  odt = list(),
  ans.all = 0,
  er = FALSE,
  resize = FALSE,
  scale.ans = FALSE,
  const.var = FALSE,
  show.warnings = FALSE,
  interactive_mode = TRUE,
  datatype = NULL,
  model_choice = NULL,
  setting_choice = NULL,
  nested_model_choice = NULL,
  indep_var_choice = NULL,
  Vyans_input = NULL,
  covariates = NULL,
  custom_CES = 0.05,
  model_selection = NULL,
  lower_dd = NULL,
  upper_dd = NULL,
  selected_model = NULL,
  adjust_CES_to_group_SD = NULL,
  model_averaging = NULL,
  num_bootstraps = NULL,
  display_plots = TRUE,
  add_nonzero_val_to_dat = FALSE,
  nonzero_val = NULL,
  detection_limit = NULL,
  seed = 125
)

Arguments

Value

Results from PROAST.

Description

This function creates a ⁠filter_mut`` column that will be read by the \code{calculate_mf} function and other downstream functions. Variants with ⁠filter_mut == TRUE“ will be excluded from group mutation counts. This function may also remove records upon on user specification. Running this function again on the same data will not overide the previous filters. To reset previous filters, set the filter_mut column values to FALSE.

Usage

filter_mut(
  mutation_data,
  vaf_cutoff = 1,
  snv_in_germ_mnv = FALSE,
  rm_abnormal_vaf = FALSE,
  custom_filter_col = NULL,
  custom_filter_val = NULL,
  custom_filter_rm = FALSE,
  regions = NULL,
  regions_filter,
  allow_half_overlap = FALSE,
  rg_sep = "\t",
  is_0_based_rg = TRUE,
  rm_filtered_mut_from_depth = FALSE,
  return_filtered_rows = FALSE
)

Arguments

Value

A data frame or a list of two data frames, depending on the value of return_filtered_rows. If return_filtered_rows is FALSE (default), a data frame of the same structure as mutation_data is returned, with an additional column, filter_mut, indicating whether each record has been flagged for filtering (TRUE) or not (FALSE). If return_filtered_rows is TRUE, a list containing two data frames is returned. The first data frame, named mutation_data, is the filtered mutation data as described above. The second data frame, named filtered_rows, contains all records that were either removed from mutation_data or flagged with filter_mut == TRUE.

Examples

# Mutation data is just for example purposes. It does not reflect real data.
mutation_data <- readRDS(system.file("extdata", "Example_files",
                                     "simple_mutation_data.rds",
                                     package = "MutSeqR"))
  # In this example, we will apply the following filters:
  # 1) Filter out putative germline variants using a VAF cutoff of 0.01
  # 2) Flag snv variants that overlap with germline mnv variants and
  # 3) Subtract the alt_depth of these variants from their total_depth
  #    (treat them as No-calls).
  filter_example <- filter_mut(
    mutation_data = mutation_data,
    vaf_cutoff = 0.01,
    snv_in_germ_mnv = TRUE,
    rm_filtered_mut_from_depth = TRUE,
    return_filtered_rows = FALSE
  )
 # Flagging germline mutations...
 # Found 15 germline mutations.
 # Flagging SNVs overlapping with germline MNVs...
 # Found 1 SNVs overlapping with germline MNVs.
 # Removing filtered mutations from the total_depth...
 # Filtering complete.

Description

This function will browse available BSgenomes, indicating which one should be installed for the specified organism and genome assembly version. If you cannot specify both organism and genome, the function can return a list of available genomes for a specified species.

Usage

find_BS_genome(organism, genome, masked = FALSE)

Arguments

Value

a BSgenome package name or a dataframe of possibilities

Examples

# Find the reference genome for Mouse, mm10 assembly:
mouse_mm10 <- find_BS_genome("mouse", "mm10")

Description

Uses the binomial distribution to create confidence intervals for mutation frequencies calculated from a single point estimate. Calculating binomial confidence intervals for mutation frequencies is not part of MutSeqR's recommended workflow, but is provided here for users who wish to use it.

Usage

get_binom_ci(
  mf_data,
  sum_col = "sum_min",
  depth_col = "group_depth",
  conf_level = 0.95,
  method = "wilson"
)

Arguments

Value

A mf data frame with added columns indicating the confidence intervals.

Examples

# Example data  consists of 24 mouse bone marrow
# samples exposed to three doses of BaP alongside vehicle controls.
# Libraries were sequenced with Duplex Sequencing using
# the TwinStrand Mouse Mutagenesis Panel which consists of 20 2.4kb
# targets = 48kb of sequence. Example data can be retrieved from
# MutSeqRData, an ExperimentHub data package:
## library(ExperimentHub)
## eh <- ExperimentHub()
## query(eh, "MutSeqRData")
# Mutation frequency data was precalculated using
## mf_data_global <- calculate_mf(mutation_data = eh[["EH9861"]],
##   cols_to_group = "sample",
##   retain_metadata_cols = c("dose_group", "dose"))

 mf <- readRDS(system.file("extdata", "Example_files",
                           "mf_data_global.rds",
                           package = "MutSeqR"))
 confint <- get_binom_ci(
   mf_data = mf,
    sum_col = "sum_min",
    depth_col = "group_depth"
  )

Description

Needs to be reworked for variants >1bp. Subset the mutation data and return only mutations that are found at positions with a specific motif. The default is CpG sites, but can be customizable.

Usage

get_cpg_mutations(
  mutation_data,
  regions,
  variant_types = c("-no_variant"),
  motif = "CG",
  filter_mut = TRUE
)

Arguments

Value

A GRanges object where each range is a mutation at a CpG site (a subset of mutations from the larger object provided to the function).

Description

Filters the ranges of your genomic regions to find all positions with a specific motif. The default is CpG sites, but can be customizable.

Usage

get_cpg_regions(regions, motif = "CG")

Arguments

Value

A GRanges object where each range is a CpG site (a subset of ranges from the larger object provided to the function).

Description

A utility function that will return the reference context of a mutation

Usage

get_ref_of_mut(mut_string)

Arguments

Value

the reference context of the mutation

Description

Create a GRanges object from the genomic target ranges and import raw nucleotide sequences.

Usage

get_seq(
  regions,
  rg_sep = "\t",
  is_0_based_rg = TRUE,
  padding = 0,
  BS_genome = NULL,
  ucsc = FALSE,
  species = NULL,
  genome = NULL
)

Arguments

Details

Consult available.genomes(splitNameParts=FALSE, type=getOption("pkgType")) for a full list of the available BS genomes and their associated species/genome/masked values. The BSgenome package will be installed if not already available. If using the UCSC API, the function will retrieve the sequences from the UCSC genome browser using the DAS API. See the UCSC website for available genomes: https://genome.ucsc.edu.

Value

a GRanges object with sequences of targeted regions.

Examples

#  Retrieve the sequences for custom regions
# We will load the TSpanel_human regions file as an example
# and supply it to the function as a GRanges object.
human <- load_regions_file("TSpanel_human")
regions_seq <- get_seq(
  regions = human,
  is_0_based_rg = FALSE,
  BS_genome = "BSgenome.Hsapiens.UCSC.hg38",
  padding = 0
)

Description

VariantAnnotation::rowRanges() defines the range width from REF, which is appropriate for small variants but does not honor structural-variant span annotations such as INFO/END. This helper prefers explicit END values from the INFO field and falls back to SVLEN for symbolic non-insertion alleles when END is absent.

Usage

get_vcf_end_positions(vcf)

Arguments

Value

An integer vector of end coordinates aligned to the rows of vcf.

Description

Imports tabular mutation file into the local R environment.

Usage

import_mut_data(
  mut_file,
  mut_sep = "\t",
  is_0_based_mut = TRUE,
  sample_data = NULL,
  sd_sep = "\t",
  regions = NULL,
  rg_sep = "\t",
  is_0_based_rg = TRUE,
  padding = 0,
  BS_genome = NULL,
  custom_column_names = NULL,
  output_granges = FALSE,
  add_chr = FALSE
)

Arguments

Details

Required columns for mut files are:

• contig: The name of the reference sequence.

• start: The start position of the feature.

• end: The half-open end position of the feature.

• sample: The sample name.

• ref: The reference allele at this position

• alt: The left-aligned, normalized, alternate allele at this position. Multiple alt alleles called for a single position should be represented as separate rows in the table.

The following columns are not required, but are recommended for full package functionality:

• alt_depth: The read depth supporting the alternate allele. If not included, the function will add this column, assuming a value of 1.

• total_depth: The total read depth at this position, excluding no-calls (N calls). If not present, the function will attempt to calculate the total_depth as depth - no_calls. If no_calls is not present, the function will use depth as the total_depth.

• depth: The total read depth at this position, including no-calls.

• no_calls: The number of no-calls (N-calls) at this position.

We recommend that files include a record for every sequenced position, regardless of whether a variant was called, along with the total_depth for each record. This enables site-specific depth calculations required for some downstream analyses.

Value

A table where each row is a mutation, and columns indicate the location, type, and other data. If output_granges is set to TRUE, the mutation data will be returned as a GRanges object, otherwise mutation data is returned as a dataframe.

Output Column Definitions:

• short_ref: The reference base at the start position.

• normalized_ref: The short_ref in C/T-base notation for this position (e.g. A -> T, G -> C).

• context The trinucleotide context at this position. Consists of the reference base and the two flanking bases (e.g. TAC).

• normalized_context: The trinucleotide context in C/T base notation for this position (e.g. TAG -> CTA).

• variation_type The type of variant (snv, mnv, insertion, deletion, complex, sv, no_variant, ambiguous, uncategorized).

• subtype The substitution type for the snv variant (12-base spectrum; e.g. A>C).

• normalized_subtype The C/T-based substitution type for the snv variant (6-base spectrum; e.g. A>C -> T>G).

• context_with_mutation: The substitution type for the snv variant including the two flanking nucleotides (192-trinucleotide spectrum; e.g. ⁠T[A>C]G⁠)

• normalized_context_with_mutation: The C/T-based substitution type for the snv variant including the two flanking nucleotides (96-base spectrum e.g. ⁠T[A>C]G⁠ -> ⁠C[T>G]A⁠).

• nchar_ref: The length (in bp) of the reference allele.

• nchar_alt: The length (in bp) of the alternate allele.

• varlen: The length (in bp) of the variant.

• ref_depth: The depth of the reference allele. Calculated as total_depth - alt_depth, if applicable.

• vaf : The variant allele fraction. Calculated as alt_depth/total_depth.

• gc_content: % GC of the trinucleotide context at this position.

• is_known: TRUE or FALSE. Flags known variants (ID != ".").

• row_has_duplicate: TRUE or FALSE. Flags rows whose position is the same as that of at least one other row for the same sample.

• filter_mut : A logical value, initially set to FALSE that indicates to calculte_mf() if the variant should be excluded from mutation counts. See the filter_mut function for more detail.

Examples

# Mutation data is just for example purposes. It does not reflect real data
file <- system.file("extdata", "Example_files",
                   "simple_mut_import.txt", package = "MutSeqR")
# Import the data
imported_example_data <- import_mut_data(mut_file = file)

Description

This function imports the regions metadata and joins it with the mutation data.

Usage

import_regions_metadata(
  mutation_granges,
  regions,
  rg_sep,
  is_0_based_rg,
  padding
)

Arguments

Value

A GRanges object that combines the mutation data with the regions metadata.

Description

This function imports sample metadata from a file or accepts a data frame directly, and performs basic validation checks.

Usage

import_sample_data(sample_data, sd_sep = "\t")

Arguments

Value

A validated data frame containing sample metadata, including a required column named sample.

Description

The function reads VCF file(s) and extracts the data into a dataframe.

Usage

import_vcf_data(
  vcf_file,
  sample_data = NULL,
  sd_sep = "\t",
  regions = NULL,
  rg_sep = "\t",
  is_0_based_rg = FALSE,
  padding = 0,
  BS_genome = NULL,
  output_granges = FALSE,
  remove_sample_suffix = NULL,
  add_chr = FALSE
)

Arguments

Details

The required fields are:

FIXED FIELDS

• CHROM: The name of the reference sequence. Equivalent to contig.

• POS: The 1-based start position of the feature. Equivalent to start.

• REF: The reference allele at this position.

• ALT: The left-aligned, normalized, alternate allele at this position. Multiple alt alleles called for a single position should be represented as separate rows in the table.

INFO FIELDS

• sample: An identifying field for your samples; either in the INFO field or as the header to the FORMAT field.

SUGGESTED FIELDS

The following FORMAT fields are not required, but are recommended for full package functionality:

• AD: The allelic depths for the reference and alternate allele in the order listed. The sum of AD is equivalent to the total_depth (read depth at this position excluding N-calls).

• DP: The read depth at this position (including N-calls). Equivalent to depth. Note that in many VCF files, the DP field is defined as total_depth. However, in most cases, the DP field includes N-calls.

• VD: The read depth supporting the alternate allele. If not included, the function will add this column, assuming a value of 1. Equivalent to alt_depth.

We recommend that files include a record for every sequenced position, regardless of whether a variant was called, along with the AD for each record. This enables site-specific depth calculations required for some downstream analyses. AD is used to calculate the total_depth (the read depth excluding No-calls). If AD is not available, the DP field will be used as the total_depth.

Value

A table where each row is a mutation, and columns indicate the location, type, and other data. If output_granges is set to TRUE, the mutation data will be returned as a GRanges object, otherwise mutation data is returned as a dataframe.

Output Column Definitions:

• 'end': The half-open end position of the feature. This is calculated as the 1-based start position + the length of the reference allele - 1. For structural variants, if the INFO field contains SVLEN, end is calculated as the 1-based start position + SVLEN - 1. End is only calculated if not already present in the VCF INFO fields.

• short_ref: The reference base at the start position.

• normalized_ref: The short_ref in C/T-base notation for this position (e.g. A -> T, G -> C).

• context The trinucleotide context at this position. Consists of the reference base and the two flanking bases (e.g. TAC).

• normalized_context: The trinucleotide context in C/T base notation for this position (e.g. TAG -> CTA).

• variation_type The type of variant (snv, mnv, insertion, deletion, complex, sv, no_variant, ambiguous, uncategorized).

• subtype The substitution type for the snv variant (12-base spectrum; e.g. A>C).

• normalized_subtype The C/T-based substitution type for the snv variant (6-base spectrum; e.g. A>C -> T>G).

• context_with_mutation: The substitution type for the snv variant including the two flanking nucleotides (192-trinucleotide spectrum; e.g. ⁠T[A>C]G⁠)

• normalized_context_with_mutation: The C/T-based substitution type for the snv variant including the two flanking nucleotides (96-base spectrum e.g. ⁠T[A>C]G⁠ -> ⁠C[T>G]A⁠).

• nchar_ref: The length (in bp) of the reference allele.

• nchar_alt: The length (in bp) of the alternate allele.

• varlen: The length (in bp) of the variant.

• ref_depth: The depth of the reference allele. Calculated as total_depth - alt_depth, if applicable.

• vaf : The variant allele fraction. Calculated as alt_depth/total_depth.

• gc_content: % GC of the trinucleotide context at this position.

• is_known: TRUE or FALSE. Flags known variants (ID != ".").

• row_has_duplicate: TRUE or FALSE. Flags rows whose position is the same as that of at least one other row for the same sample.

• filter_mut : A logical value, initially set to FALSE that indicates to calculte_mf() if the variant should be excluded from mutation counts. See the filter_mut function for more detail.

Examples

# Mutation data is just for example purposes. It does not reflect real data
file <- system.file("extdata", "Example_files",
                   "simple_vcf_data.vcf", package = "MutSeqR")
# Import the data
imported_example_data <- import_vcf_data(
 vcf_file = file,
BS_genome = find_BS_genome("mouse", "mm10"))

Description

A helper function to import the regions metadata file and return a GRanges object.

Usage

load_regions_file(regions, rg_sep = "\t", is_0_based_rg = TRUE)

Arguments

Value

a GRanges object of the imported regions metadata file.

Examples

#' # Example 1: Load built-in TwinStrand's Human Mutagenesis
human_rg <- load_regions_file(regions = "TSpanel_human")
human_rg
# Load a custom regions file from an interval list
# We will use the human TSpanel system file for this example,
# but any file can be imported.
file <- system.file("extdata",
  "inputs",
  "metadata",
  "human_mutagenesis_panel_hg38.txt",
  package = "MutSeqR"
)
custom_rg <- load_regions_file(regions = file, rg_sep = "\t", is_0_based_rg = TRUE)
custom_rg

Description

model_mf will fit a linear model to analyse the effect(s) of given factor(s) on mutation frequency and perform specified pairwise comparisons. This function will fit either a generalized linear model (glm) or, if supplied random effects, a generalized linear mixed-effects model (glmer). Pairwise comparisons are conducted using the doBy library (esticon) and estimates are then back-transformed. The delta method is employed to approximate the back-transformed standard-errors. A Sidak correction is applied to adjust p-values for multiple comparisons.

Usage

model_mf(
  mf_data,
  fixed_effects,
  test_interaction = TRUE,
  random_effects = NULL,
  reference_level,
  muts = "sum_min",
  total_count = "group_depth",
  contrasts = NULL,
  cont_sep = "\t",
  ...
)

Arguments

Details

fixed_effects are variables that have a direct and constant effect on the dependent variable (ie mutation frequency).They are typically the experimental factors or covariates of interest for their impact on the dependent variable. One or more fixed_effect may be provided. If you are providing more than one fixed effect, avoid using correlated variables; each fixed effect must independently predict the dependent variable. Ex. fixed_effects = c("dose", "genomic_target", "tissue", "age", etc).

Interaction terms enable you to examine whether the relationship between the dependent and independent variable changes based on the value of another independent variable. In other words, if an interaction is significant, then the relationship between the fixed effects is not constant across all levels of each variable. Ex. Consider investigating the effect of dose group and tissue on mutation frequency. An interaction between dose and tissue would capture whether the dose response differs between tissues.

random_effects account for the unmeasured sources of statistical variance that affect certain groups in the data. They help account for unobserved heterogeneity or correlation within groups. Ex. If your model uses repeated measures within a sample, random_effects = "sample".

Setting a reference_level for your fixed effects enhances the interpretability of the model. Ex. Consider a fixed_effect "dose" with levels 0, 25, 50, and 100 mg/kg. Intuitively, the reference_level would refer to the negative control dose, "0" since we are interested in testing how the treatment might change mutation frequency relative to the control.

Examples of contrasts:

If you have a fixed_effect "dose" with dose groups 0, 25, 50, 100, then the first column would contain the treated groups (25, 50, 100), while the second column would be 0, thus comparing each treated group to the control group.

25 0

50 0

100 0

Alternatively, if you would like to compare mutation frequency between treated dose groups, then the contrast table would look as follows, with the lower dose always in the second column, as we expect it to have a lower mutation frequency. Keeping this format aids in interpretability of the estimates for the pairwise comparisons. Should the columns be reversed, with the higher group in the second column, then the model will compute the fold-decrease instead of the fold-increase.

100 25

100 50

50 25

Ex. Consider the scenario where the fixed_effects are "dose" (0, 25, 50, 100) and "genomic_target" ("chr1", "chr2"). To compare the three treated dose groups to the control for each genomic target, the contrast table would look like:

25:chr1 0:chr1

50:chr1 0:chr1

100:chr1 0:chr1

25:chr2 0:chr2

50:chr2 0:chr2

100:chr2 0:chr2

Troubleshooting: If you are having issues with convergence for your generalized linear mixed-effects model, it may be advisable to increase the tolerance level for convergence checking during model fitting. This is done through the control argument for the lme4::glmer function. The default tolerance is tol = 0.002. Add this argument as an extra argument in the model_mf function. Ex. control = lme4::glmerControl(check.conv.grad = lme4::.makeCC("warning", tol = 3e-3, relTol = NULL)) Alternate approach: control = lme4::glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e5)) Similar approaches may be taken for glm models.

Value

Model results are output as a list. Included are:

• model_data: the supplied mf_data with added column for the Pearson's residuals of the model.

• summary: the summary of the model.

• anova: the analysis of variance for models with two or more effects. Anova(model)

• residuals_histogram: the Pearson's residuals plotted as a histogram. This is used to check whether the variance is normally distributed. A symmetric bell-shaped histogram, evenly distributed around zero indicates that the normality assumption is likely to be true.

• residuals_qq_plot: the Pearson's residuals plotted in a quantile-quantile plot. For a normal distribution, we expect points to roughly follow the y=x line.

• point_estimates_matrix: the contrast matrix used to generate point-estimates for the fixed effects.

• point_estimates: the point estimates for the fixed effects.

• pairwise_comparisons_matrix: the contrast matrix used to conduct the pairwise comparisons specified in the contrasts.

• pairwise_comparisons: the results of pairwise comparisons specified in the contrasts.

Examples

# Example data  consists of 24 mouse bone marrow
# samples exposed to three doses of BaP alongside vehicle controls.
# Libraries were sequenced with Duplex Sequencing using
# the TwinStrand Mouse Mutagenesis Panel which consists of 20 2.4kb
# targets = 48kb of sequence. Example data can be retrieved from
# MutSeqRData, an ExperimentHub data package:
## library(ExperimentHub)
## eh <- ExperimentHub()
## query(eh, "MutSeqRData")
# Mutation frequency data was precalculated using
## mf_data_global <- calculate_mf(mutation_data = eh[["EH9861"]],
##   cols_to_group = "sample",
##   retain_metadata_cols = c("dose_group", "dose"))

# We will model the effect of dose on mutation frequency min.
mf_example <- readRDS(system.file("extdata/Example_files/mf_data_global.rds",
  package = "MutSeqR"
))
# We will compare all BaP dose groups to the control group
# using pairwise comparisons.
contrasts <- data.frame(
  col1 = c("12.5", "25", "50"),
  col2 = c("0", "0", "0")
)
# Fit the model
model <- model_mf(
  mf_data = mf_example,
  fixed_effects = "dose",
  reference_level = "0",
  muts = "sum_min",
  total_count = "group_depth",
  contrasts = contrasts
)

Description

A list of column specifications

Usage

op

Format

A list with potential variable column names

Examples

# examples of "alt" synonyms
op$column$alternate

Description

Produces a ggplot object of bubble plots from given mutation data. Optionally, bubble plots can be facetted and coloured by a specified column.

Usage

plot_bubbles(
  mutation_data,
  size_by = "alt_depth",
  facet_col = NULL,
  color_by = "normalized_subtype",
  circle_spacing = 1,
  circle_outline = "none",
  circle_resolution = 50,
  custom_palette = NULL
)

Arguments

Details

The function will plot a circle for each mutation in mutation_data. Mutations flagged by the filter_mut column will be excluded from the plot. The size of the circle is determined by the size_by parameter. Sizing by the "alt_depth" or the "vaf" will give users the ability to visualize the the distribution of recurrent mutations within their data with large multiplets having a large circle.

Value

A ggplot object with the bubble plot, facetted if specified.

Examples

# The example data is a subset of variants from the target chr1
# from samples of the high dose group (50mg).
example_data <- readRDS(system.file("extdata", "Example_files",
                                    "variants_subset_d50_chr1.rds",
                                     package = "MutSeqR")
)
 plot <- plot_bubbles(
    mutation_data = example_data
  )

Description

Plot confidence intervals

Usage

plot_ci(
  data,
  order = "none",
  custom_order = NULL,
  nudge = 0.3,
  log_scale = FALSE,
  x_lab = NULL,
  y_lab = NULL,
  title = NULL
)

Arguments

Value

a ggplot object

Examples

if (requireNamespace("MutSeqRData", quietly = TRUE)) {
  # Plot results from PROAST
  dat <- data.frame(
    Response = c("PROAST MF Min", "PROAST MF Max"),
    BMD = c(NA, NA),
    BMDL = c(7.38, 2.98),
    BMDU = c(10.9, 7.68)
  )
  plot <- plot_ci(dat)
}

Description

Plot recurrent mutations in a lollipop plot using ggplot2

Usage

plot_lollipop(
  mutation_data,
  min_recurrence = 2,
  group_col = "contig",
  subtype_resolution = "base_6",
  custom_palette = NULL
)

Arguments

Value

A list of ggplot objects.

Examples

# For this example, we will use a subset of the example mutation data.
# The subset contains mutations from target chr1 in samples from the high
# dose group (50mg).
 example_data <- readRDS(system.file("extdata", "Example_files",
     "variants_subset_d50_chr1.rds", package = "MutSeqR"))
 # We will plot mutations that recoccur in at least two samples, grouped
 # by the "label" column, which signifies the target region (chr1).
 # Mutations will be grouped and coloured by their base 6 subtype (default) 
 plot <- plot_lollipop(
     mutation_data = example_data,
     min_recurrence = 2,
     group_col = "label"
 )

Description

This function calculates the mean mutation frequency across samples for given groups and plots the results.

Usage

plot_mean_mf(
  mf_data,
  group_col = "dose",
  fill_col = NULL,
  mf_type = "both",
  plot_type = "line",
  plot_error_bars = TRUE,
  plot_indiv_vals = TRUE,
  group_order = "none",
  group_order_input = NULL,
  add_labels = "mean_count",
  scale_y_axis = "linear",
  x_lab = NULL,
  y_lab = NULL,
  plot_title = NULL,
  custom_palette = NULL,
  plot_legend = TRUE,
  rotate_labels = FALSE,
  label_size = 3
)

Arguments

Value

a ggplot object

Examples

# Example data  consists of 24 mouse bone marrow
# samples exposed to three doses of BaP alongside vehicle controls.
# Libraries were sequenced with Duplex Sequencing using
# the TwinStrand Mouse Mutagenesis Panel which consists of 20 2.4kb
# targets = 48kb of sequence. Example data can be retrieved from
# MutSeqRData, an ExperimentHub data package:
## library(ExperimentHub)
## eh <- ExperimentHub()
## query(eh, "MutSeqRData")
# Mutation frequency data was precalculated using
## mf_data_global <- calculate_mf(mutation_data = eh[["EH9861"]],
##   cols_to_group = "sample",
##   retain_metadata_cols = c("dose_group", "dose"))

mf_example <- readRDS(system.file("extdata/Example_files/mf_data_global.rds",
  package = "MutSeqR"
))
# Specify the order of the dose groups along the x-axis
mf_example$dose_group <- factor(mf_example$dose_group,
  levels = c(
    "Control", "Low",
    "Medium", "High"
  )
)
# Plot the mean min MF per dose group as a bar plot with error bars
plot <- plot_mean_mf(
  mf_data = mf_example,
  group_col = "dose_group",
  mf_type = "min",
  plot_type = "line",
  fill_col = "dose_group",
  plot_error_bars = TRUE,
  plot_indiv_vals = TRUE,
  add_labels = "none"
)

Description

This function creates a plot of the mutation frequency.

Usage

plot_mf(
  mf_data,
  group_col,
  plot_type = "bar",
  mf_type = "min",
  fill_col = NULL,
  custom_palette = NULL,
  group_order = "none",
  group_order_input = NULL,
  labels = "count",
  scale_y_axis = "linear",
  x_lab = NULL,
  y_lab = NULL,
  title = NULL,
  rotate_labels = FALSE,
  label_size = 3
)

Arguments

Value

A ggplot object

Examples

# Example data  consists of 24 mouse bone marrow
# samples exposed to three doses of BaP alongside vehicle controls.
# Libraries were sequenced with Duplex Sequencing using
# the TwinStrand Mouse Mutagenesis Panel which consists of 20 2.4kb
# targets = 48kb of sequence. Example data can be retrieved from
# MutSeqRData, an ExperimentHub data package:
## library(ExperimentHub)
## eh <- ExperimentHub()
## query(eh, "MutSeqRData")
# Mutation frequency data was precalculated using
## mf_data_global <- calculate_mf(mutation_data = eh[["EH9861"]],
##   cols_to_group = "sample",
##   retain_metadata_cols = c("dose_group", "dose"))
# Load Example MF data
mf_example <- readRDS(system.file("extdata/Example_files/mf_data_global.rds",
  package = "MutSeqR"
))
# Specify the order of the dose groups along the x-axis
mf_example$dose_group <- factor(mf_example$dose_group,
  levels = c(
    "Control", "Low",
    "Medium", "High"
  )
)
# Plot the min MF per sample as a bar plot with count labels
plot <- plot_mf(
  mf_data = mf_example,
  group_col = "sample",
  plot_type = "bar",
  mf_type = "min",
  fill_col = "dose_group",
  group_order = "arranged",
  group_order_input = "dose_group",
  labels = "count",
  title = "Mutation Frequency per Sample"
)

Description

Provide a visualization of the point estimates derived using model_mf()

Usage

plot_model_mf(
  model,
  plot_type = "point",
  x_effect = NULL,
  plot_error_bars = TRUE,
  plot_signif = TRUE,
  ref_effect = NULL,
  x_order = NULL,
  fill_order = NULL,
  x_label = NULL,
  y_label = NULL,
  plot_title = NULL,
  fill_label = NULL,
  custom_palette = NULL
)

Arguments

Details

See model_mf() for examples.

Value

A ggplot object.

Examples

# Example data consists of 24 mouse bone marrow DNA samples imported
# using import_mut_data() and filtered with filter_mut.
# Data was summarized per sample using calculate_mf() and modeled using
# model_mf() (see example).
file <- system.file("extdata/Example_files/mf_model_global.rds",
     package = "MutSeqR")
model <- readRDS(file)
# Plot the results using plot_model_mf()
plot <- plot_model_mf(model,
  plot_type = "bar",
  x_effect = "dose",
  plot_error_bars = TRUE,
  plot_signif = TRUE,
  x_order = c("0", "12.5", "25", "50"),
  x_label = "Dose (mg/kg-bw/d)",
  y_label = "Estimated Mean MF (mutations/bp)",
  plot_title = ""
)

Description

Create a radar plot

Usage

plot_radar(mf_data, response_col, label_col, facet_col, indiv_y = TRUE)

Arguments

Value

A radar plot

Examples

# Plot the MFmin for each variation type, including the 6 SNV subtypes
# Facet the plots by dose group
mf_ex <- readRDS(system.file("extdata", "Example_files", "mf_data_6.rds",
     package = "MutSeqR")
)
plot <- plot_radar(
     mf_data = mf_ex,
     response_col = "mf_min",
     label_col = "normalized_subtype",
     facet_col = "dose_group",
     indiv_y = TRUE
)

Description

Given mf data, construct a plot displaying the mutation subtypes observed in a cohort.

Usage

plot_spectra(
  mf_data,
  group_col = "sample",
  subtype_resolution = "base_6",
  response = "proportion",
  mf_type = "min",
  group_order = "none",
  group_order_input = NULL,
  dist = "cosine",
  cluster_method = "ward.D",
  custom_palette = NULL,
  x_lab = NULL,
  y_lab = NULL,
  rotate_xlabs = FALSE
)

Arguments

Value

A ggplot object representing the mutation spectra plot.

Examples

# Example data consists of 24 mouse bone marrow DNA samples imported
# using import_mut_data() and filtered with filter_mut. Filtered
# mutation data is available in the MutSeqRData ExperimentHub package:
# eh <- ExperimentHub::ExperimentHub()
# Example 1: Visualized the 6-base mutation proportions per dose group.
# Data was summarized per dose_group using:
# calculate_mf(mutation_data = eh[["EH9861"]],
#              cols_to_group = "dose_group",
#              subtype_resolution = "base_6")
# Load the example data
mf_example <- readRDS(system.file("extdata", "Example_files", "mf_data_6.rds",
  package = "MutSeqR"
))
# Convert dose_group to a factor with the desired order.
mf_example$dose_group <- factor(mf_example$dose_group,
  levels = c("Control", "Low", "Medium", "High")
)
# Plot the mutation spectra
plot <- plot_spectra(
  mf_data = mf_example,
  group_col = "dose_group",
  subtype_resolution = "base_6",
  response = "proportion",
  group_order = "arranged",
  group_order_input = "dose_group"
)

# Example 2: plot the proportion of 6-based mutation subtypes
# for each sample, ordered by hierarchical clustering:
# Data was summarized per dose_group using:
# calculate_mf(mutation_data = eh[["EH9861"]],
#              cols_to_group = "sample",
#              subtype_resolution = "base_6")
# Load the example data
mf_example2 <- readRDS(system.file("extdata", "Example_files", "mf_data_6_sample.rds",
  package = "MutSeqR"
))
plot <- plot_spectra(
  mf_data = mf_example2,
  group_col = "sample",
  subtype_resolution = "base_6",
  response = "proportion",
  group_order = "clustered"
)

Description

Creates barplots of the trinucleotide spectrum for all levels of a given group.

Usage

plot_trinucleotide(
  mf_96,
  response = "proportion",
  mf_type = "min",
  group_col = "dose",
  indiv_y = FALSE,
  sum_totals = TRUE,
  output_path = NULL,
  output_type = "svg"
)

Arguments

Details

The function plots the trinucleotide spectrum for all levels of a given group from the provided mf_96 data; the output of calculate_mf with subtype_resolution = "base_96".

Value

A named list containing ggplots.

Examples

# Calculate the mutation frequency data at the 96-base resolution
 mf_96 <- readRDS(system.file("extdata", "Example_files", "mf_data_96.rds",
    package = "MutSeqR"))
# Plot the trinucleotide proportions for the control and high dose groups
mf_96 <- dplyr::filter(mf_96, dose_group %in% c("Control", "High"))

# Scale y-axis the same for all groups
  plots <- plot_trinucleotide(
    mf_96 = mf_96,
    response = "proportion",
    mf_type = "min",
    group_col = "dose_group",
    indiv_y = FALSE,
    output_path = NULL
  )

Description

This function creates a heatmap plot of subtype proportions for a given grouping variable. The groups may be facetted by a second variable. Mutation sums for each facet group and normalized subtype are calculated and displayed.

Usage

plot_trinucleotide_heatmap(
  mf_data,
  group_col = "sample",
  facet_col = "dose",
  mf_type = "min",
  mut_proportion_scale = "turbo",
  max = 0.2,
  rescale_data = FALSE,
  condensed = FALSE
)

Arguments

Value

A ggplot object representing the heatmap plot.

Examples

mf_96 <- readRDS(system.file("extdata/Example_files/mf_data_96_sample.rds",
package = "MutSeqR"))
# define dose_group order
mf_96$dose_group <- factor(mf_96$dose_group,
     levels = c("Control", "Low","Medium", "High")
)
plot <- plot_trinucleotide_heatmap(mf_96,
  group_col = "sample",
  facet_col = "dose_group"
)

Description

This function populates the trinucleotide context for each mutation in the mutation data.

Usage

populate_sequence_context(mutation_granges, BS_genome, n = 1)

Arguments

Value

A GRanges object with an additional column for the trinucleotide context.

Description

This function prints ASCII art when the package is loaded

Usage

print_ascii_art()

Value

None

Description

Rapidly simulates count matrices drawn from a Dirichlet-Multinomial distribution. This function is highly optimized for parametric bootstrapping, utilizing a Gamma-Multinomial approximation for speed.

Usage

rdm_fast(depths, p, theta)

Arguments

Value

A numeric matrix of simulated mutation counts, with subtypes as rows and simulated samples as columns.

Description

Map column names of mutation data to default column names. A utility function that renames columns of mutation data to default columns names.

Usage

rename_columns(data, column_map = op$column)

Arguments

Value

the mutation data with column names changed to match default.

Examples

df <- data.frame(
  chromosome = c("chr1", "chr2", "chr3"),
  pos = c(100, 200, 300),
  end = c(100, 200, 300),
  sample_id = c("S1", "S2", "S3"),
  reference = c("G", "C", "T"),
  alternate = c("A", "T", "G")
)
renamed_data <- rename_columns(df, column_map = op$column)

Description

This function reads a configuration file in YAML format, extracts the parameters, and renders an R Markdown document using the specified parameters.

Usage

render_report(
  config_filepath,
  output_file = "./MutSeqR_Summary_Report.html",
  output_format = "html_document"
)

Arguments

Value

A rendered R Markdown document.

Examples

# Step 1: Copy the example configuration file to your working directory
## config <- system.file("extdata", "inputs", "summary_config.yaml", package = "MutSeqR")
## file.copy(from = config, to = "your/working/directory/summary_config.yaml")
# Step 2: Edit the configuration file with your inputs
# Step 3: Render the report
## render_report(config_filepath = "your/working/directory/summary_config.yaml",
##  output_file = "MutSeqR_Summary_Report.html",
##  output_format = "html_document")

Description

Get the reverse complement of a DNA or RNA sequence.

Usage

reverseComplement(
  x,
  content = c("dna", "rna"),
  case = c("lower", "upper", "as is")
)

Arguments

Details

This file is part of the source code for SPGS: an R package for identifying statistical patterns in genomic sequences. Copyright (C) 2015 Universidad de Chile and INRIA-Chile A copy of Version 2 of the GNU Public License is available in the share/licenses/gpl-2 file in the R installation directory or from http://www.R-project.org/Licenses/GPL-2. reverseComplement.R

Value

A character vector of the reverse complement sequences.

Description

Core statistical engine for the RP-DMM test. Evaluates whether two count matrices represent significantly different mutational spectra. It utilizes a "hybrid logic": robust parameters (thetas) are estimated utilizing a ridge penalty (QA phase), while the test statistic (LRT) is evaluated using the unpenalized likelihood. Significance is evaluated via parametric bootstrap.

Usage

run_penalized_comparison(g1, g2, lambda = 1, n_boot = 500)

Arguments

Value

A list containing 5 items:

• p_value: Numeric. The bootstrapped p-value.

• lrt: Numeric. The observed unpenalized Likelihood Ratio Test statistic.

• theta_g1: Numeric. The optimized robust theta (overdispersion) for Group 1.

• theta_g2: Numeric. The optimized robust theta for Group 2.

• theta_shared: Numeric. The optimized robust theta under the Null hypothesis.

Description

Performs leave-one-out cross-validation (LOOCV) to empirically select the optimal ridge penalty (lambda) for a given count matrix. Models are fit on $N-1$ samples using a specified lambda, and unpenalized likelihood is evaluated on the held-out sample.

Usage

select_optimal_lambda(count_matrix, lambdas = c(0.1, 0.5, 1, 5, 10, 50, 100))

Arguments

Value

A list containing two items:

• best_lambda: Numeric scalar. The lambda value from the provided vector that yielded the highest overall cross-validation log-likelihood.

• scores: A named numeric vector detailing the cumulative log-likelihood score achieved by each tested lambda.

Description

This function initializes the Python environment used by MutSeqR. It is not run automatically to avoid issues during installation and checks.

Usage

setup_mutseqr_python(force = FALSE)

Arguments

Value

None

Description

Correct p-values for multiple comparisons

Usage

sidak(vecP)

Arguments

Details

This function corrects a vector of probabilities for multiple testing using the Bonferroni (1935) and Sidak (1967) corrections. References: Bonferroni (1935), Sidak (1967), Wright (1992). Bonferroni, C. E. 1935. Il calcolo delle assicurazioni su gruppi di teste. Pp. 13-60 in: Studi in onore del Professore Salvatore Ortu Carboni. Roma. Sidak, Z. 1967. Rectangular confidence regions for the means of multivariate normal distributions. Journal of the American Statistical Association 62:626-633. Wright, S. P. 1992. Adjusted P-values for simultaneous inference. Biometrics 48: 1005-1013. Pierre Legendre, May 2007

Value

adjusted p-values

Examples

p_values <- c(0.01, 0.04, 0.03, 0.08, 0.05)
adjusted_p <- sidak(p_values)
adjusted_p$SidakP

Description

Run COSMIC signatures comparison using SigProfilerAssignment

Usage

signature_fitting(
  mutation_data,
  project_name = "Default",
  project_genome = "GRCh38",
  env_name = "MutSeqR",
  group = "sample",
  output_path = NULL,
  python_version
)

Arguments

Details

Assign COSMIC SBS signatures to mutation data using SigProfilerAssignment. Data is cleaned and formatted for input into SigProfiler tools. This function will create a virtual environment using reticulate to run python, as this is a requirement for the SigProfiler suite of tools. Note that it will also install several python dependencies using a conda virtual environment on first use. Please be aware of the implications of this. For advanced use, it is suggested to use the SigProfiler python tools directly in python as described in their respective documentation. Users must have python installed on their computer to use this function.

Mutation data will be filtered to only include SNVs. Variants flagged by the filter_mut column will be excluded.

Value

Creates a subfolder "SigProfiler" in the output directory with SigProfiler tools results. For a complete breakdown of the results, see the Readme file for MutSeqR. Most relevant results are stored in SigProfiler > ggplot2::group > matrices > output > Assignment_Solution > Activities > SampleReconstruction > WebPNGs. These plots show a summary of the signature assignment results for each group. In each plot, the top left panel represents the base_96 mutation count for the group. The bottom left panel represents the reconstructed profile. Below the reconstruction are the solution statistics that indicate the goodness of fit of the reconstructed profile to the observed profile. (Recommended cosine similarity > 0.9). The panels on the right represent the SBS signatures that contribute to the reconstructed profile. The signature name and its contribution % are shown in the panel. A high contribution means a high association of the signature with the group's mutation spectra.

Examples

if (requireNamespace("MutSeqRData", quietly = TRUE)) {
  # Example data consists of 24 mouse bone marrow DNA samples imported
  # using import_mut_data() and filtered with filter_mut as in Example 4.
  # Sequenced on TS Mouse Mutagenesis Panel. Example data is
  # retrieved from MutSeqRData, an ExperimentHub data package.
  library(ExperimentHub)
  eh <- ExperimentHub()
  example_data <- eh[["EH9861"]]
  output_path <- tempdir()

  signature_fitting(
    mutation_data = example_data,
    project_name = "Example",
    project_genome = "mm10",
    env_name = "MutSeqR",
    group = "dose",
    python_version = "3.11",
    output_path = output_path
  )
}

Description

spectra_comparison compares the mutation spectra of groups using a modified contingency table approach.

Usage

spectra_comparison(
  mf_data,
  exp_variable,
  mf_type = "min",
  contrasts,
  cont_sep = "\t"
)

Arguments

Details

This function creates an R * 2 contigency table of the subtype counts, where R is the number of subtypes for the 2 groups being compared. The G2 likelihood ratio statistic is used to evaluate whether the proportion (count/group total) of each mutation subtype equals that of the other group.

The G2 statistic refers to a chi-squared distribution to compute the p-value for large sample sizes. When N / (R-1) < 20, where N is the total mutation counts across both groups, the function will use a F-distribution to compute the p-value in order to reduce false positive rates.

The comparison assumes independance among the observations, as such, it is highly recommended to use mf_type = "min".

Examples of contrasts: For 'exp_variable = "dose"' with dose groups 0, 12.5, 25, 50, compare each treated dose to the control:

12.5 0

25 0

50 0

Ex. Consider two 'exp_variables = c("dose", "tissue")'; with levels dose (0, 12.5, 25, 50) and tissue("bone_marrow", "liver"). To compare the mutation spectra between tissues for each dose group, the contrast table would look like:

0:bone_marrow 0:liver

12.5:bone_marrow 12.5:liver

25:bone_marrow 25:liver

50:bone_marrow 50:liver

Value

the log-likelihood statistic G2 for the specified comparisons with the p-value adjusted for multiple-comparisons.

Examples

# Example data consists of 24 mouse bone marrow DNA samples imported
# using import_mut_data() and filtered with filter_mut. Filtered
# mutation data is available in the MutSeqRData ExperimentHub package:
# eh <- ExperimentHub::ExperimentHub()
# Data was summarized per sample using:
# calculate_mf(mutation_data = eh[["EH9861"]],
#              cols_to_group = "dose_group",
#              subtype_resolution = "base_6")

# Example: compare 6-base mutation spectra between dose groups
# Load the example data
mf_example <- readRDS(
     system.file("extdata", "Example_files", "mf_data_6.rds",
         package = "MutSeqR"
     )
)
# Create the contrasts table
contrasts <- data.frame(
  col1 = c("Low", "Medium", "High"),
  col2 = rep("Control", 3)
)
# Run the comparison
spectra_comparison(
  mf_data = mf_example,
  exp_variable = "dose_group",
  mf_type = "min",
  contrasts = contrasts
)

Description

Ridge-penalized Dirichlet-Multinomial test for comparing mutation spectra between two groups. Robust for large mutation counts.

Usage

spectra_comparison_rpdmm(
  mf_data,
  exp_variable,
  contrasts,
  cont_sep = "\t",
  mf_type = "min",
  lambda = 1,
  n_boot = 500
)

Arguments

Value

A data frame containing one row per specified contrast. Columns include the group names, comparison string, bootstrap p-value, observed Likelihood Ratio Test (LRT) statistic, and robust overdispersion metrics (theta) for group 1, group 2, and the shared null model.

Description

Associates the subtype resolution names with their corresponding column names in the mutation data.

Usage

subtype_dict

Format

A vector with corresponding values

Examples

subtype_dict["base_96"]
subtype_dict["type"]

Description

A comprehensive list of mutation subtypes at different resolutions

Usage

subtype_list

Format

A list with corresponding values

Examples

subtype_list[["type"]]

Description

Checks to find the sample name of the vcf in the INFO field or in the FORMAT header. Can also handle sample name synonyms.

Usage

vcf_sample_fix(vcf)

Arguments

Value

The vcf with sample column name corrected

Description

Writes data to an Excel file.

Usage

write_excel(data, output_path = NULL, workbook_name, model_results = FALSE)

Arguments

Value

A saved Excel workbook.

Examples

# Example data consists of 24 mouse bone marrow DNA samples imported
# using import_mut_data(), filtered with filter_mut, and summarized
# using calculate_mf().
outputpath <- tempdir()
mf_example <- readRDS(system.file("extdata/Example_files/mf_data_global.rds",
  package = "MutSeqR"
))
mf_example2 <- readRDS(system.file("extdata/Example_files/mf_data_6.rds",
  package = "MutSeqR"
))
mf_example3 <- readRDS(
 system.file("extdata/Example_files/mf_data_6_sample.rds",
     package = "MutSeqR"
 )
)
list <- list(mf_example, mf_example2, mf_example3)
names(list) <- c("Global MF", "Base 6 Spectra", "Base 6 Sample Spectra")

# save a single data frame to an Excel file
write_excel(
 mf_example,
 output_path = outputpath,
 workbook_name = "test_single"
)

# save a list of data frames to an Excel file
write_excel(list, output_path = outputpath, workbook_name = "test_list")

Description

Creates a .txt file from mutation data that can be used for mutational signatures analysis using the SigProfiler Assignment web application. Currently only supports SBS analysis i.e. snvs.

Usage

write_mutation_calling_file(
  mutation_data,
  project_name = "Example",
  project_genome = "GRCm38",
  output_path = NULL
)

Arguments

Details

Mutations will be be filtered for SNVs. Mutations flagged in filter_mut will be excluded from the output.

Value

a .txt file that can be uploaded to the SigProfiler Assignment web application (https://cancer.sanger.ac.uk/signatures/assignment/) as a "Mutational calling file".

Examples

# Example data  consists of 24 mouse bone marrow
# samples exposed to three doses of BaP alongside vehicle controls.
# Libraries were sequenced with Duplex Sequencing using
# the TwinStrand Mouse Mutagenesis Panel which consists of 20 2.4kb
# targets = 48kb of sequence. Example data can be retrieved from
# MutSeqRData, an ExperimentHub data package:
## library(ExperimentHub)
## eh <- ExperimentHub()
## query(eh, "MutSeqRData")
# The data is a subset of variants from the target chr1
# from samples of the high dose group (50mg).
example_data <- readRDS(system.file("extdata", "Example_files",
                                    "variants_subset_d50_chr1.rds",
                                     package = "MutSeqR")
)
 write_mutation_calling_file(
    mutation_data = example_data,
    project_name = "Example",
    project_genome = "GRCm38",
    output_path = tempdir()
  )
  list.files(tempdir())
  # The file is saved in the temporary directory
  # To view the file, use the following code:
  ## output_file <- file.path(tempdir(), "mutation_calling_file.txt")
  ## file.show(output_file)

Description

Creates a .txt file from mutation data that can be used for mutational signatures analysis using the SigProfiler web application. Can handle group analyses (ex dose, tissue, etc). Currently only supports SBS matrices i.e. snvs.

Usage

write_mutational_matrix(
  mutation_data,
  group = "dose",
  subtype_resolution = "base_96",
  mf_type = "min",
  output_path = NULL
)

Arguments

Details

Mutations will be be filtered for SNVs. Mutations flagged in filter_mut will be excluded from the output. Mutations will be summed across the groups specified in the group argument.

Value

a .txt file that can be uploaded to the SigProfiler Assignment web application (https://cancer.sanger.ac.uk/signatures/assignment/) as a "Mutational Matrix".

Examples

# Example data  consists of 24 mouse bone marrow
# samples exposed to three doses of BaP alongside vehicle controls.
# Libraries were sequenced with Duplex Sequencing using
# the TwinStrand Mouse Mutagenesis Panel which consists of 20 2.4kb
# targets = 48kb of sequence. Example data can be retrieved from
# MutSeqRData, an ExperimentHub data package:
## library(ExperimentHub)
## eh <- ExperimentHub()
## query(eh, "MutSeqRData")
# The data is a subset of variants from the target chr1
# from samples of the high dose group (50mg).
example_data <- readRDS(system.file("extdata", "Example_files",
                                    "variants_subset_d50_chr1.rds",
                                     package = "MutSeqR")
)

  write_mutational_matrix(
    mutation_data = example_data,
    group = "sample",
    subtype_resolution = "base_96",
    mf_type = "min",
    output_path = tempdir()
  )
  list.files(tempdir())
  # The file is saved in the temporary directory
  # To view the file, use the following code:
  ## output_file <- file.path(tempdir(), "sample_base_96_mutational_matrix.txt")
  ## file.show(output_file)

Description

Write FASTA file of reference sequences.

Usage

write_reference_fasta(regions_gr, output_path = NULL)

Arguments

Details

Generate an arbitrary multi-sequence FASTA file from GRanges including the reference sequences.

Value

Writes a FASTA reference file "reference_output.fasta". If multiple ranges are included in the GRanges object, the sequences will be written to a single FASTA file. Sequences names will be the seqnames (contig) of the range.

Examples

# Write FASTA files for the 20 genomic target sequences
# of TwinStrand's Mouse Mutagenesis Panel.
output_path <- tempdir()
rg <- get_seq("TSpanel_mouse")
write_reference_fasta(rg, output_path = output_path)

Description

Export your mutation_data to a VCF file for downstream applications.

Usage

write_vcf_from_mut(mutation_data, output_path = NULL)

Arguments

Value

Writes a VCF file of mutations "mutation_output.vcf".

Examples

# Example data  consists of 24 mouse bone marrow
# samples exposed to three doses of BaP alongside vehicle controls.
# Libraries were sequenced with Duplex Sequencing using
# the TwinStrand Mouse Mutagenesis Panel which consists of 20 2.4kb
# targets = 48kb of sequence. Example data can be retrieved from
# MutSeqRData, an ExperimentHub data package:
## library(ExperimentHub)
## eh <- ExperimentHub()
## query(eh, "MutSeqRData")
# The data is a subset of variants from the target chr1
# from samples of the high dose group (50mg).
example_data <- readRDS(system.file("extdata", "Example_files",
                                    "variants_subset_d50_chr1.rds",
                                     package = "MutSeqR")
)
# Export mutation data of the four samples to a multi-sample VCF file.
write_vcf_from_mut(mutation_data = example_data, output_path = tempdir())