Title:	Tandem Repeat Analysis by Capillary Electrophoresis
Version:	1.0.0
Description:	A pipeline for short tandem repeat instability analysis from fragment analysis data. Inputs of fsa files or peak tables, and a user supplied metadata data-frame. The package identifies ladders, calls peaks, identifies the modal peaks, calls repeats, then calculates repeat instability metrics (e.g. expansion index from Lee et al. (2010) <doi:10.1186/1752-0509-4-29>).
License:	MIT + file LICENSE
Encoding:	UTF-8
RoxygenNote:	7.3.3
Suggests:	knitr, dplyr, ggplot2, shinytest2, testthat (≥ 3.0.0)
Config/testthat/edition:	3
Imports:	graphics, grDevices, lme4, methods, mgcv, plotly, pracma, seqinr, shiny, stats, utils, yaml
VignetteBuilder:	knitr
Depends:	R (≥ 3.5)
LazyData:	true
URL:	https://zachariahmclean.github.io/trace/
NeedsCompilation:	no
Packaged:	2025-12-22 16:26:20 UTC; gusella_lab
Author:	Zachariah McLean [aut, cre, cph], Kevin Correia [aut], Andrew Jiang [ctb]
Maintainer:	Zachariah McLean <zachariah.louis.mclean@gmail.com>
Repository:	CRAN
Date/Publication:	2025-12-22 17:00:02 UTC

trace: Tandem Repeat Analysis by Capillary Electrophoresis

Description

A pipeline for short tandem repeat instability analysis from fragment analysis data. Inputs of fsa files or peak tables, and a user supplied metadata data-frame. The package identifies ladders, calls peaks, identifies the modal peaks, calls repeats, then calculates repeat instability metrics (e.g. expansion index from Lee et al. (2010) doi:10.1186/1752-0509-4-29).

Author(s)

Maintainer: Zachariah McLean zachariah.louis.mclean@gmail.com (ORCID) [copyright holder]

Authors:

• Kevin Correia

Other contributors:

• Andrew Jiang [contributor]

See Also

Useful links:

• https://zachariahmclean.github.io/trace/

Add Metadata to Fragments List

Description

This function adds metadata information to a list of fragments.

Usage

add_metadata(fragments_list, metadata_data.frame)

Arguments

Details

This function adds specified metadata attributes to each fragment in the list. It matches the unique sample identifiers from the fragments list with those in the metadata data frame. This metadata is only required if using the functionality described below.

There are two key things metadata are required for. First is the grouping of samples (metrics_group_id & metrics_baseline_control) for the calculation of metrics and is used in assign_index_peaks(). For example, specifying a sample where the modal allele is the inherited repeat length (eg a mouse tail sample) or sample(s) at the start of a time-course experiment. This is indicated with a TRUE in the metrics_baseline_control column of the metadata. Samples are then grouped together with the metrics_group_id column of the metadata. Multiple samples can be metrics_baseline_control, which can be helpful for the average repeat gain metric to have a more accurate representation of the average repeat at the start of the experiment.

The second key thing metadata can be used for is corrections in call_repeats(). There are two main correction approaches in call_repeats() that are somewhat related: either 'batch' or 'repeat'. Batch correction is relatively simple and just requires you to link samples across batches to correct batch-batch variation in repeat sizes. However, even though the repeat size that is return will be precise, it will not be accurate and underestimates the real repeat length. By contrast, repeat correction can be used to accurately call repeat lengths (which also corrects the batch effects). However, the repeat correction will only be as good as your sample used to call the repeat length so this is a challenging and advanced feature. You need to use a sample that reliably returns the same peak as the modal peak, or you need to be willing to understand the shape of the distribution and manually validate the repeat length of each batch_sample_id for each run.

Batch correction uses common sample(s) across fragment analysis runs to correct systematic batch effects that occur with repeat-containing amplicons in capillary electrophoresis. There are slight fluctuations of size across runs for amplicons containing repeats that result in systematic differences around 1-3 base pairs. So, if samples are to be analyzed for different runs, the absolute bp size is not comparable unless this batch effect is corrected. This is only relevant when the absolute size of a amplicons are compared for grouping metrics as described above (otherwise instability metrics are all relative and it doesn’t matter that there’s systematic batch effects across runs) or when plotting traces from different runs. This correction can be achieved by running a couple of samples in every fragment analysis run, or having a single run that takes a couple of samples from every run together, thereby linking them. These samples are then indicated in the metadata with batch_run_id (to group samples by fragment analysis run) and batch_sample_id (to enable linking samples across batches).

Finally, samples with known and validated repeat size can be used to accurately call the repeat length (and therefore also correct batch effects) in call_repeats(). Similar to batch correction, batch_run_id (to group samples by fragment analysis run) and batch_sample_id (to enable linking samples across batches) are used, but importantly batch_sample_modal_repeat is also set. The batch_sample_modal_repeat is the validated repeat length of the modal repeat of the sample. This validated repeat length is then used to call the repeat length of the modal repeat for each sample (by each batch_run_id). Importantly, this correction requires you to know with confidence the repeat length of the modal peak of the sample. Therefore it's important that the sample used for repeat correction has a clear and prominent modal peak. If the repeat length is very long, it's common for the modal peak of a sample to change so if you use this feature you're going to have to understand the shape of the distribution of your sample and double check that the correct peak has been called as the modal peak after you have used find_alleles(). If a different peak is selected as the modal peak than usual, you need to go back to the metadata and adjust the repeat size of the size standard (For example, your size standard sample has been validated to have 120 repeats. You run find_alleles() and look at the distribution of peaks and notice that the peak one repeat unit higher is the modal peak this time. Therefore, you're going to need to set the batch_sample_modal_repeat as 121 in the metadata just for that batch_run_id. In the other runs you would keep the batch_sample_modal_repeat as 120.).

Value

This function modifies list of fragments objects in place with metadata added.

