--- title: "MutSeqR: Error-Corrected Sequencing (ECS) Analysis For Mutagenicity Assessment" author: - name: Annette E. Dodge affiliation: Environmental Health Science and Research Bureau, Health Canada, Ottawa, ON, Canada. - name: Matthew J. Meier affilitation: Environmental Health Science and Research Bureau, Health Canada, Ottawa, ON, Canada. email: matthew.meier@hc-sc.gc.ca package: MutSeqR output: rmarkdown::html_document: toc: true toc_float: yes code_folding: show code_download: yes css: custom.css bibliography: references.bib link-citations: true vignette: > %\VignetteIndexEntry{MutSeqR: Error-Corrected Sequencing (ECS) Analysis For Mutagenicity Assessment} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE) library(DT) library(htmltools) ``` # Introduction In this module we will provide instructions on how to import, filter, and summarize mutation data generated by Error-correct next-generation sequencing (ECS). For information on how to use MutSeqR for statistical analysis of the data, see the relevant modules: - Modelling Mutation Frequencies (Generalized Linear Modelling and BMD) - Analyzing Mutation Spectra (Comparison of Spectra and Mutation Signature Analysis) - Target Analysis - Other Features (visualization of clonal mutations, retrieve target sequences) ## What is ECS? Error-corrected next-generation sequencing (ECS) uses various methods to combine multiple independent raw sequence reads derived from an original starting molecule, thereby subtracting out artifacts introduced during sequencing or library preparation. This results in a highly accurate representation of the original molecule. ECS is particularly useful for detecting rare somatic mutations (or mutations induced in germ cells), such as those that arise from mutagen exposure or other sources of DNA damage. ECS is a powerful tool for assessing the mutagenicity of chemicals, drugs, or other agents, and can be used to identify the mutational signatures of these agents. ECS can also be used to detect rare mutations in cancer or other diseases, and to track the clonal evolution of these diseases over time. For more background on how ECS works and its context in regulatory toxicology testing and genetic toxicology, see the following articles: - @menon-2023 - @marchetti-2023a - @marchetti-2023b - @kennedy-2014 This R package is meant to facilitate the import, cleaning, and analysis of ECS data, beginning with a table of variant calls or a variant call file (VCF). The package is designed to be flexible and enable users to perform common statistical analyses and visualisations. # Installation Install the package from Bioconductor: ```{r install-bioc, eval=FALSE} if (!require("BiocManager", quietly = TRUE)) { install.packages("BiocManager") } BiocManager::install("MutSeqR") ``` Load the package ```{r load-lib} library(MutSeqR) ``` # Data import The main goal of MutSeqR is to generate summary statistics, visualizations, exploratory analyses, and other post-processing tasks such as mutational signature analysis or generalized linear modeling. Mutation data should be supplied as a table of variants with their genomic positions. Mutation data can be imported as either VCF files or as tabular data using the functions `import_vcf_data` and `import_mut_data`, respectively. It is reccomended that files include a record for every sequenced position, regardless of whether a variant was called or not, along with the `total_depth` for each record. This enables site-specific depth calculations that are required for the calculation of mutation subtype frequencies and other site-specific frequencies. The data set can be pared down later to include only mutations of interest (SNVs, indels, SVs, or any combination). Table: (\#tab:required-columns) Required columns for mutation data import. | **Column** | **VCF Specification** | **Definition** | | -------------------- | ------------------------------- | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ | | contig | CHROM | The name of the reference sequence. | | start | POS | The start position of the feature. 0-based coordinates are accepted but will be changed to 1-based during import. | | end | INFO(END) | The half-open end position of the feature in contig. | | ref | REF | The reference allele at this position. | | alt | ALT | The left-aligned, normalized, alternate allele at this position. | | sample | INFO(sample) or Genotype Header | A unique identifier for the sample library. For VCF files, this field may be provided in either the INFO field, or as the header to the GENOTYPE field. | | _SUGGESTED FIELDS_ | | | | alt_depth | VD | The read depth supporting the alternate allele. If not included, the function will assume an alt_depth of 1 at variant sites. | | total_depth or depth | AD or DP | The total read depth at this position. This column can be “total_depth” which excludes N-calls, or “depth”, which includes N-calls, if “total_depth” is not available. For VCF files, the total_depth is calculated as the sum of AD. DP is equivalent to "depth". | VCF files should follow the VCF specification (version 4.5; @danecek-2011). VCF files may be bg/g-zipped. Each individual VCF file should contain the mutation data for a single sample library. Multiple alt alleles called for a single position should be represented as separate rows in the data. All extra columns, INFO fields, and FORMAT fields will be retained upon import. Upon import, records are categorized within the `variation_type` column based on their REF and ALT. Categories are listed below. Table: (\#tab:variation-types) Definitions of Variation types. | variation_type | Definition | | -------------- | ----------------------------------------------------------------- | | no_variant | No variation, the null-case. | | snv | Single nucleotide variant. | | mnv | Multiple nucleotide variant. | | insertion | Insertion. | | deletion | Deletion. | | complex | REF and ALT are of different lengths and nucleotide compositions. | | symbolic | Structural variant | | ambiguous | ALT contains IUPAC ambiguity codes. | | uncategorized | The record does not fall into any of the preceding categories. | Additional columns are created to further characterise variants. Table: (\#tab:mut-columns) Definitions of Mutation data columns. | Column Name | Definition | | -------------------------------- | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------ | | short_ref | The reference base at the start position. | | normalized_ref | The short_ref in C/T (pyrimidine) notation for this position. Ex. `A` -> `T`, `G` -> `C` | | context | The trinucleotide context at this position. Consists of the reference base and the two flanking bases. Sequences are retrieved from the appropriate BS genome. Ex. `TAC` | | normalized_context | The trinucleotide context in C/T (pyrimidine) notation for this position (Ex. `TAG` -> `CTA`) | | variation_type | The type of variant (no_variant, snv, mnv, insertion, deletion, complex, sv, ambiguous, uncategorized) | | subtype | The substitution type of the snv variant (12-base spectrum; Ex. `A>C`) | | normalized_subtype | The snv subtype in C/T (pyrimidine) notation (6-base spectrum; Ex. `A>C` -> `T>G`) | | context_with_mutation | The snv subtype including the two flanking nucleotides (192-base spectrum; Ex. `T[A>C]G`) | | normalized_context_with_mutation | The snv subtype in C/T (pyrimidine) notation including the two flanking nucleotides (96-base spectrum; Ex. `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 | A logical value indicating if the record is a known variant; i.e. ID field is not NULL. | | row_has_duplicate | A logical value that flags rows whose position is the same as that of at least one other row for the same sample. | ## General Usage Indicate the filepath to your mutation data file(s) using the `vcf_file`/ `mut_file` parameter. This can be either a single file or a directory containing multiple files. Should you provide a directory, then all files within will be bound into a single data frame. An important component of importing your data for proper use is to assign each mutation to a biological sample, and also make sure that some additional information about each sample is present (ex., a chemical treatment, a dose, etc.). This is done by providing `sample_data`. This parameter can take a data frame, or it can read in a file if provided it with a filepath. If using a filepath, specify the proper delimiter using the `sd_sep` parameter. Sample metadata will be joined with the mutation data using the "sample" column to capture that information and associate it with each mutation. Specify the appropriate BS genome with which to populate the context column. Prior to data import, ensure that the BS genome is installed. Supply the pkgname to the BS_genome parameter in the import function. If you are not sure which BS genome you should use, you can use the find_BS_genome() function. This function will browse [BSgenome::available.genomes](https://www.rdocumentation.org/packages/BSgenome/versions/1.40.1/topics/available.genomes) for given organisms and genome assembly versions. Once you've identified the appropriate BS genome, it can be installed using BiocManager(). Context information will be extracted from the BS genome. Note, BSgenome offers genomes with masked sequences. If you wish to use the masked version of the genome, you can find the right one with find_BS_genome(masked = TRUE). ## Examples {.tabset} We provide an example data set taken from @leblanc-2022. This data consists of 24 mouse bone marrow samples sequenced with Duplex Sequencing using Twinstrand's Mouse Mutagenesis Panel of twenty 2.4kb targeted genomic loci. Mice were exposed to three doses of benzo[a]pyrene (BaP) alongside vehicle controls, n = 6. Example data is retrieved from MutSeqRData, an ExperimentHub data package available via Bioconductor. Example data can be retrieved from the ExperimentHub index (eh) through specific accessors. ```{r load-MutSeqRData} library(ExperimentHub) # load the index eh <- ExperimentHub() # Query the index for MutSeqRData records and their associated accessors. query(eh, "MutSeqRData") ``` Identify the BS genome: Sequencing data was aligned to the mm10 mouse genome. ```{r find_BS_genome} mouse_BS_genome <- find_BS_genome(organism = "mouse", genome = "mm10") print(mouse_BS_genome) ``` Install the appropriate BS genome: ```{r install_BS_genome, eval=FALSE} BiocManager("BSgenome.Mmusculus.UCSC.mm10") ``` ### VCF _Example 1.1. Import the example .vcf.bgz file "Example import_vcf_data"._ _Provided is the genomic vcf.gz file for sample dna00996.1. It is comprised_ _of a record for all 48K positions sequenced for the Mouse Mutagenesis Panel_ _with the alt_depth and the tota_depth values for each record._ ```{r import-vcf} example_file <- eh[["EH9859"]] sample_metadata <- data.frame( sample = "dna00996.1", dose = "50", dose_group = "High" ) # Import the data imported_example_data <- import_vcf_data( vcf_file = example_file, sample_data = sample_metadata, BS_genome = mouse_BS_genome ) ``` **Data Table 1.** Mutation data imported from VCF for sample dna00996.1. Displays the first 6 rows. ```{r import-vcf-table, echo=FALSE} DT::datatable(head(imported_example_data), options = list(scrollX = TRUE)) ``` ### Tabular _Example 1.2. Import example tabular data "Example import_mut_data". This is_ _the equivalent file to the example vcf file. It is stored as an .rds file._ _We will load the dataframe and supply it to `import_mut_data`. The mut_file_ _parameter can accept file paths or data frames as input._ ```{r import-mut} example_data <- eh[["EH9857"]] sample_metadata <- data.frame( sample = "dna00996.1", dose = "50", dose_group = "High" ) # Import the data imported_example_data <- import_mut_data( mut_file = example_data, sample_data = sample_metadata, BS_genome = find_BS_genome("mouse", "mm10"), # Must be installed!! is_0_based_mut = TRUE # indicates that the genomic coordinates are 0-based. # Coordinates will be changed to 1-based upon import. ) ``` **Data Table 2.** Mutation data imported from a data.frame object for sample dna00996.1. Displays the first 6 rows. ```{r import-mut-table, echo=FALSE} DT::datatable(head(imported_example_data), options = list(scrollX = TRUE)) ``` ## Region Metadata Similar to sample metadata, you may supply a file containing the metadata of genomic regions to the `regions` parameter. In the file, the genomic regions are defined by the name of their reference sequence plus their start and end coordinates. Additional columns are added to describe features of the genomic regions (Ex. transcription status, GC content, etc.). Region metadata will be joined with mutation data by checking for overlap between the target region ranges and the position of the record. The `regions` parameter can be either a filepath, a data frame, or a [GRanges](https://bioconductor.org/packages/devel/bioc/vignettes/GenomicRanges/inst/doc/GenomicRangesIntroduction.html) object. Filepaths will be read using the `rg_sep`. Users can also choose from the built-in TwinStrand DuplexSeq™ Mutagenesis Panels by inputting "TSpanel_human", "TSpanel_mouse", or "TSpanel_rat". Required columns for the regions file are "contig", "start", and "end". In a GRanges object, the required columns are "seqnames", "start", and "end". Users must indicate whether the region coordinates are 0-based or 1-based with `is_0_based_rg`. Mutation data and region coordinates will be converted to 1-based. If you do not wish to specify regions, then set the `regions` parameter to NULL (default). _Example 1.3. Add the metadata for TwinStrand's Mouse Mutagenesis panel to our_ _example tabular file._ ```{r import-mut-rg} imported_example_data <- import_mut_data( mut_file = example_data, sample_data = sample_metadata, BS_genome = find_BS_genome("mouse", "mm10"), is_0_based_mut = TRUE, regions = "TSpanel_mouse" ) ``` **Data Table 3.** Mutation data imported from a data.frame object for sample dna00996.1 with added metadata for the Mouse Mutagenesis Panels sequencing targets. Displays the first 6 rows. Added columns are 'target_size', 'label', 'genic_context', 'region_GC_content', 'genome', and 'in_regions'. ```{r import-mut-table-rg, echo=FALSE} DT::datatable(head(imported_example_data), options = list(scrollX = TRUE)) ``` Users can load the example TwinStrand Mutagenesis Panels with load_regions(). This function will output a GRanges object. Custom panels may also be loaded with this function by providing their file path to the regions parameter. ```{r load-rg} region_example <- load_regions_file("TSpanel_mouse") region_example ``` ## Custom Column Names We recognize that column names may differ from those specified in Table \@ref(tab:required-columns). Therefore, we have implemented some default column name synonyms. If your column name matches one of our listed synonyms, it will automatically be changed to match our set values. For example, your `contig` column may be named `chr` or `chromosome`. After importing your data, this synonymous column name will be changed to `contig`. Column names are case-insensitive. A list of column name synonyms are listed below. Table: (\#tab:name-sym) Predefined column name synonyms. Synonyms will be automatically changed to the default Column value upon import. | Column | Synonyms | | ----------- | ---------------------------------------------------------- | | alt | alternate | | alt_depth | alt_read_depth, var_depth, variant_depth, vd | | context | sequence_context, trinucleotide_context, flanking_sequence | | contig | chr, chromosome, seqnames | | depth | dp | | end | end_position | | no_calls | n_calls, n_depth, no_depth | | ref | reference, ref_value | | sample | sample_name, sample_id | | start | pos, position | | total_depth | informative_somatic_depth | If your data contains a column that is synonymous to one of the required columns, but the name is not included in our synonyms list, your column name may be substituted using the `custom_column_names` parameter. Provide this parameter with a list of names to specify the meaning of column headers. ```{r import-mut-custom-cols, message=TRUE} mut_data <- eh[["EH9858"]] imported_example_data_custom <- import_mut_data( mut_file = mut_data, custom_column_names = list(my_contig_name = "contig", my_sample_name = "sample") ) ``` **Data Table 4.** Mutation data imported from tabular file with custom column names. Displays the first 6 rows. Custom column names are changed to default values. ```{r import-mut-cust-table, echo=FALSE} DT::datatable(head(imported_example_data_custom), options = list(scrollX = TRUE)) ``` ## Output Mutation data can be output as either a data frame or a _GRanges_ object for downstream analysis. Use the `output_granges` parameter to specify the output. _GRanges_ may faciliate use in other packages and makes doing genomic based analyses on the ranges significantly easier. Most downstream analyses provided by MutSeqR will use a data frame. # Variant Filtering Following data import, the mutation data may be filtered using the `filter_mut()` function. This function will flag variants based on various parameters in the `filter_mut` column. Variants flagged as TRUE in the `filter_mut` column will be automatically excluded from downstream analyses such as `calculate_mf()`, `plot_bubbles`, and `signature_fitting`. When specified, this function may also remove records from the mutation data. Flagging the variants in the `filter_mut` column will not remove them from the mutation data, however, they will be excluded from mutation counts in all downstream analyses. Filtered variants are retained in the data so that their `total_depth` values may still be used in frequency calculations. By default, all filtering parameters are disabled. Users should be mindful of the filters that they use, ensuring first that they are applicable to their data. ## Germline Variants Germline variants may be identified and flagged for filtering by setting the `vaf_cutoff` parameter. The variant allele fraction (vaf) is the fraction of haploid genomes in the original sample that harbor a specific mutation at a specific base-pair coordinate of the reference genome. Specifically, is it calculated by dividing the number of variant reads by the total sequencing depth at a specific base pair coordinate (`alt_depth` / `total_depth`). In a typical diploid cell, a homozygous germline variant will appear on both alleles, in every cell. As such, we expect this variant to occur on every read, giving us a vaf = 1. A heterozygous germline variant occurs on one of the two alleles in every cell, as such we expect this variant to occur on about half of the reads, giving a vaf = 0.5. Somatic variants occur in only a small portion of the cells, thus we expect them to appear in only a small percentage of the reads. Typical vaf values for somatic variants are less than 0.01 - 0.1. Setting the `vaf_cutoff` parameter will flag all variants that have a vaf **greater than this value as germline** within the `is_germline` column. It will also flag these variants in the `filter_mut` so as to exclude them from downstream analyses. vaf germline filtering is only applicable to users whose data was sequenced to a sufficient depth. High coverage/low depth sequencing cannot be used to calculate the vaf, thus it is recommended to filter germline mutations by contrasting against a database of known polymorphisms, or by using conventional whole-genome sequencing to identify germline variants for each sample. ## Quality Assurance The `filter_mut()` function offers some filtering options to ensure the quality of the mutation data. `snv_in_germ_mnv` = TRUE will flag all snvs that overlap with germline mnvs. Germline mnvs are defined as having a vaf > vaf_cutoff. These snvs are often artifacts from variant calling. No-calls in reads supporting the germline mnv will create false minor-haplotypes from the original mnv that can appear as sub-clonal snvs, thus such variants are excluded from downstream analyses. `rm_abnormal_vaf` = TRUE This parameter identifies rows with abnormal vaf values and **removes** them from the mutation data. Abnormal vaf values are defined as being between 0.05-0.45 or between 0.55-0.95. In a typical diploid organism, we expect variants to have a vaf ~0, 0.5, or 1, reflecting rare somatic variants, heterozygous variants, and homozygous variants respectively. Users should be aware of the ploidy of their model system when using this filter. Non-diploid organisms may exhibit different vafs. `rm_filtered_mut_from_depth = TRUE` Variants flagged in the `filter_mut` column will have their `alt_depth` subtracted from the `total_depth`. When TRUE, this parameter treats flagged variants as No-calls. This does not apply to variants that were idenfied as germline variants. ## Custom Filtering Users may use the `filter_mut()` function to flag or remove variants based on their own custom column. Any record that contains the `custom_filter_val` value within the `custom_filter_col` column of the mutation data will be either flagged in the `filter_mut` column or, if specified by the `custom_filter_rm` parameter, removed from the mutation data. ## Filtering by Regions Users may remove rows that are either within or outside of specified genomic regions. Provide the region ranges to the `regions` parameter. This may be provided as either a filepath, data frame, or a GRanges object. `regions` must contain `contig` (or `seqnames` for GRanges), `start`, and `end`. The function will check whether each record falls within the given regions. Users can define how this filter should be used with `regions_filter`. `region_filter = "remove_within"` will remove all rows whose positions overlap with the provided regions. `region_filter = "keep_within"` will remove all rows whose positions are outside of the provided regions. By default, records that are 1bp+ must start and end within the regions to count as being within the region. `allow_half_overlap = TRUE` allow records that only start or end within the regions but extend outside of them to be counted as being within the region. Twinstrand Mutagenesis Panels may be used by setting `regions` to one of "TSpanel_human", "TSpanel_mouse", or "TSpanel_rat". ## Output The function will return the mutation data as a data.frame object with rows removed or flagged depending on input parameters. Added columns are 'filter_mut', 'filter_reason', 'is_germline', and 'snv_mnv_overlaps'. If `return_filtered_rows = TRUE` The function will return both the data.frame containing the processed mutation data and a data.frame containing the rows that were removed/flagged during the filtering process. The two dataframes will be returned inside of a list, with names `mutation_data` and `filtered_rows`. The function will also print the number of mutations/rows that were filtered according to each parameter. ## Example _Example 2. The example file "Example mutation data" is the output of_ _import_mut_data() run on all 24 mouse libraries from @leblanc-2022._ _Filters used:_ - Filter germline variants: vaf < 0.01 - Filter snvs overlapping with germline variants and have their `alt_depth` removed from their `total_depth`. - Remove records outside of the TwinStrand Mouse Mutagenesis Panel. - Filter variants that contain "EndRepairFillInArtifact" in the "filter" column. Their `alt_depth` will be removed from their `total_depth`. ```{r filter-mut, message=TRUE} # load the example data example_data <- eh[["EH9860"]] # Filter filtered_example_mutation_data <- filter_mut( mutation_data = example_data, vaf_cutoff = 0.01, regions = "TSpanel_mouse", regions_filter = "keep_within", custom_filter_col = "filter", custom_filter_val = "EndRepairFillInArtifact", custom_filter_rm = FALSE, snv_in_germ_mnv = TRUE, rm_filtered_mut_from_depth = TRUE, return_filtered_rows = FALSE ) ``` **Data Table 5.** Filtered Mutation data for all 24 samples. Displayed are 6 variants that were flagged during filtering. ```{r filter-mut-table, echo=FALSE} filter_table <- filtered_example_mutation_data %>% dplyr::filter(filter_mut == TRUE) DT::datatable(head(filter_table), options = list(scrollX = TRUE)) ``` # Calculating MF Mutation Frequency (MF) is calculated by dividing the sum of mutations by the sum of the `total_depth` for a given group (units mutations/bp). The function `calculate_mf()` sums mutation counts and `total_depth` across user-defined groupings within the mutation data to calculate the MF. Mutations can be summarised across samples, experimental groups, and mutation subtypes for later statistical analyses. Variants flagged as `TRUE` in the `filter_mut` column will be excluded from the mutation sums; however, the _total_depth_ of these variants will be counted in the _total_depth_ sum. If the total_depth is not available, the function will sum the mutations across groups, without calculating the mutation frequencies. ## Mutation Counting Mutations are counted based on two opposing assumptions [@dodge-2023]: The **Minimum Independent Mutation Counting Method** (min) Each mutation is counted once, regardless of the number of reads that contain the non-reference allele. This method assumes that multiple instances of the same mutation within a sample are the result of clonal expansion of a single mutational event. When summing mutations across groups using the Min method (`sum_min`), the `alt_depth` of each variant is set to 1. _Ex. For 3 variants with_ _`alt_depth` values of 1, 2, and 10, the `sum_min `= 3._ The **Maximum Independent Mutation Counting Method** (max) Multiple identical mutations at the same position within a sample are counted as independent mutation events. When summing mutations across groups using the Max method (`sum_max`), the `alt_depth` of each variant is summed unchanged. _Ex. for 3_ _variants with `alt_depth` values of 1, 2, and 10, the `sum_max` = 13._ The Min and Max mutation counting methods undercount and overcount the mutations, respectively. We expect some recurrent mutations to be a result of clonal expansion. We also expect some recurrent mutations to arise independently of each other. As we expect the true number of independent mutations to be somewhere in the middle of the two counting methods, we calculate frequencies for both methods. However, the Min mutation counting method is generally recommended for statistical analyses to ensure independance of values and because the Max method tends to increase the sample variance by a significant degree. ## Grouping Mutations Mutation counts and `total_depth` are summed across groups that can be designated using the `cols_to_group` parameter. This parameter can be set to one or more columns in the mutation data that represent experimental variables of interest. _Example 3.1. Calculate mutation sums and frequencies per sample._ _The file "Example mutation data filtered" is the output of filter_mut()_ _from Example 2_ ```{r mf-global} # load example data: example_data <- eh[["EH9861"]] mf_data_global <- calculate_mf( mutation_data = example_data, cols_to_group = "sample", subtype_resolution = "none", retain_metadata_cols = c("dose_group", "dose") ) ``` **Table 6.** Global Mutation Frequency per Sample ```{r mf-global-table, echo=FALSE} #| tbl-cap: "datatable {#tbl-datatable}" DT::datatable(mf_data_global, options = list( columnDefs = list( list(targets = c("mf_min", "mf_max"), render = DT::JS("function(data, type, row, meta) { return data.toExponential(2); }")) ), rowCallback = DT::JS( "function(row, data, dataIndex) { $('td:eq(2)', row).css('text-align', 'right'); // Align the content of column 2 to the right }" ) ) ) ``` Alternatively, you can sum mutations by experimental groups such as `label` (the genomic target). Counts and frequencies will be returned for every level of the designated groups. _Example 3.2. Calculate mutation sums and frequencies per sample and genomic target._ ```{r mf-rg} mf_data_rg <- calculate_mf( mutation_data = example_data, cols_to_group = c("sample", "label"), subtype_resolution = "none", retain_metadata_cols = "dose_group" ) ``` **Table 7.** Mutation Frequency per Sample and genomic target. ```{r mf-rg-table, echo=FALSE} DT::datatable(mf_data_rg, options = list( columnDefs = list( list(targets = c("mf_min", "mf_max"), render = DT::JS("function(data, type, row, meta) { return data.toExponential(2); }")) ), rowCallback = DT::JS( "function(row, data, dataIndex) { $('td:eq(2)', row).css('text-align', 'right'); // Align the content of column 2 to the right }" ) ) ) ``` **calculate_mf() does not calculate the mean MF for any given group** If you want to calculate the mean MF for a given experimental variable, you may group by "sample" and retain the experimental variable in the summary table for averaging using the `retain_metadata_cols` parameter. _Example 3.3. Calculate the mean MF per dose_ ```{r mean-mf} mf_data_global <- calculate_mf( mutation_data = example_data, cols_to_group = "sample", subtype_resolution = "none", retain_metadata_cols = c("dose_group", "dose") ) mean_mf <- mf_data_global %>% dplyr::group_by(dose_group) %>% dplyr::summarise(mean_mf_min = mean(mf_min), SE = sd(mf_min) / sqrt(dplyr::n())) ``` **Table 8.** Mean Mutation Frequency per Dose Group. ```{r mf-mean-table, echo=FALSE} DT::datatable(mean_mf, options = list( columnDefs = list( list(targets = c("mean_mf_min", "SE"), render = DT::JS("function(data, type, row, meta) { return data.toExponential(2); }")) ), rowCallback = DT::JS( "function(row, data, dataIndex) { $('td:eq(2)', row).css('text-align', 'right'); // Align the content of column 2 to the right }" ) ) ) ``` ## Mutation Subtypes Mutations can also be grouped by mutation subtype at varying degrees of resolution using the `subtype_resolution` parameter. Table: (\#tab:subtypes) Definitions of mutation subtype resolutions. | Subtype resolutions | Definition | Example | | ------------------- | ----------------------------------------------------------------------------------------------------------------- | ------------------------------------------------------------- | | `type` | The variation type | snv, mnv, insertion, deletion, complex, and symbolic variants | | `base_6` | The simple spectrum; snv subtypes in their pyrimidine context | C>A, C>G, C>T, T>A, T>C, T>G | | `base_12` | The snv subtypes | A>C, A>G, A>T, C>A, C>G, C>T, G>A, G>C, G>T, T>A, T>C, T>G | | `base_96` | The trinculeotide spectrum; the snv subtypes in their pyrimidine context alongside their two flanking nucleotides | A[C>T]A | | `base_192` | The snv subtypes reported alongside their two flanking nucleotides | A[G>A]A | Mutations and `total_depth` will be summed across groups for each mutation subtype to calculate frequencies. For SNV subtypes, the `total_depth` is summed based on the sequence context in which the SNV subtype occurs (`subtype_depth`). In the simplest example, for the `base_6` SNV subtypes, the two possible reference bases are C or T; hence, the `total_depth` is summed separately for C:G positions and T:A positions. Thus, the MF for C>T mutations is calculated as the total number of C>T mutations divided by total_depth for C:G positions within the group: `sum` / `subtype_depth`. Non-snv mutation subtypes, such as mnvs, insertions, deletions, complex variants, and structural variants, will be calculated as their `sum` / `group_depth`, since they can occur in the context of any nucleotide. Upon import of mutation data, columns are created that facilitate the grouping of SNV subtypes and their associated sequence context for the various resolutions. Below the columns associated with each subtype_resolution are defined: Table: (\#tab:subtypes-cols) Subtype Resolutions and their associated subtype/context columns. | Subtype resolution | Subtype Column | Context Column | | ------------------ | ---------------------------------- | -------------------- | | `type` | `variation_type` | NA | | `base_6` | `normalized_subtype` | `normalized_ref` | | `base_12` | `subtype` | `short_ref` | | `base_96` | `normalized_context_with_mutation` | `normalized_context` | | `base_192` | `context_with_mutation` | `context` | The function will also calculate the proportion of mutations for each subtype, normalized to the `total_depth`: $$P_s = \frac{\left(\frac{M_s}{D_s}\right)}{\sum_s \left(\frac{M_s}{D_s}\right)}$$ Where $P_s$ is the normalized mutation proportion for subtype $s$. $M_s$ is the group mutation sum for subtype $s$. $D_s$ is the group sum of the `subtype_depth` for subtype $s$. If total*depth is not available for the mutation data, `calculate_mf()` will return the subtype mutation counts per group. It will also calculate subtype proportions, without normalizing to the total_depth: $$P'\_s = \frac{M_s}{M*{total}}$$ Where, $P'_s$ is the non-normalized mutation proportion of subtype $s$. $M_s$ is the group mutation sum for subtype $s$. $M_{total}$ is the total mutation sum for the group. _Example 3.4. The following code will return the base_6 mutation spectra for_ _all samples with mutation proportions normalized to depth._ ```{r mf-6} mf_data_6 <- calculate_mf( mutation_data = example_data, cols_to_group = "sample", subtype_resolution = "base_6" ) ``` **Table 9.** Frequency and Proportions of 6-base Mutation Subtypes per Sample. ```{r mf-6-table, echo=FALSE} DT::datatable(mf_data_6, options = list( scrollX = TRUE, columnDefs = list( list(targets = c("mf_min", "mf_max"), render = DT::JS("function(data, type, row, meta) { return data.toExponential(2); }")) ), rowCallback = DT::JS( "function(row, data, dataIndex) { $('td:eq(2)', row).css('text-align', 'right'); // Align the content of column 2 to the right }" ) ) ) ``` ### Selecting Variation Types `calculate_mf()` can be used on a user-defined subset of `variation_type` values. The `variant_types` parameter can be set to a character string of `variation_type` values that the user wants to include in the mutation counts. When calculating group mutation sums, only variants of the specified `variation_type`s will be counted. The `total_depth` for records of excluded `variation_types` will still be included in the `group_depth` and the `subtype_depth`, if applicable. By default the function will calculate summary values based on all mutation types. _Example 3.5. Calculate global mutation frequencies per sample including only_ _insertion and deletion mutations in the count._ ```{r mf-indels} mf_data_global_indels <- calculate_mf( mutation_data = example_data, cols_to_group = "sample", subtype_resolution = "none", variant_types = c("insertion", "deletion") ) ``` **Table 10.** Indel Mutation Frequency per Sample. ```{r mf-indels-table, echo=FALSE} DT::datatable(mf_data_global_indels, options = list( columnDefs = list( list(targets = c("mf_min", "mf_max"), render = DT::JS("function(data, type, row, meta) { return data.toExponential(2); }")) ), rowCallback = DT::JS( "function(row, data, dataIndex) { $('td:eq(2)', row).css('text-align', 'right'); // Align the content of column 2 to the right }" ) ) ) ``` Users may also supply a list of variation_types to exclude to the `variant_types` parameter, so long as the value is preceeded by a "-". _Example 3.6. Calculate mutation frequencies at the "type" subtype resolution,_ _excluding ambiguous and uncategorized mutations._ ```{r mf-types} mf_data_types <- calculate_mf( mutation_data = example_data, cols_to_group = "sample", subtype_resolution = "type", variant_types = c("-ambiguous", "-uncategorized") ) ``` **Table 11.** Frequency and Proportions of Variation Types per Sample. ```{r mf-types-table, echo=FALSE} DT::datatable(mf_data_types, options = list( scrollX = TRUE, columnDefs = list( list(targets = c("mf_min", "mf_max"), render = DT::JS("function(data, type, row, meta) { return data.toExponential(2); }")) ), rowCallback = DT::JS( "function(row, data, dataIndex) { $('td:eq(2)', row).css('text-align', 'right'); // Align the content of column 2 to the right }" ) ) ) ``` _Example 3.7. Include only snv mutations at the base_96 resolution_ ```{r mf-96} mf_data_96 <- calculate_mf( mutation_data = example_data, cols_to_group = "sample", subtype_resolution = "base_96", variant_types = "snv" ) ``` **Table 12.** Frequency and Proportions of 96-base Trinucleotide Subtypes per Sample. ```{r mf-96_tables, echo=FALSE} DT::datatable(mf_data_types, options = list( scrollX = TRUE, columnDefs = list( list(targets = c("mf_min", "mf_max"), render = DT::JS("function(data, type, row, meta) { return data.toExponential(2); }")) ), rowCallback = DT::JS( "function(row, data, dataIndex) { $('td:eq(2)', row).css('text-align', 'right'); // Align the content of column 2 to the right }" ) ) ) ``` ## Correcting Depth There may be cases when the same genomic position is represented multiple times within the mutation data of a single sample. For instance, we require that multiple different alternate alleles for a single position be reported as seperate rows. Another common example of this that we have observed are deletion calls. Often, for each deletion called, a no_variant call at the same start position will also be present. For data that includes the `total_depth` of each record, we want to prevent double-counting the depth at these positions during when summing the depth across groups. When set to TRUE, `correct_depth` parameter will correct the depth at these positions. For positions with the same sample, contig, and start values, the `total_depth` will be retained for only one row. All other rows in the group will have their `total_depth` set to 0 before beign summed. The import functions will automatically check for duplicated rows and return a warning advising to correct the depth should it find any in the mutation data. This is recommended for all users whose data contain `total_depth` values. By default, correct_depth is set to TRUE. _NOTE: The case of deletions and no_variants: When identifying a deletion,_ _some variant callers will re-align the sequences. As such a second_ _total_depth value will be calculated, specifically for the deletion. The_ _variant caller will then provide both the no_variant and the deletion call,_ _the former with the initial depth, and the latter with the re-aligned depth_ _value. Generally, when correcting the depth, the function will retain the_ _total_depth of the first row in the group, and set the rest to 0. However, in_ _the case of deletions, the function may prioritize retaining the re-aligned_ _total_depth over the total_depth assigned to the no_variant. This_ _prioritization can be activated with the `correct_depth_by_indel_priority`_ _parameter._ ## Precalculated Depth For mutation data that does not include a `total_depth` value for each sequenced site, precalculated depths for the specified groups can be supplied separately to `calculate_mf()` through the `precalc_depth_data` parameter. Depth should be provided at the correct subtype resolution for each level of the specified grouping variable. This parameter will accept either a data frame or a filepath to import the depth data. Required columns are: 1. The user-defined cols_to_group variable(s) 2. `group_depth`: the `total_depth` summed across the cols_to_group. 3. The context column for the specified `subtype_resolution`. Only applicable if using the SNV resolutions (base_6, base_12, base_96, base_192). Column names are listed in the table above. 4. `subtype_depth`: the `total_depth` summed across the cols_to_group for each context. Only applicable if using the SNV resolutions. _Example 3.8. Use precalculated depth values to calculate the global per sample MF._ **Table 13. ** 'sample_depth' - Precalculated Depth Data per sample. ```{r precalc1-depth} sample_depth <- data.frame( sample = unique(example_data$sample), group_depth = c(565395266, 755574283, 639909215, 675090988, 598104021, 611295330, 648531765, 713240735, 669734626, 684951248, 716913381, 692323218, 297661400, 172863681, 672259724, 740901132, 558051386, 733727643, 703349287, 884821671, 743311822, 799605045, 677693752, 701163532) ) DT::datatable(sample_depth) ``` ```{r mf-precalc-global} mf_data_global_precalc <- calculate_mf( mutation_data = example_data, cols_to_group = "sample", subtype_resolution = "none", calculate_depth = FALSE, precalc_depth_data = sample_depth ) ``` **Table 14.** Mutation Frequency per Sample, Precalculated Depth Example. ```{r mf-precalc1-table, echo=FALSE} DT::datatable(mf_data_global_precalc, options = list( columnDefs = list( list(targets = c("mf_min", "mf_max"), render = DT::JS("function(data, type, row, meta) { return data.toExponential(2); }")) ), rowCallback = DT::JS( "function(row, data, dataIndex) { $('td:eq(2)', row).css('text-align', 'right'); // Align the content of column 2 to the right }" ) ) ) ``` _Example 3.9. Use precalculated depth values to calculate the base_6 per sample MF._ ```{r mf-precalc-6} mf_data_6_precalc <- calculate_mf( mutation_data = example_data, cols_to_group = "sample", subtype_resolution = "base_6", calculate_depth = FALSE, precalc_depth = eh[["EH9862"]] ) ``` **Table 16.** Precalculated Depth Data per sample for base_6 subtype resolution. ```{r precalc2-depth, echo=FALSE} depth <- eh[["EH9862"]] DT::datatable(depth) ``` **Table 17.** Frequency and Proportions of 6-base Subtypes per Sample, Precalculated Depth Example. ```{r mf-precalc2-table, echo=FALSE} DT::datatable(mf_data_6_precalc, options = list( scrollX = TRUE, columnDefs = list( list(targets = c("mf_min", "mf_max"), render = DT::JS("function(data, type, row, meta) { return data.toExponential(2); }")) ), rowCallback = DT::JS( "function(row, data, dataIndex) { $('td:eq(2)', row).css('text-align', 'right'); // Align the content of column 2 to the right }" ) ) ) ``` View examples for base_12, base_96, and base_192 below: ```{r pre-calc-other-ex} base_12 <- eh[["EH9863"]] base_96 <- eh[["EH9864"]] base_192 <- eh[["EH9865"]] ``` ## Summary Table The function will output the resulting `mf_data` as a data frame with the MF and proportion calculated. If the `summary` parameter is set to `TRUE`, the data frame will be a summary table with the MF calculated for each group. If `summary` is set to `FALSE`, the MF will be appended to each row of the original `mutation_data`. The summary table will include: - `cols_to_group`: all columns used to group the data. - Subtype column for the given resolution, if applicable. - Context column for the given resolution, if applicable. - `sum_min` & `sum_max`: the min/max mutation counts for the group(s)/subtypes. - `group_depth`: the total_depth summed across the group(s). - `subtype_depth`: the total_depth summed across the group(s), for the given context, if applicable. - `MF_min` & `MF_max`: the min/max MF for the group(s)/subtype. - `proportion_min` & `proportion_max`: the min/max subtype proportion, if applicable. Additional columns from the orginal mutation data can be retained using the `retain_metadata_cols` parameter. Retaining higher-order experimental groups may be useful for later statistical analyses or plotting. See above for example _calculating mean MF per dose_. **NOTE** 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. ## MF Plots {.tabset} You can visualize the results of `calculate_mf` using the `plot_mf()` and `plot_mean_mf`() functions. These functions offer variable aesthetics for easy visualization. The output is a ggplot object that can be modified using ggplot2. ### plot_mf _Example 3.10. Plot the Min and Max MF per sample, coloured and ordered by_ _dose group. See example 3.1 for calculating mf_data_global_ ```{r plot-mf-caption, include=FALSE} caption <- paste( "Mutation Frequency (MF) Minimum and Maximum (mutations/bp) per Sample.", "Light colored bars represent MFmin and dark coloured bars represent MFmax.", "Bars are coloured and grouped by dose.", "Data labels are the number of mutations per sample." ) ``` ```{r plot-mf, fig.width=12, fig.height=6, fig.cap=caption} # Define the order for dose groups mf_data_global$dose_group <- factor( mf_data_global$dose_group, levels = c("Control", "Low", "Medium", "High") ) plot <- plot_mf( mf_data = mf_data_global, group_col = "sample", plot_type = "bar", mf_type = "both", fill_col = "dose_group", group_order = "arranged", group_order_input = "dose_group" ) plot ``` ### plot_mean_mf() Calculate and plot the mean MF for a user-defined group using `plot_mean_mf()`. _Example 3.11. Plot the mean MF min per dose, including SEM_ _and individual values coloured by dose. See example 3.1 for calculating_ _mf_data_global_ ```{r plot-mean-mf-caption, include=FALSE} caption <- paste( "Mean Mutation Frequency (MF) Minimum per Dose. Lines are mean ± S.E.M.", "Points are individual samples, coloured by dose." ) ``` ```{r plot-mean-mf, fig.cap=caption} plot_mean <- plot_mean_mf( mf_data = mf_data_global, group_col = "dose_group", mf_type = "min", fill_col = "dose_group", add_labels = "none", group_order = "arranged", group_order_input = "dose_group", plot_legend = FALSE, x_lab = "Dose Group" ) plot_mean ``` # Exporting Results Users can easily output data frames to an Excel workbook with `write_excel()`. This function can write single data frames or it can take a list of dataframes and write each one to a separate Excel sheet in a workbook. In addition to data frames, `write_excel()` will also extract the mf_data, point_estimates, and pairwise_comparisons from `model_mf()` output to write to an excel workbook. Set `model_results` to TRUE if supplying the function with the output to model_mf(). See Modelling Mutation Frequencies for more information on using `model_mf()`. _Example 4.1. Write MF data to excel workbook._ ```{r write-excel, eval=FALSE} # save a single data frame to an Excel file write_excel(mf_data_global, workbook_name = "example_mf_data") # Write multiple data frames to a list to export all at once. list <- list(mf_data_global, mf_data_rg, mf_data_6) names(list) <- c("mf_per_sample", "mf_per_region", "mf_6spectra") #save a list of data frames to an Excel file write_excel(list, workbook_name = "example_mf_data_list") ``` Mutation data can be written to a VCF file for downstream applications with `write_vcf_from_mut()`. _Examples 4.2: export our example data as a vcf file._ ```{r write-vcf, eval=FALSE} write_vcf_from_mut(example_data) ``` # Summary Report To facilitate reproducible analysis of error-corrected sequencing data, the MutSeqR package includes a pre-built R Markdown workflow: Summary_Report.Rmd. This report template allows users to efficiently evaluate the effects of an experimental variable (e.g., dose, age, tissue) on mutation frequency and mutational spectra. The workflow guides users through a standardized analysis pipeline that includes: - Mutation data Import - Mutation data Filtering - Calculation of Global Mutation Frequencies per Sample - Visualization of mutation frequencies per Sample and Mean mutation frequencies per experimental group. - Generalized Linear Modeling of Mutation Frequency (min) as an effect of the experimental variable. - Optional Benchmark Dose Modeling of mutation frequency across a dose or concentration variable. - Calculation of Mutation Subtype Frequencies and Proportions at the base-6 and base-96 resolutions. - Comparison of base-6 mutation spectra between experimental groups. - Visualization of base-6 subtype proportions per sample with hierarchical clustering - Visualization of base-6 subtype proportions per experimental group. - Visualization of base-96 subtype proportions per experimental group. - Mutation signature analysis per experimental group. - Generalized Linear Mixed Modeling of mutation frequency as an effect of specified genomic regions and the experimental group. - Visualization of multiplet mutations per experimental group. All figures will be included in the report. Data frames included in the report will also be exported to excel workbooks within a predefined output path. ## General Usage Users must download the config file "summary_config.yaml". This file defines the parameters to be passed to the Summary_Report.rmd. Users must fill out the .yaml file and save it. ```{r download-yaml, eval=FALSE} config <- system.file("extdata", "inputs", "summary_config.yaml", package = "MutSeqR") file.copy(from = config, to = "path/to/save/the/summary_config.yaml") ``` Users can then render the report using MutSeqR's render_report(). Provide the function with the .yaml file using the `config_file_path` parameter. Provide the name and the file path of the output file to the `output_file` parameter. Set the `output_format` to one of "html_document" (recommended), "pdf_document", or "all". ```{r render-report, eval=FALSE} render_report(config_file_path = "path/to/summary_config.yaml", output_file = "path/to/output_file.html", output_format = "html_document") ``` ## Parameters Below is a description of the parameters set in summary_config.yaml. Parameters are catagerised by their function. System Set Up - _projectdir_: The working directory for the project. All filepaths will be imported from the projectdir. If NULL, the current working directory will be used. - _outputdir_: The output directory where output files will be saved. If NULL, the projectdir will be used. Project Details - _project_title_: The project title. - _research_name_: The name of the principle investigator for the project. - _user_name_: The name of the one running the analysis. - _project_description_: Optional description of the project. Workflow Profile - _config_profile_ : Choose a workflow profile to load pre-defined parameters. Current options are "Duplex Sequencing Mouse Mutagenesis Panel", "Duplex Sequencing Human Mutagenesis Panel", "Duplex Sequencing Rat Mutagenesis Panel", "CODEC"^1, and "None". Pre-defined parameters are set to handle data import and filtering. More profiles coming soon. If not using a pre-defined profile ("None"), users must fill out the Custom_Profile_Params (see below). ^1 CODEC .maf files output from the CODECsuite analysis pipeline may be imported provided the user 1. Adds a "sample" column, 2. Adds an "end" column equal to the Start_Position + length(Reference_Allele) - 1, 3. Ensure Chromsome entries are consistent with notation of BSgenome seqnames for the given reference genome. Ex. GRCh38 requires users to remove "chr" from the Chromosome value. Data Import - _BS_genome_: The BS genome to use to populate the sequence context. The given BS genome must be installed prior to running the render_report() function. To identify an appropriate BS genome based on species and genome assembly, use find_BS_genome(). Ex. "BSgenome.Mmusculus.UCSC.mm10" - _file_type_: The format of the mutation data to be imported. Options are table or vcf. - _mutation_file_: The file path to the mutation data file(s). This can be an individual file or a directory containing multiple files. The file path will be read as projectdir/mutation_file. - _mut_sep_: The delimiter for importing tabular mutation data. Ex. "\t" is tab-delimited. - _sample_data_: An optional file path to the sample metadata that will be joined with the mutation data. Set to NULL if not importing metadata. The file path will be read as projectdir/sample_data. - _sd_sep_: The delimiter for importing the sample_data file. Ex. "\t" is tab-delimited. Calculating MF Precalculated-depth files: The Summary_Report will automatically check during import whether the mutation data contains a depth column (total_depth or depth). If a depth column does not exist, per-sample precalculated depth can be provided for any subtype resolution to calculate mutation frequencies. If an exp_variable is provided, the Summary_Report will automatically sum the per-sample depths to obtain the depth per experimental group as needed. _In order for MutSeqR to calculate the depth from the mutation data, the_ _data must have the depth-value for **every sequenced site**. It is not_ _sufficient to calculate the depth from a list of variants. No-variant sites_ _must be included. The Summary Report will check for a total_depth column, but it will_ _not check whether no_variant sites are included. If it finds a total_depth_ _column, it will proceed to calculate the depth from the mutation data_. - _precalc_depth_data_global_: Optional file path to the precalculated per-sample total_depth data. This is the total number of bases sequenced per sample, used for calculating mutation frequencies. Columns are "sample" and "group_depth". If using an exp_variable (see below), please also include it in this table. The file path will be read as projectdir/precalc_depth_data_global. - _precalc_depth_data_base6_: Optional file path to the precalculated per-sample total_depth data in the base_6 context. This is the total number of C and T bases sequenced for each sample. Columns are "sample", "normalized_ref", and "subtype_depth". If using an exp_variable (see below), please also include it in this table. The file path will be read as projectdir/precalc_depth_data_base6. - _precalc_depth_data_base96_: Optional file path to the precalculated per-sample total_depth data in the base_96 context. This is the total number bases sequenced per sample for each of the 32 possible trinucleotide contexts in their pyrimidine notation. Columns are "sample", "normalized_context", and "subtype_depth". If using an exp_variable, please also include it in this table. The file path will be read as projectdir/precalc_depth_data_base96. - _precalc_depth_data_rg_: Optional file path to the precalculated per-sample total_depth data for each target region. This is the total number of bases sequenced per sample for each region. Columns are "sample", 'region_col', "group_depth". Only applicable if performing regions analysis. The file path will be read as projectdir/precalc_depth_data_rg. - _d_sep_: The delimiter for importing the precalc_depth_data files. Statistical Analysis The Summary Report will run several analyses to investigate the effect of a specified experimental variable on mutation frequency and spectra. - _exp_variable_: Optional. The experimental variable of interest. This argument should refer to a column in your mutation data or sample_data. The Summary Report does not support the analysis of multiple exp_variables; supply only one column name. Ex. "dose". If NULL, statistical analyses will be skipped. - _exp_variable_order_: A vector specifying the unique levels of the exp_variable in the order desired for plotting. - _reference_level_: The reference level of the exp_variable. Ex. the vehicle control group for a chemical dose. This value must match one of the levels within the exp_variable column. If your experimental variable does not have a reference level, then any level of the variable is fine. - _contrasts_: An optional file path to the contrasts table specifying pairwise comparisons between levels of the experimental variable. Required 2 columns, no header. The level in the first column will be compared to the level in the second column for each row. P-values will be adjusted for multiple comparisons if more than one comparison is specified. This contrasts table will be used for Generalized Linear Modeling of MFmin by the exp_variable (model_mf), Generalized Linear Mixed Modeling of MFmin by the exp_variable and target regions, if using a targeted panel (model_mf), and comparing the base_6 mutation spectra between exp_variable levels (spectra_comparison). The file path will be read as projectdir/contrasts. If NULL, spectra_comparison will be skipped, and the GLM/GLMM will return only the model estimates per exp_variable. - _cont_sep_: the delimiter for importing the contrasts file. - _bmd_: A logical variable indicating whether to run a BMD analysis on the MF data. If TRUE, the exp_variable must refer to numeric dose or concentration values. - _bmr_: A numeric value indicating the benchmark response for the bmd analysis. The Summary Report defines the bmr as a bmr-% increase in MF from the reference level. Default is 0.5. Custom Profile Parameters Users who do not wish to use one of the pre-defined parameter profiles may still run the Summary Report by filling out the Custom_Profile_Params. If using one of the config_profiles, users should skip this section. - _ecs_technology_: The technology used to generate the data. - _is_0_based_mut_: A logical variable indicating whether the genomic coordinates of the tabular mutation data are 0-based (TRUE), or 1-based (FALSE). - _regions_: An optional file path to the regions metadata file. The file must contain columns "contig", "start", "end" in addition to the metadata. - _is_0_based_rg_: A logical variable indicating whether the genomic coordinates of the regions file are 0-based (TRUE), or 1-based (FALSE). - _rg_sep_: The delimiter for importing the regions file. - _vaf_cutoff_: A numeric value from 0-1. Will flag variants with a VAF > vaf_cutoff as germline variants and filter them from mutation counts. - _snv_in_germ_mnv_: A logical variable indicating whether filter_mut should flag SNVs that overlap with germline MNVs for filtering. - _rm_abnormal_vaf_: A logical variable indicaring whether filter_mut should remove records where the VAF is between 0.05-0.45 or between 0.55-0.95. - _custom_filter_col_: A column name used by filter_mut to apply a custom filter to the specified column. - _custom_filter_val_: The value in custom_filter_col by which to filter. - _custom_filter_rm_: A logical variable indicating whether records that contain the custom_filter_val within the custom_filter_col should be removed (TRUE) or flagged (FALSE). - _filtering_regions_: Optional file path to regions file by which to filter variants. Must contain the contig, start, and end of each region. The file path will be read as projectdir/filtering_regions. - _filtering_rg_sep_: The delimiter for importing the filtering_regions file. - _regions_filter_: "keep_within" will remove any records that fall outside of the filtering_regions. "remove_within" will remove records that fall inside of the filtering_regions. - _allow_half_overlap_: A logical variable indicating whether to include records that half overlap with the regions. If FALSE, the start and end position of the record must fall within the region interval to be counted as "falling in the region". If TRUE, records that start/end within the region interval, but extend outside of it will be counted as "falling inside the region". - _filtering_is_0_based_rg_: A logical variable indicating whether the genomic coordinates of the filtering_regions are 0-based (TRUE) or 1-based (FALSE). - _rm_filtered_mut_from_depth_: A logical variable indicating whether the alt_depth of variants flagged during the filtering process should be removed from their total_depth values. This does not apply to records flagged as germline variants. - _do_regions_analysis_: A logical variable indicating whether to perform Generalized Linear Mixed Modeling of the data by sequencing target. This is applicable to data sets that used a targeted panel or has specific regions of interest. - _region_col_: The column name that uniquely identifies each target region for the regions_analysis. This column must be present in the mutation_data or in the regions metadata file. ## Output render_report() will output the summary report of results as an html file or a pdf file. All output will be saved to the output_directory specified in the summary_config.yaml. ### Report The Summary Report consists of the following sections. 1. **Import the Mutation Data:** Statement of number of imported rows and the total number of samples detected. 2. **Filter the Mutation Data**: Summary of applied filters with total counts for the number of records that were flagged/removed. 3. **Calculate Per-Sample Mutation Frequency**: Mutation frequency calculated per sample. Includes plot of mutation frequency per sample and plot of mean mutation frequency per experimental group (given experimental variable was provided). 4. **Generalized Linear Modelling**: Model results of mutation frequency as an effect of the experimental variable. Includes model summary, residual plots, model estimates, pairwise comparisons as specified in "contrasts", and a plot of the model estimates and significant results. 5. **Benchmark Dose Analysis**: BMD results for the experimental variable on mutation frequency. BMD estimates per individual model and model-averaged result. Includes model plots, cleveland plots, and confidence interval plots. 6. **6-Base Spectrum**: Mutation frequency calculated at the base_6 resolution per-sample and averaged across the experimental variable. Includes hierarchical clustering of samples by subtype proportion, spectra comparison between experimental groups, and plot of subtype proportions per experimental group. 7. **96-Base Trinucleotide Spectrum**: Mutation frequency calculated at the base_96 resolution per-sample and averaged across the experimental variable. Includes trinucleotide proportion plots for all experimental groups or by sample if no experimental variable is provided. 8. **Regions Analysis**: Mutation frequency calculated per sample for genomic regions of interest. Includes model results of mutation frequency as an effect of the experimental variable and genomic target. Pairwise comparisons specified in contrasts are conducted for each genomic target individually. If no experimental variable is provided, model results reflect mutation frequency as an effect of genomic target. 9. **Mutation HotSpots**: Statement # of mutations that occur in more than one sample. Also includes recurrent mutations within experimental groups. 10. **Lollipop Plots**: Given genomic regions of interest, plots recurrent mutations across genomic position. 11. **Bubble Plot: Multiple Mutations**: Bubble plot depicting multiplet mutations per experimental group. If an experimental variable is not provided, statistical analysis will be skipped. # References # Appendix ## Session Info ```{r session-info, echo=FALSE} sessionInfo() ```