Examples

gm_raw <- trace::example_data
metadata <- trace::metadata

test_fragments <- genemapper_table_to_fragments(gm_raw,
  dye_channel = "B",
  min_size_bp = 300
)

trace:::add_metadata(
  fragments_list = test_fragments,
  metadata_data.frame = metadata
)

Assign index peaks

Description

Assign index peaks in preparation for calculation of instability metrics

Usage

assign_index_peaks(
  fragments_list,
  config,
  index_override_dataframe = NULL,
  ...
)

Arguments

Details

A key part of instability metrics is the index peak. This is the repeat length used as the reference peak for relative instability metrics calculations, like expansion index. This is usually the the inherited repeat length of a mouse, or the modal repeat length for the cell line at a starting time point.

If grouped is set to TRUE, this function groups the samples by their metrics_group_id and uses the samples set as metrics_baseline_control to set the index peak. Use add_metadata() to set these variables. This is useful for cases like inferring repeat size of inherited alleles from mouse tail data. If the samples that are going to be used to assign index peak are from different fragment analysis runs, use correction = "batch" in call_repeats() to make sure the systematic differences between runs are corrected and the correct index peak is assigned. If there are multiple samples used as baseline control, the median value will be used to assign index peak to corresponding samples.

For mice, if just a few samples have the inherited repeat signal shorter than the expanded population, you could not worry about this and instead use the index_override_dataframe. This can be used to manually override these assigned index repeat values (irrespective of whether grouped is TRUE or FALSE).

As a final option, the index peak could be manually assigned directly to a fragments class using the internal setter function fragments$set_index_peak().

Value

This function modifies list of fragments objects in place with index repeat added.

Examples

fsa_list <- lapply(cell_line_fsa_list, function(x) x$clone())
config <- load_config()

trace:::add_metadata(
  fsa_list,
  metadata_data.frame = trace::metadata
)

trace:::find_ladders(fsa_list, config, show_progress_bar = FALSE)

trace:::find_fragments(fsa_list, config,
  min_bp_size = 300,
  show_progress_bar = FALSE
)

trace:::find_alleles(
  fsa_list,
  config
)
trace:::call_repeats(
  fsa_list,
  config
)

trace:::assign_index_peaks(
  fsa_list,
  config,
  grouped = TRUE
)

plot_traces(fsa_list[1], xlim = c(100,150))

Calculate Repeat Instability Metrics

Description

This function computes instability metrics from a list of fragments data objects.

Usage

calculate_instability_metrics(
  fragments_list,
  peak_threshold = 0.05,
  window_around_index_peak = c(NA_real_, NA_real_),
  percentile_range = c(0.5, 0.75, 0.9, 0.95),
  repeat_range = c(2, 5, 10, 20),
  index_modal_signal_threshold = NA_real_,
  index_signal_sum_threshold = NA_real_
)

Arguments

Details

Each of the columns in the supplied dataframe are explained below:

General Information

• unique_id: A unique identifier for the sample (usually the fsa file name).

Quality Control

• QC_comments: Quality control comments.

• QC_modal_peak_signal: Quality control status based on the modal peak signal (Low < 500, very low < 100).

• QC_peak_number: Quality control status based on the number of peaks (Low < 20, very low < 10).

• QC_off_scale: Quality control comments for off-scale peaks. Potential peaks that are off-scale are given. However, a caveat is that this could be from any of the channels (ie it could be from the ladder channel but is the same scan as the given repeat).

settings used

• peak_threshold: THe peak_threshold parameter used.

• lower_repeat_threshold: The lower repeat limit based of the index repeat of each sample.

• upper_repeat_threshold: The upper repeat limit based of the index repeat of each sample.

• index_modal_signal_threshold: The index_modal_signal_threshold parameter used.

• index_signal_sum_threshold: The index_signal_sum_threshold parameter used.

General sample metrics

• modal_peak_repeat: The repeat size of the modal peak.

• modal_peak_signal: The signal of the modal peak.

• index_peak_repeat: The repeat size of the index peak (the repeat value closest to the modal peak of the index sample).

• index_weighted_mean_repeat: The weighted mean repeat size (weighted on the signal of the peaks) of the index sample.

• n_peaks_total: The total number of peaks in the repeat table.

• n_peaks_analysis_subset: The number of peaks in the analysis subset.

• n_peaks_analysis_subset_expansions: The number of expansion peaks in the analysis subset.

• min_repeat: The minimum repeat size in the analysis subset.

• max_repeat: The maximum repeat size in the analysis subset.

• mean_repeat: The mean repeat size in the analysis subset.

• weighted_mean_repeat: The weighted mean repeat size (weight on peak signal) in the analysis subset.

• median_repeat: The median repeat size in the analysis subset.

• max_signal: The maximum peak signal in the analysis subset.

• sum_signal: The sum of the peak signal in the analysis subset.

• max_delta_neg: The maximum negative delta to the index peak.

• max_delta_pos: The maximum positive delta to the index peak.

• skewness: The skewness of the repeat size distribution.

• kurtosis: The kurtosis of the repeat size distribution.

Repeat instability metrics

• modal_repeat_change: The difference between the modal repeat and the index repeat.

• average_repeat_change: The weighted mean of the sample (weighted by peak signal) subtracted by the weighted mean repeat of the index sample(s).

• instability_index_change: The instability index of the sample subtracted by the instability index of the index sample(s). This will be very similar to the average_repeat_change, with the key difference of instability_index_change being that it is an internally calculated metric for each sample, and therefore the random slight fluctuations of bp size (or systematic if across plates for example) will be removed. However, it requires the index peak to be correctly set for each sample, and if set incorrectly, can produce large arbitrary differences.

• instability_index: The instability index based on peak signal and distance to the index peak. (See Lee et al., 2010, doi:10.1186/1752-0509-4-29).

• instability_index_abs: The absolute instability index. The absolute value is taken for the "Change from the main allele".

• expansion_index: The instability index for expansion peaks only.

• contraction_index: The instability index for contraction peaks only.

• expansion_ratio: The ratio of expansion peaks' signals to the main peak signal. Also known as "peak proportional sum" (See Genetic Modifiers of Huntington’s Disease (GeM-HD) Consortium, 2019, doi:10.1016/j.cell.2019.06.036).

• contraction_ratio: The ratio of contraction peaks' signals to the main peak signal.

• ⁠expansion_percentile_*⁠: The repeat size at specified percentiles of the cumulative distribution of expansion peaks.

• ⁠expansion_percentile_for_repeat_*⁠: The percentile rank of specified repeat sizes in the distribution of expansion peaks.

Value

A data.frame with calculated instability metrics for each sample.

Examples

fsa_list <- lapply(cell_line_fsa_list, function(x) x$clone())
# import data with read_fsa() to generate an equivalent list to cell_line_fsa_list
test_fragments <- trace(fsa_list, grouped = TRUE, metadata_data.frame = metadata)

test_metrics_grouped <- calculate_instability_metrics(
  fragments_list = test_fragments,
  peak_threshold = 0.05,
  window_around_index_peak = c(-40, 40)
)

Call Repeats for Fragments

Description

This function calls the repeat lengths for a list of fragments.

Usage

call_repeats(fragments_list, config, ...)

Arguments

Details

This function has a lot of different options features for determining the repeat length of your samples. This includes i) an option to force the peaks to be whole repeat units apart, ii) corrections to correct batch effects or accurately call repeat length by comparing to samples of known length, and iii) algorithms or re-calling the peaks to remove any contaminating peaks or shoulder-peaks.

———— correction ————

There are two main correction approaches that are somewhat related: either 'batch' or 'repeat'. Batch correction is relatively simple and just requires you to link samples across batches to correct batch-batch variation in repeat sizes. However, even though the repeat size that is return will be precise, it will not be accurate and underestimates the real repeat length. By contrast, repeat correction can be used to accurately call repeat lengths (which also corrects the batch effects). However, the repeat correction will only be as good as your sample used to call the repeat length so this is a challenging and advanced feature. You need to use a sample that reliably returns the same peak as the modal peak, or you need to be willing to understand the shape of the distribution and manually validate the repeat length of each batch_sample_id for each run.

• Batch correction uses common sample(s) across fragment analysis runs to correct systematic batch effects that occur with repeat-containing amplicons in capillary electrophoresis. There are slight fluctuations of size across runs for amplicons containing repeats that result in systematic differences around 1-3 base pairs. So, if samples are to be analyzed for different runs, the absolute bp size is not comparable unless this batch effect is corrected. This is only relevant when the absolute size of a amplicons are compared for grouping metrics as described above (otherwise instability metrics are all relative and it doesn’t matter that there’s systematic batch effects across runs) or when plotting traces from different runs. This correction can be achieved by running a couple of samples in every fragment analysis run, or having a single run that takes a couple of samples from every run together, thereby linking them. These samples are then indicated in the metadata with batch_run_id (to group samples by fragment analysis run) and batch_sample_id (to enable linking samples across batches) (see add_metadata()). Use plot_batch_correction_samples() to plot the samples before and after correction to make sure that is has worked as expected.

• Samples with known and validated repeat size can be used to accurately call the repeat length (and therefore also correct batch effects). Similar to batch correction, batch_run_id (to group samples by fragment analysis run) and batch_sample_id (to enable linking samples across batches) are used, but importantly batch_sample_modal_repeat is also set (see add_metadata()). The batch_sample_modal_repeat is the validated repeat length of the modal repeat of the sample. This validated repeat length is then used to call the repeat length of the modal repeat for each sample (by each batch_run_id). Importantly, this correction requires you to know with confidence the repeat length of the modal peak of the sample. Therefore it's important that the sample used for repeat correction has a clear and prominent modal peak. If the repeat length is very long, it's common for the modal peak of a sample to change so if you use this feature you're going to have to understand the shape of the distribution of your sample and double check that the correct peak has been called as the modal peak after you have used find_alleles(). If a different peak is selected as the modal peak than usual, you need to go back to the metadata and adjust the repeat size of the size standard (For example, your size standard sample has been validated to have 120 repeats. You run find_alleles() and look at the distribution of peaks and notice that the peak one repeat unit higher is the modal peak this time. Therefore, you're going to need to set the batch_sample_modal_repeat as 121 in the metadata just for that batch_run_id. In the other runs you would keep the batch_sample_modal_repeat as 120.). For repeat correction, there are several functions to help visualize and summarize the correction:

  • Use plot_batch_correction_samples() to visualize the same sample across different batches. This can be helpful to make sure that the correction has worked the same across different runs.

  • Use plot_repeat_correction_model() to visualize the linear model use to correct repeat length for each batch_run_id. This can be helpful to make sure the supplied repeat length of different samples are lining up within each run.

  • Generate a summary table of the predicted repeat length for each sample and the average residuals using extract_repeat_correction_summary(). This can be helpful to pinpoint the sample(s) that need adjusting.

———— force_whole_repeat_units ————

The force_whole_repeat_units option aims to correct for the systematic underestimation in fragment sizes that occurs in capillary electrophoresis. It is independent to the algorithms described above and can be used in conjunction. It modifies repeat lengths in a way that helps align peaks with the underlying repeat pattern, making the repeat lengths whole units (rather than ~0.9 repeats). The calculated repeat lengths start from the main peak's repeat length and increases in increments of the specified repeat_size in either direction. This option basically enables you to get exactly the same result as expansion_index values calculated from data from Genemapper.

———— force_repeat_pattern ————

This parameter re-calls the peaks based on specified (force_repeat_pattern_size_period) periodicity of the peaks. The main application of this algorithm is to solve the issue of contaminating peaks in the expected regular pattern of peaks. We can use the periodicity to jump between peaks and crack open a window (force_repeat_pattern_size_window) to then pick out the tallest scan in the window.

Value

This function modifies list of fragments objects in place with repeats added.

See Also

find_alleles(), add_metadata(), plot_batch_correction_samples(), plot_repeat_correction_model(), extract_repeat_correction_summary()

Examples

fsa_list <- lapply(cell_line_fsa_list[c(16:19)], function(x) x$clone())
config <- load_config()

trace:::find_ladders(fsa_list, config, show_progress_bar = FALSE)

trace:::find_fragments(
  fsa_list,
  config,
  min_bp_size = 300,
  show_progress_bar = FALSE
)

trace:::find_alleles(fsa_list, config)

trace:::add_metadata(fsa_list,
   metadata[c(16:19), ]
)

# Simple conversion from bp size to repeat size
trace:::call_repeats(
  fsa_list,
  config,
  assay_size_without_repeat = 87,
  repeat_size = 3
)

plot_traces(fsa_list[1], xlim = c(120, 170))

# Use force_whole_repeat_units algorithm to make sure called
# repeats are the exact number of bp apart

trace:::call_repeats(
  fsa_list,
  config,
  force_whole_repeat_units = TRUE,
  assay_size_without_repeat = 87,
  repeat_size = 3
)

plot_traces(fsa_list[1], xlim = c(120, 170))


# apply batch correction
trace:::call_repeats(
  fsa_list,
  config,
  correction = "batch",
  assay_size_without_repeat = 87,
  repeat_size = 3
)

plot_traces(fsa_list[1], xlim = c(120, 170))

# apply repeat correction
trace:::call_repeats(
  fsa_list,
  config,
  correction = "repeat",
  assay_size_without_repeat = 87,
  repeat_size = 3
)

plot_traces(fsa_list[1], xlim = c(120, 170))

#ensure only periodic peaks are called
trace:::call_repeats(
  fsa_list,
  config,
  force_repeat_pattern_size_period = 2.75,
  assay_size_without_repeat = 87,
  repeat_size = 3
)

plot_traces(fsa_list[1], xlim = c(120, 170))

A list of fsa files

Description

A list of fsa files read into R using trace::read_fsa() that for example data

Usage

cell_line_fsa_list

Format

cell_line_fsa_list

A list with 92 elements, each one being the contents of an fsa file:

Source

doi:10.1038/s41467-024-47485-0

example_data

Description

example_data is genemapper output peak table

Usage

example_data

Format

example_data

A genemapper output dataframe

Source

doi:10.1038/s41467-024-47485-0

example_data_repeat_table

Description

example_data_repeat_table is data with repeats called

Usage

example_data_repeat_table

Format

example_data_repeat_table

A genemapper output dataframe

Source

doi:10.1038/s41467-024-47485-0

Extract Modal Peaks

Description

Extracts modal peak information from each sample in a list of fragments.

Usage

extract_alleles(fragments_list)

Arguments

Value

A dataframe containing modal peak information for each sample

Examples

fsa_list <- lapply(cell_line_fsa_list, function(x) x$clone())
# import data with read_fsa() to generate an equivalent list to cell_line_fsa_list
test_fragments <- trace(fsa_list, grouped = TRUE, metadata_data.frame = metadata)

extract_alleles(test_fragments)

Extract All Fragments

Description

Extracts peak data from each sample in a list of fragments.

Usage

extract_fragments(fragments_list)

Arguments

Value

A dataframe containing peak data for each sample

Examples

fsa_list <- lapply(cell_line_fsa_list, function(x) x$clone())
# import data with read_fsa() to generate an equivalent list to cell_line_fsa_list
test_fragments <- trace(fsa_list, grouped = TRUE, metadata_data.frame = metadata)

extract_fragments(test_fragments)

Extract ladder summary

Description

Extract a table summarizing the ladder models

Usage

extract_ladder_summary(fragments_list, sort = FALSE)

Arguments

Details

The ladder peaks are assigned using a custom algorithm that maximizes the fit of detected ladder peaks and given base-pair sizes. This function summarizes the R-squared values of these individual correlations.

Value

a dataframe of ladder quality information

Examples

fsa_list <- lapply(cell_line_fsa_list, function(x) x$clone())
# import data with read_fsa() to generate an equivalent list to cell_line_fsa_list
test_fragments <- trace(fsa_list, grouped = TRUE, metadata_data.frame = metadata)

  extract_ladder_summary(test_fragments, sort = TRUE)

Extract repeat correction summary

Description

Extracts a table summarizing the model used to correct repeat length

Usage

extract_repeat_correction_summary(fragments_list)

Arguments

Details

For each of the samples used for repeat correction, this table pulls out the modal repeat length called by the model (allele_repeat), how far that sample is on average from the linear model in repeat units by finding the average residuals (avg_residual), and the absolute value of the avg_residual (abs_avg_residual)

Value

A data.frame

Examples

fsa_list <- lapply(cell_line_fsa_list, function(x) x$clone())
# import data with read_fsa() to generate an equivalent list to cell_line_fsa_list
test_fragments <- trace(
   fsa_list, 
   grouped = TRUE, 
   metadata_data.frame = metadata, 
   correction = "repeat",
   show_progress_bar = FALSE
)

# finally extract repeat correction summary
extract_repeat_correction_summary(test_fragments)

Extract traces

Description

Extract the raw trace from a list of fragments objects

Usage

extract_trace_table(fragments_list)

Arguments

Value

A dataframe of the raw trace data. Each row representing a single scan.

Examples

fsa_list <- lapply(cell_line_fsa_list, function(x) x$clone())
# import data with read_fsa() to generate an equivalent list to cell_line_fsa_list
test_fragments <- trace(fsa_list, grouped = TRUE, metadata_data.frame = metadata)

extracted_traces <- extract_trace_table(test_fragments)

Find Alleles

Description

This function identifies main allele within each fragment object.

Usage

find_alleles(fragments_list, config, ...)

Arguments

Details

This function finds the main alleles for each fragment in the list by identifying clusters of peaks ("peak regions") with the highest signal intensities. This is based on the idea that PCR amplicons of repeats have clusters of peaks (from somatic mosaicism and PCR artifacts) that help differentiate the main allele of interest from capillary electrophoresis noise/contamination.

If number_of_alleles = 1, the tallest of peaks will be selected as the allele. This means that if your sample has multiple alleles, you have two options i) make sure that your data is subsetted to only include the allele of interest (using min_bp_size in find_fragments() to make sure that the smaller allele is excluded), or ii) setting number_of_alleles = 2, which will pick the two tallest peaks in their respective peak regions and set the main allele as the larger repeat size, and allele_2 as the shorter repeat size. We recommend the subsetting approach since that is far simpler and less likely to fail, and the second option only if you're doing an experiment analysis a large number of human samples where both the normal and expanded allele repeat lengths vary, which makes it very difficult to find a common bp size that excludes the normal allele.

The parameters peak_region_signal_threshold_multiplier and peak_region_size_gap_threshold will only need to be adjusted in rare cases if peaks are not being found for some reason. They influence the criteria for identifying peak regions. peak_region_signal_threshold_multiplier is multiplied to the mean height of all the peaks to create a hight threshold for inclusion into the peak region, so most of the time it's already a very low value and probably only needs to be changed if you have very few peaks. peak_region_size_gap_threshold is the distance between the peaks, either bp size, or repeats if repeats have already been called.

Value

This function modifies list of fragments objects in place with alleles added.

Examples

fsa_list <- lapply(cell_line_fsa_list[1], function(x) x$clone())
config <- load_config()

trace:::find_ladders(fsa_list, config, show_progress_bar = FALSE)

trace:::find_fragments(fsa_list,
  config,
  min_bp_size = 300,
  show_progress_bar = FALSE
)


trace:::find_alleles(
  fsa_list,
  config,
  peak_region_signal_threshold_multiplier = 1
)

Find fragment peaks

Description

Find fragment peaks in continuous trace data and convert to fragments class.

Usage

find_fragments(fragments_list, config, ...)

Arguments

Details

This takes in a list of fragments objects and returns a list of new fragments objects.

This function is basically a wrapper around pracma::findpeaks(). If your amplicon is large, there may be fewer scans that make up individual peak. So for example you may want to set peak_scan_ramp as a smaller value.

If too many and inappropriate peaks are being called, this may also be solved with the different repeat calling algorithms in call_repeats().

Value

a list of fragments objects.

Examples

fsa_list <- lapply(cell_line_fsa_list[1], function(x) x$clone())
config <- load_config()

trace:::find_ladders(fsa_list, config)

trace:::find_fragments(fsa_list,
  config,
  min_bp_size = 300
)


# Manually inspect the ladders
plot_traces(fsa_list,
  show_peaks = TRUE, n_facet_col = 1,
  xlim = c(400, 550), ylim = c(0, 1200)
)

Ladder and bp sizing

Description

Find the ladder peaks in and use that to call bp size

Usage

find_ladders(fragments_list, config, ...)

Arguments

Details

This function takes a list of fragments files (the output from read_fsa) and identifies the ladders in the ladder channel which is used to call the bp size. The output is a list of fragments.

In this package, base pair (bp) sizes are assigned using a generalized additive model (GAM) with cubic regression splines. The model is fit to known ladder fragment sizes and their corresponding scan positions, capturing the relationship between scan number and bp size. Once trained, the model predicts bp sizes for all scans by interpolating between the known ladder points. This approach provides a flexible and accurate assignment of bp sizes, accommodating the slightly non-linear relationship.

Use plot_data_channels() to plot the raw data on the fsa file to identify which channel the ladder and data are in.

Each ladder should be manually inspected to make sure that is has been correctly assigned.

Value

This function modifies list of fragments objects in place with the ladder assigned and base pair calculated.

See Also

plot_data_channels() to plot the raw data in all channels. plot_ladders() to plot the assigned ladder peaks onto the raw ladder signal. fix_ladders_interactive() to fix ladders with incorrectly assigned peaks.

Examples

fsa_list <- lapply(cell_line_fsa_list[1], function(x) x$clone())
config <- load_config()

trace:::find_ladders(fsa_list, config, show_progress_bar = FALSE)

# Manually inspect the ladders
plot_ladders(fsa_list[1])

Fix ladders interactively

Description

An app for fixing ladders

Usage

fix_ladders_interactive(fragment_trace_list)

Arguments

Details

This function helps you fix ladders that are incorrectly assigned. Run fix_ladders_interactive() and provide output from find_ladders. In the app, for each sample, click on line for the incorrect ladder size and drag it to the correct peak.

Once you are satisfied with the ladders for all the broken samples, click the download button to generate a file that has the ladder correction data. Read this file back into R using readRDS, then use fix_ladders_manual() and supply the ladder correction data as ladder_df_list. This allows the manually corrected data to be saved and used within a script so that the correct does not need to be done every time. An example of what you would need to do:

ladder_df_list <- readRDS('path/to/exported/data.rds') test_ladders_fixed <- fix_ladders_manual(test_ladders_broken, ladder_df_list)

Value

interactive shiny app for fixing ladders

See Also

fix_ladders_manual(), find_ladders()

Examples

# to create an example, lets brake one of the ladders
brake_ladder_list <-  list(
   "20230413_A08.fsa" = data.frame(
     size = c(35, 50, 75, 100, 139, 150, 160, 200, 250, 300, 340, 350, 400, 450, 490, 500),
     scan = c(1544, 1621, 1850, 1912, 2143, 2201, 2261, 2506, 2805, 3135, 3380, 3442, 3760, 
              4050, 4284, 4332)
   )
 )

fsa_list <- lapply(cell_line_fsa_list, function(x) x$clone())
# import data with read_fsa() to generate an equivalent list to cell_line_fsa_list
test_fragments <- trace(fsa_list, ladder_df_list = brake_ladder_list)

if (interactive()) {
  fix_ladders_interactive(test_fragments)
}

# once you have corrected your ladders in the app, 
# export the data for incorporation into the script.
# You can then re-import the data and fix ladders as described in the help details.

Fix ladders manually

Description

Manually assign the ladder peaks for samples in a fragments_list

Usage

fix_ladders_manual(
  fragments_list,
  ladder_df_list,
  warning_rsq_threshold = 0.998
)

Arguments

Details

This function returns a fragments list the same length as was supplied. It goes through each sample and either just returns the same fragments if the unique id doesn't match the samples that need the ladder fixed, or if it is one to fix, it will use the supplied dataframe in the ladder_df_list as the ladder. It then reruns the bp sizing methods on those samples.

This is best used with fix_ladders_interactive() that can generate a ladder_df_list.

Value

This function modifies list of fragments objects in place with the selected ladders fixed.

Examples

config <- load_config()
fsa_list <- lapply(cell_line_fsa_list[1], function(x) x$clone())

trace:::find_ladders(fsa_list, config, show_progress_bar = FALSE)

# first manually determine the real ladder peaks using your judgment
# the raw ladder signal can be extracted
raw_ladder <- fsa_list[1]$raw_ladder

# or we can look at the "trace_bp_df" to see a dataframe that includes the scan and ladder signal
raw_ladder_df <- fsa_list[[1]]$trace_bp_df[, c("unique_id", "scan", "ladder_signal")]
plot(raw_ladder_df$scan, raw_ladder_df$ladder_signal)

# once you have figured what sizes align with which peak, make a dataframe. The
# fix_ladders_manual() function takes a list as an input so that multiple ladders
# can be fixed. Each sample would have the the list element name as it's unique id.

example_list <- list(
 "20230413_A07.fsa" = data.frame(
   size = c(100, 139, 150, 160, 200, 250, 300, 340, 350, 400, 450, 490, 500),
   scan = c(1909, 2139, 2198, 2257, 2502, 2802, 3131, 3376, 3438, 3756, 4046, 4280, 4328)
 )
)

trace:::fix_ladders_manual(
  fsa_list,
  example_list
)

fragments object

Description

An R6 Class representing a fragments object.

Details

This is the parent class of both fragments and fragments object. The idea is that shared fields and methods are both inherited from this object, but it is not itself directly used.

Public fields

Methods

Public methods

• fragments$new()

• fragments$print()

• fragments$plot_trace()

• fragments$plot_ladder()

• fragments$plot_data_channels()

• fragments$get_allele_peak()

• fragments$set_allele_peak()

• fragments$get_index_peak()

• fragments$set_index_peak()

• fragments$plot_fragments()

• fragments$clone()

Method new()

initialization function that is not used since the child classes are the main object of this package.

Usage

fragments$new(unique_id, input_method, object)

Arguments

Method print()

A function to print informative information to the console

Usage

fragments$print()

Method plot_trace()

plot the trace data

Usage

fragments$plot_trace(
  show_peaks = TRUE,
  x_axis = NULL,
  ylim = NULL,
  xlim = NULL,
  signal_color_threshold = 0.05,
  plot_title = NULL
)

Arguments

Returns

A base R plot

Method plot_ladder()

plot the ladder data

Usage

fragments$plot_ladder(xlim = NULL, ylim = NULL, plot_title = NULL)

Arguments

Returns

A base R plot

Method plot_data_channels()

plot the raw data channels in the fsa file. It identifies every channel that has "DATA" in its name.

Usage

fragments$plot_data_channels()

Returns

A base R plot

Method get_allele_peak()

This returns a list with the allele information for this object.

Usage

fragments$get_allele_peak()

Method set_allele_peak()

This sets a single allele size/repeat. It searches through the appropriate peak table and finds the closest peak to the value that's provided.

Usage

fragments$set_allele_peak(allele, unit, value)

Arguments

Method get_index_peak()

This returns a list with the index peak information for this object.

Usage

fragments$get_index_peak()

Method set_index_peak()

This sets the index repeat length. It searches through the repeat table and finds the closest peak to the value that's provided.

Usage

fragments$set_index_peak(value)

Arguments

Method plot_fragments()

This plots the peak/repeat table as a histogram

Usage

fragments$plot_fragments(ylim = NULL, xlim = NULL, plot_title = NULL)

Arguments

Method clone()

The objects of this class are cloneable with this method.

Usage

fragments$clone(deep = FALSE)

Arguments

Convert Genemapper Peak Table to fragments class

Description

This function converts a genemapper peak table data frame into a list of fragments objects.

Usage

genemapper_table_to_fragments(
  df,
  dye_channel,
  min_size_bp = 200,
  max_size_bp = 1000
)

Arguments

Details

This function takes a peak table data frame (eg. Genemapper 5 output) and converts it into a list of fragment objects. It uses the "Sample.File.Name" as the unique id, so make sure that each file as a unique name. Column names should contain: "Dye.Sample.Peak", "Sample.File.Name", "Allele", "Size", "Height".

Value

A list of fragments objects.

See Also

repeat_table_to_fragments(), size_table_to_fragments(), read_fsa()

Examples

gm_raw <- trace::example_data

test_fragments <- genemapper_table_to_fragments(
  gm_raw,
  dye_channel = "B",
  min_size_bp = 400
)

Load Trace Config File

Description

Load configuration file that can be passed to individual trace modules.

Usage

load_config(config_file = NULL)

Arguments

Details

Use this function to load in a config file for passing to individual trace modules. This is only required if creating a custom pipeline rather than the main function. This provides a central place to adjust parameters for the pipeline. Either edit the YAML directly, or modify the object (which is basically just a list).

Use the following command to make a copy of the YAML file: file.copy(system.file("extdata/trace_config.yaml", package = "trace"), ".").

Value

A trace_config object.

Examples

config <- load_config()

metadata

Description

This is a dataframe containing the metadata information for the example data

Usage

metadata

Format

metadata

A genemapper output dataframe

Source

doi:10.1038/s41467-024-47485-0

Plot correction samples

Description

Plot the overlapping traces of the batch control samples

Usage

plot_batch_correction_samples(fragments_list, selected_sample, xlim = NULL)

Arguments

Details

A plot of the raw signal by bp size or repeats for the batch correction samples.

When plotting the traces before repeat correction, we do not expect the samples to be closely overlapping due to run-to-run variation. After repeat correction, the traces should be basically overlapping.

These plots are made using base R plotting. Sometimes these fail to render in the viewing panes of IDEs (eg you get the error 'Error in plot.new(): figure margins too large)'. If this happens, try saving the plot as a pdf using traditional approaches (see grDevices::pdf).

Value

plot of batch corrected samples

See Also

call_repeats() for more info on batch correction.

Examples

fsa_list <- lapply(cell_line_fsa_list, function(x) x$clone())
# import data with read_fsa() to generate an equivalent list to cell_line_fsa_list
fragments_list <- trace(fsa_list, metadata_data.frame = metadata, correction = "batch")


# traces of bp size shows traces at different sizes
plot_batch_correction_samples(
  fragments_list,
  selected_sample = "S-21-212", xlim = c(100, 120)
)

plot_data_channels

Description

Plot the raw data from the fsa file

Usage

plot_data_channels(fragments_list, sample_subset = NULL, n_facet_col = 1)

Arguments

Details

A plot of the raw data channels in the fsa file.

These plots are made using base R plotting. Sometimes these fail to render in the viewing panes of IDEs (eg you get the error 'Error in plot.new(): figure margins too large)'. If this happens, try saving the plot as a pdf using traditional approaches (see grDevices::pdf). To get it to render in the IDE pane, trying matching n_facet_col to the number of samples you're attempting to plot, or using sample_subset to limit it to a single sample.

Value

a plot of the raw data channels

Examples

plot_data_channels(cell_line_fsa_list[1])

Plot Peak Data

Description

Plots peak data from a list of fragments.

Usage

plot_fragments(
  fragments_list,
  n_facet_col = 1,
  sample_subset = NULL,
  xlim = NULL,
  ylim = NULL
)

Arguments

Value

A plot object displaying the peak data.

Examples

fsa_list <- lapply(cell_line_fsa_list, function(x) x$clone())
# import data with read_fsa() to generate an equivalent list to cell_line_fsa_list
fragments_list <- trace(fsa_list)

plot_fragments(fragments_list[1])

Plot ladder

Description

Plot the ladder signal

Usage

plot_ladders(
  fragments_list,
  n_facet_col = 1,
  sample_subset = NULL,
  xlim = NULL,
  ylim = NULL
)

Arguments

Value

a plot of ladders

Examples

fsa_list <- lapply(cell_line_fsa_list, function(x) x$clone())
# import data with read_fsa() to generate an equivalent list to cell_line_fsa_list
fragments_list <- trace(fsa_list)

# Manually inspect the ladders
plot_ladders(fragments_list[1])

Plot Repeat Correction Model

Description

Plots the results of the repeat correction model for a list of fragments.

Usage

plot_repeat_correction_model(
  fragments_list,
  batch_run_id_subset = NULL,
  n_facet_col = 1
)

Arguments

Details

This function makes plots for the model used to correct samples for each batch_run_id. The repeat correction algorithm assigns the user supplied repeat length to the modal peak of the sample, then pulls out a set of robust neighboring peaks to help get enough data to build an accurate linear model for the relationship between base-pair size and repeat length. So on this plot, each dot is an individual peak, with the colour indicating each sample, with the y-axis is the repeat length called from the user-supplied value in the metadata and the value assigned to each peak, with the x-axis showing the corresponding base-pair size.

Value

A base R graphic object displaying the repeat correction model results.

Examples

fsa_list <- lapply(cell_line_fsa_list, function(x) x$clone())
# import data with read_fsa() to generate an equivalent list to cell_line_fsa_list
fragments_list <- trace(fsa_list, metadata_data.frame = metadata, correction = "repeat")

# traces of bp size shows traces at different sizes
plot_repeat_correction_model(
  fragments_list,
  batch_run_id_subset = "20230414"
)

Plot sample traces

Description

Plot the raw trace data

Usage

plot_traces(
  fragments_list,
  show_peaks = TRUE,
  n_facet_col = 1,
  sample_subset = NULL,
  xlim = NULL,
  ylim = NULL,
  x_axis = NULL,
  signal_color_threshold = 0.05
)

Arguments

Details

A plot of the raw signal by bp size. Red vertical line indicates the scan was flagged as off-scale. This is in any channel, so use your best judgment to determine if it's from the sample or ladder channel.

If peaks are called, green is the tallest peak, blue is peaks above the signal threshold (default 5%), purple is below the signal threshold. If force_whole_repeat_units is used within call_repeats(), the called repeat will be connected to the peak in the trace with a horizontal dashed line.

The index peak will be plotted as a vertical dashed line when it has been set using assign_index_peaks().

Value

plot traces from fragments object

Examples

fsa_list <- lapply(cell_line_fsa_list, function(x) x$clone())
# import data with read_fsa() to generate an equivalent list to cell_line_fsa_list
fragments_list <- trace(fsa_list)

plot_traces(fragments_list, xlim = c(105, 150))

Read fsa file

Description

Read fsa file into memory and create fragments object

Usage

read_fsa(files)

Arguments

Details

read_fsa is just a wrapper around seqinr::read.abif() that reads the fsa file into memory and stores it inside a fragments object. That enables you to use the next function find_ladders().

Value

A list of fragments objects

See Also

find_ladders(), plot_data_channels()

Examples

fsa_file <- read_fsa(system.file("abif/2_FAC321_0000205983_B02_004.fsa", package = "seqinr"))
plot_data_channels(fsa_file)

Remove Samples from List

Description

A convenient function to remove specific samples from a list of fragments.

Usage

remove_fragments(fragments_list, samples_to_remove)

Arguments

Value

A modified list of fragments with the specified samples removed.

Examples

gm_raw <- trace::example_data

test_fragments <- genemapper_table_to_fragments(
  gm_raw,
  dye_channel = "B",
  min_size_bp = 300
)

all_fragment_names <- names(test_fragments)

# pull out unique ids of samples to remove
samples_to_remove <- all_fragment_names[c(1, 5, 10)]

samples_removed <- remove_fragments(test_fragments, samples_to_remove)

Convert Repeat Table to Repeats Fragments

Description

This function converts a repeat table data frame into a list of fragments. class.

Usage

repeat_table_to_fragments(df, min_repeat = 0, max_repeat = 1000)

Arguments

Details

This function takes a repeat table data frame and converts it into a list of fragments objects. The column names must be "unique_id" (unique sample id), "repeats" (number of repeats), and "signal" (the signal associated with the number of repeats. for example a frequency count of reads.). The dataframe should be long (eg bind rows if the data are separate).

Value

A list of fragments objects.

Examples

repeat_table <- trace::example_data_repeat_table
test_fragments <- repeat_table_to_fragments(repeat_table)

Convert Size Table to Fragments

Description

This function converts a size table data frame into a list of fragments. class.

Usage

size_table_to_fragments(df, min_size_bp = 200, max_size_bp = 1000)

Arguments

Details

This function takes a size table data frame and converts it into a list of fragments objects. The column names must be "unique_id" (unique sample id), "size" (base pair size), and "signal" (the signal associated with fragment). The dataframe should be long (eg bind rows if the data are separate).

Value

A list of fragments objects.

See Also

repeat_table_to_fragments(), size_table_to_fragments(), read_fsa()

Examples

size_table <- trace::example_data
colnames(size_table)[c(2, 5, 6)] <- c("unique_id", "size", "signal")
fragments_list <- size_table_to_fragments(size_table)

Main function for sample processing

Description

The main function for the trace package that handles processing of samples through the pipeline ready for the calculation of repeat instability metrics.

Usage

trace(
  fragments_list,
  metadata_data.frame = NULL,
  index_override_dataframe = NULL,
  ladder_df_list = NULL,
  config_file = NULL,
  ...
)

Arguments

Details

This function processes samples through the full pipeline, applying the library of functions within this package. Parameters can be adjusted either by passing them directly to this function or by editing a configuration file and supplying it to config_file. Below is a breakdown of the pipeline stages and their key functionalities:

Ladder-related parameters: The ladder is used to calibrate the base pair (bp) sizes of fragments. The ladder peaks are identified in the ladder channel using find_ladders, and a generalized additive model (GAM) with cubic regression splines is used to fit the relationship between scan numbers and bp sizes. This allows for accurate interpolation of bp sizes for all scans. Manual inspection of ladder assignments is recommended to ensure correctness. If ladders are broken, first try and adjust parameters. In a last resort, ladders may be manually set by using the ladder_df_list parameter. To generate this list, use the helper shiny app fix_ladders_interactive.

Peak-finding parameters: Fragment peaks are identified in the continuous trace data using find_fragments. The signal is smoothed using a Savitzky-Golay filter, and peaks are detected based on their signal intensity and bp size. Parameters such as min_bp_size and max_bp_size allow filtering peaks outside the desired range, while peak_scan_ramp controls the number of scans around the peak maxima.

Allele-calling parameters: The main alleles within each fragment are identified by clustering peaks into "peak regions" using find_alleles. The tallest peak in each region is selected as the main allele. If number_of_alleles is set to 2, the two tallest peaks in their respective regions are selected, with the larger repeat size assigned as the main allele.

Repeat-calling parameters: Repeat lengths are calculated using call_repeats, based on the assay size without the repeat and the specified repeat size. Options include batch correction to account for run-to-run variability and repeat correction using samples with known modal repeat lengths. The force_whole_repeat_units option ensures repeat lengths are whole numbers, while force_repeat_pattern re-calls peaks to fit the expected repeat pattern.

Index peak assignment parameters: The index peak is the reference repeat length used for instability metrics calculations. It is typically the inherited repeat length or the modal repeat length at a starting time point. Samples can be grouped to share a common index peak using assign_index_peaks, or the index peak can be manually overridden using index_override_dataframe.

Metadata: Metadata is used to group samples for metrics calculations, batch correction, or repeat length validation. Use add_metadata to add metadata to the fragments list. Required (even if the column is left blank) metadata fields include:

• batch_run_id: Groups samples by fragment analysis run.

• batch_sample_id: Links samples across batches for batch or repeat correction.

• batch_sample_modal_repeat: Specifies the validated repeat length for samples used in repeat correction.

• metrics_group_id: Groups samples for shared index peak assignment.

• metrics_baseline_control: Identifies samples used as baseline controls for index peak assignment.

Value

A list of fragments objects ready for calculation of instability metrics using calculate_instability_metrics()

Examples

fsa_list <- lapply(cell_line_fsa_list, function(x) x$clone())

# import data with read_fsa() to generate an equivalent list to cell_line_fsa_list
fragments_list <- trace(fsa_list, grouped = TRUE, metadata_data.frame = metadata)

metrics <- calculate_instability_metrics(
  fragments_list,
  peak_threshold = 0.05,
  window_around_index_peak = c(-40, 40),
  percentile_range = c(0.5, 0.75, 0.9, 0.95),
  repeat_range = c(2, 5, 10, 20)
)