Package 'TmCalculator'

Description

result$df is computed on demand from result$gr instead of being stored eagerly, so genome-scale results are not converted to a data.frame unless actually requested. All other elements behave as plain list access.

Usage

## S3 method for class 'TmCalculator'
x$name

Arguments

Value

The element, or a data.frame built from x$gr when name == "df" and no stored df exists.

Description

This function processes nucleotide sequences by converting characters to uppercase and replacing invalid bases with "". based on the specified method. The function preserves the sequence length and attributes (name and Tm) of each sequence.

Usage

check_filter_seq(seq_list, method)

Arguments

Value

Returns a list of sequences with the same structure as input, where invalid bases are replaced with ""

Author(s)

Junhui Li

References

citation("TmCalculator")

Description

Apply corrections to melting temperature calculations based on the presence of DMSO and formamide. These corrections are rough approximations and should be used with caution.

Usage

chem_correct(
  DMSO = 0,
  formamide_unit = list(value = 0, unit = "percent"),
  dmso_factor = 0.75,
  formamide_factor = 0.65,
  pt_gc = NULL
)

Arguments

Details

When formamide_unit$unit = "percent": Correction = - factor * percentage_of_formamide

When formamide_unit$unit = "molar": Correction = (0.453 * GC/100 - 2.88) * formamide, from Blake & Delcourt (1996).

Author(s)

Junhui Li

References

von Ahsen N, Wittwer CT, Schutz E. Oligonucleotide melting temperatures under PCR conditions: nearest-neighbor corrections for Mg2+, deoxynucleotide triphosphate, and dimethyl sulfoxide concentrations with comparison to alternative empirical formulas. Clin Chem 2001, 47:1956-1961.

Blake RD, Delcourt SG (1996). Thermodynamic effects of formamide on DNA stability. NAR 24:2095-2103. <doi:10.1093/nar/24.11.2095>

Hutton JR (1977). Renaturation kinetics and thermal stability of DNA in aqueous solutions of formamide and urea. <doi:10.1093/nar/4.10.3537>

Examples

# DMSO correction
chem_correct(DMSO = 3)

# Formamide correction (percent)
chem_correct(formamide_unit = list(value = 1.25, unit = "percent"), pt_gc = 50)

# Formamide correction (molar)
chem_correct(formamide_unit = list(value = 1.25, unit = "molar"), pt_gc = 50)

Description

Test whether one or more numeric metadata columns (e.g. Tm, GC) differ between levels of a categorical grouping column in a GRanges object.

Usage

compare_groups(
  gr,
  target = "Tm",
  method = c("wilcoxon", "t.test"),
  group,
  min_n_per_group = 2L,
  paired = FALSE,
  alternative = c("two.sided", "less", "greater"),
  posthoc = TRUE,
  p.adjust.method = c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr",
    "none")
)

Arguments

Details

The test used depends on method and the number of group levels:

When there are three or more groups and posthoc = TRUE, pairwise follow-up tests are run (Wilcoxon or Welch $t$-test) with p.adjust multiple-testing correction.

Value

A list with:

results

Data frame with one row per target: omnibus test name, statistic, p.value, and group information.

summary

Per-group descriptive statistics (n, mean, sd, median).

pairwise

Pairwise comparisons when posthoc = TRUE and there are $\ge$ 3 groups; otherwise NULL.

Author(s)

Junhui Li

Examples

## Not run: 
library(GenomicRanges)
set.seed(1)
gr <- GRanges(
  seqnames = rep("chr1", 90),
  ranges = IRanges(start = sort(sample(1e6, 90)), width = 50),
  Tm = c(rnorm(30, 74, 2), rnorm(30, 70, 2), rnorm(30, 66, 2)),
  GC = runif(90, 40, 60)
)
gr$group <- rep(c("high", "mid", "low"), each = 30)

# Two groups
gr2 <- gr[gr$group %in% c("high", "low")]
compare_groups(gr2, target = "Tm", group = "group", posthoc = FALSE)

# Three or more groups (Kruskal-Wallis + pairwise Wilcoxon)
compare_groups(gr, target = c("Tm", "GC"), group = "group",
                     method = "wilcoxon")

# Three or more groups (one-way ANOVA + pairwise t-tests)
compare_groups(gr, target = "Tm", group = "group", method = "t.test")

# Post-hoc with Benjamini-Hochberg FDR control
compare_groups(gr, target = "Tm", group = "group",
                     p.adjust.method = "BH")

## End(Not run)

Description

Uses chartr() instead of character-by-character translation. ~10x faster than the seqinr-style implementation for long sequences.

Usage

complement_fast(seq_str, rev = FALSE)

Arguments

Value

Character string(s) with complement.

Description

Fast conversion of genomic coordinate strings to a GRanges object with reference sequences fetched from installed BSgenome.* packages. Designed for large sliding-window inputs: genome packages are loaded once, coordinate strings are parsed in one pass, and sequences are extracted with vectorized getSeq (or optional chromosome preloading).

Usage

coor_to_genomic_ranges(
  input,
  complement_seq = NULL,
  method = c("vectorized", "preload_chr")
)

Arguments

Details

For genome-wide tiling with thousands of windows, pass coordinates as list(pkg_name = "BSgenome.Hsapiens.UCSC.hg38", seq = ...) so the genome package is loaded once instead of per interval.

Value

A GRanges object with metadata columns:

sequence

Reference sequence for each interval.

complement

Complementary sequence.

GC

GC percentage (0-100) per interval, computed as 100 * (G+C)/(A+C+G+T) so that N is excluded from the denominator.

region_id

Region identifier from the coordinate string.

genome_pkg

BSgenome package name used.

Author(s)

Junhui Li

See Also

to_genomic_ranges, to_genomic_ranges_fast

Examples

## Not run: 
coords <- c(
  "chr1:1000-1199:+:win1",
  "chr1:1200-1399:+:win2"
)
gr <- coor_to_genomic_ranges(
  list(pkg_name = "BSgenome.Hsapiens.UCSC.hg38", seq = coords)
)
gr

## End(Not run)

Description

A named list of genomic feature tables for *Escherichia coli* K-12 MG1655 (NCBI assembly GCF_000005845.2 / ASM584v2, chromosome U00096.3). The data support the genome-wide Tm vignette and plot_genome_track examples by providing independent multi-omics layers that can be overlaid on the same genomic axis alongside Tm profiles computed with tm_calculate.

Usage

ecoli_rep_hotspots

Format

A named list with five elements:

all_peaks_IP_mutH

A data frame (38 rows) of MutL protein ChIP-seq peaks marking mismatch-repair-associated regions (MutL-AR). Columns: chr, start, end, Sample, name, Size.

bins_rep

A data frame (4,642 rows) of tandem-repeat (microsatellite) counts in non-overlapping 1 kb bins. Columns: chr, start, end, count.

bins_cru

A data frame (4,642 rows) of cruciform-forming sequence counts in non-overlapping 1 kb bins. Columns: chr, start, end, count.

ssdna

A data frame (2,636 rows) of single-stranded DNA regions. Columns: chr, start, end, Cells..strand., Region.

bins_gatc

A data frame (4,642 rows) of GATC methylation-site counts in non-overlapping 1 kb bins. Columns: chr, start, end, count.

All coordinate-based tables use chr = "U00096.3" and are compatible with plot_genome_track and compare_groups.

Source

Curated from published *E. coli* K-12 MG1655 multi-omics datasets used in the package vignette vignette("genome_wide_tm_ecoli", package = "TmCalculator").

Examples

data(ecoli_rep_hotspots)
names(ecoli_rep_hotspots)

# MutL-AR peak coordinates
head(ecoli_rep_hotspots$all_peaks_IP_mutH)

# Microsatellite density in 1 kb bins
summary(ecoli_rep_hotspots$bins_rep$count)

Description

This function reads sequences from a FASTA file and converts them to a GenomicRanges object. If named with format ">chr2:1-10:[+|-]:[seq_name]", the name will be parsed into GRanges components.

Usage

fa_to_genomic_ranges(input_seq)

Arguments

Value

A GenomicRanges object containing: - GRanges information (seqnames, ranges, strand) - sequence data from FASTA file - Complementary sequences (if provided) - Names from FASTA headers

Examples

# Example with single FASTA file
input_seq <- system.file("extdata", "example1.fasta", package = "TmCalculator")
gr <- fa_to_genomic_ranges(input_seq)

Description

Calculate G and C content of nucleotide sequences. The function calculates the percentage of G and C bases relative to the total number of A, T, G, and C bases in the sequence.

Usage

gc_content(input_seq, ambiguous = FALSE)

Arguments

Value

Content of G and C as a percentage (range from 0 to 100

Author(s)

Junhui Li

Examples

# Calculate GC content of a DNA sequence
gc_content(c("a","t","c","t","g","g","g","c","c","a","g","t","a"))  # 53.85%

# Calculate GC content including ambiguous bases
gc_content("GCATSWSYK", ambiguous = TRUE)  # 55.56%

Description

Generate the complementary sequence of a nucleic acid sequence, with an option to reverse it.

Usage

generate_complement(input_seq, reverse = FALSE)

Arguments

  input                 5'-A T G C G-3'
  reverse = FALSE          T A C G C     (3' to 5', pairs position by position)
  reverse = TRUE        5'-C G C A T-3'  (the same strand, written 5' to 3')
  

Value

Returns the complementary sequence(s) in the specified direction.

Author(s)

Junhui Li

References

citation("TmCalculator")

Examples

# Plain complement: pairs position by position, reads 3' to 5'
generate_complement("ATGCG", reverse = FALSE)

# Reverse complement: the same strand written 5' to 3'
generate_complement("ATGCG", reverse = TRUE)

Description

Combines the output of tm_calculate (a GRanges object with Tm and GC columns) with a second GRanges carrying arbitrary multi-omic metadata (ChIP-seq peaks, ATAC-seq signal, methylation sites, gene annotations, etc.) using one of four positional strategies:

"overlap"

Each tm range is annotated with the aggregated metadata of all feature ranges it directly overlaps.

"nearest"

Each tm range is annotated with the metadata of its single closest feature range, plus an added distance column.

"window"

Each tm range is expanded symmetrically by window_size bp and annotated with aggregated metadata from all features that fall within the expanded window.

"bin"

The genomic space covered by the data is tiled into equal-width bins. Each bin is annotated with the mean tm / GC of overlapping tm ranges and the aggregated feature values - suitable for joint heatmaps and genome-wide correlation analyses.

For strategies "overlap" and "window", when a single Tm range matches multiple features the default behaviour is to summarise: numeric columns are aggregated via agg_fun (default mean), and categorical columns are collapsed to a comma-separated string of unique values.

Usage

integrate_granges(
  gr_tm,
  gr_features,
  strategy = c("overlap", "nearest", "window", "bin"),
  feature_cols = NULL,
  prefix = "",
  window_size = 1000L,
  bin_size = 1e+06,
  agg_fun = mean,
  weight = c("none", "overlap"),
  report_coverage = FALSE,
  min_overlap = 1L,
  ignore_strand = TRUE,
  keep_unmatched = TRUE,
  distance_col = "distance_to_feature"
)

Arguments

Value

• "overlap", "nearest", "window": A GRanges object with the same ranges as gr_tm (minus unmatched ranges if keep_unmatched = FALSE), with additional metadata columns from gr_features.

• "bin": A new GRanges of genomic bins. Each bin carries Tm_mean, GC_mean (if available), n_tm_ranges, n_features, and one aggregated column per requested feature column.

Aggregating continuous signals

When several features map to the same range, numeric columns are summarised by agg_fun and character columns are joined as their unique values. Which features take part is decided by min_overlap, and how much each one counts is decided by weight. The two are easy to confuse, and the distinction is what determines whether a coverage-like signal is summarised correctly at a boundary.

Let $R_i$ be the $i$-th range of gr_tm, $F_j$ the $j$-th feature, $x_j$ its value, and $w_{ij}$ the number of base pairs $R_i$ and $F_j$ share. Write $m$ for min_overlap and $f$ for agg_fun. The features entering the summary of $R_i$ are those with $w_{ij} \ge m$, and

$v_i = f(\{x_j : w_{ij} \ge m\})$

with weight = "none", or

$v_i = \frac{\sum_j w_{ij} x_j}{\sum_j w_{ij}}$

taken over the same features, with weight = "overlap".

The unweighted form gives a feature that overlaps by one base pair the same influence as one that spans the whole range. This is harmless where a signal is flat and wrong where it steps, which is to say at exon boundaries, peak edges and promoters. Raising min_overlap does not repair it: the threshold is a filter, so a short feature is either counted in full or discarded in full, and the bias changes sign rather than disappearing. The example below shows both failures against a case with a known answer.

The weighted mean normalises by the covered base pairs, not by the width of the range, so a value derived from a quarter of a window is indistinguishable from one derived from all of it. report_coverage = TRUE adds a covered_frac column giving the fraction of each range covered by at least one feature, computed after reduce-ing the features so that overlapping ones are not double counted. Multiply by it to convert a mean over covered bases into a mean over the range.

Weighting is not the default. Enabling it changes numeric output, and existing analyses should stay reproducible unless their author decides otherwise.

Author(s)

Junhui Li

Examples

## Aggregation: a coverage track with a known answer -----------------------
## Two 200 bp windows over a signal that steps sharply inside the first.
##
##   window 1  [  1 .. 200]        window 2  [401 .. 600]
##   depth     [  1 .. 190] = 5
##             [191 .. 400] = 200
##                                           [401 .. 450] = 60
library(GenomicRanges)

win <- GRanges("chr1", IRanges(start = c(1, 401), width = 200),
               Tm = c(70, 72))
cov <- GRanges("chr1", IRanges(start = c(1, 191, 401),
                               end   = c(190, 400, 450)),
               cov = c(5, 200, 60))

## Window 1 truly averages (5 * 190 + 200 * 10) / 200 = 14.75.
integrate_granges(win, cov, strategy = "overlap")$cov
## 102.5  60   the 10 bp feature counts as much as the 190 bp one

integrate_granges(win, cov, strategy = "overlap",
                  weight = "overlap")$cov
## 14.75  60   overlap-weighted mean recovers the true depth

integrate_granges(win, cov, strategy = "overlap", min_overlap = 20L)$cov
## 5      60   the threshold discards the short feature; bias reverses

## Window 2 is only a quarter covered. Weighting cannot show that, because
## it normalises by the covered bases; the coverage column can.
res <- integrate_granges(win, cov, strategy = "overlap",
                         weight = "overlap", report_coverage = TRUE)
res$cov                      # 14.75  60
res$covered_frac             # 1.00   0.25
res$cov * res$covered_frac   # 14.75  15    mean over the whole window

## Not run: 
library(GenomicRanges)

# -- Sample data ----------------------------------------------------------
set.seed(42)
gr_tm <- GRanges(
  seqnames = c(rep("chr1", 60), rep("chr2", 30)),
  ranges   = IRanges(
    start = c(sort(sample(1:249e6, 60)),
              sort(sample(1:243e6, 30))),
    width = sample(50:200, 90, replace = TRUE)
  ),
  Tm = runif(90, 55, 85),
  GC = runif(90, 30, 70)
)

gr_features <- GRanges(
  seqnames = c(rep("chr1", 40), rep("chr2", 20)),
  ranges   = IRanges(
    start = c(sort(sample(1:249e6, 40)),
              sort(sample(1:243e6, 20))),
    width = sample(500:5000, 60, replace = TRUE)
  ),
  score      = runif(60, 0, 100),
  peak_type  = sample(c("narrow", "broad"), 60, replace = TRUE),
  signal     = rnorm(60, 5, 2)
)

# Strategy 1: overlap - annotate Tm ranges with overlapping peak features
res_overlap <- integrate_granges(gr_tm, gr_features,
                                  strategy = "overlap")

# Strategy 2: nearest - every Tm range gets its closest peak + distance
res_nearest <- integrate_granges(gr_tm, gr_features,
                                  strategy = "nearest")
head(res_nearest$distance_to_feature)

# Strategy 3: window - 5 kb window around each probe
res_window <- integrate_granges(gr_tm, gr_features,
                                 strategy = "window", window_size = 5000)

# Strategy 4: bin - 500 kb genome bins with mean Tm and aggregated signal
res_bin <- integrate_granges(gr_tm, gr_features,
                              strategy = "bin", bin_size = 5e5)
as.data.frame(res_bin) |> head()

# Use a subset of feature columns and add a prefix
integrate_granges(gr_tm, gr_features,
                  strategy    = "overlap",
                  feature_cols = c("score", "peak_type"),
                  prefix       = "chip_")

## End(Not run)

Description

Tiles one or more chromosomes from a BSgenome object into overlapping windows and returns a named character vector of coordinate strings in the format "chr:start-end:strand:genome_pkg:region_id". This vector is the primary input for downstream Tm calculation functions (tm_nn, tm_gc, tm_wallace) that accept genomic coordinate strings.

Usage

make_genomiccoord(
  bsgenome,
  chromosomes = NULL,
  window = 200L,
  slide = 50L,
  start = NULL,
  end = NULL,
  strand = "+",
  trim_N = c("ends", "filter", "none"),
  max_N_frac = 0.1,
  N_scan_block = window,
  region_prefix = "region",
  genome_pkg_name = NULL,
  as_vector = TRUE,
  verbose = TRUE
)

Arguments

Value

When as_vector = TRUE (default): a named character vector of coordinate strings, one element per window. Names are the region IDs (region1, region2, ...).

When as_vector = FALSE: a data.frame with columns coord, chr, win_start, win_end, strand, genome, region_id, chr_start_used, chr_end_used.

Coordinate string format

Each element of the returned vector follows the pattern:

  chr1:10001-10200:+:BSgenome.Hsapiens.UCSC.hg38:region1
  chr1:10051-10250:+:BSgenome.Hsapiens.UCSC.hg38:region2
  ...

Fields (colon-separated):

1. chromosome – e.g. chr1

2. start-end – 1-based, inclusive coordinates

3. strand – + or -

4. genome – BSgenome package name (character)

5. region_id – unique label regionN

N-base trimming

Human (and most mammalian) chromosomes begin and end with long stretches of N bases that represent assembly gaps. Windows that consist entirely or predominantly of Ns produce meaningless Tm values. The function offers two trimming strategies controlled by trim_N:

"none"

No trimming. Windows start at start and end at end (or chromosome length).

"ends"

Detect the first and last positions of non-N bases on each chromosome using Biostrings::letterFrequency() on a coarse block scan, then trim start and end to those positions. Efficient: reads the chromosome sequence only once in blocks.

"filter"

Generate windows across the full start-end range but then remove any window whose N-base fraction exceeds max_N_frac (default 0.1, i.e. 10%). More granular than "ends" but slower because it reads each window sequence.

Author(s)

Junhui Li

See Also

tm_nn, tm_gc, tm_wallace, BSgenome, letterFrequency

Examples

## Not run: 
library(BSgenome.Hsapiens.UCSC.hg38)

## -- Basic usage: tile chr1 with 200 bp windows, 50 bp slide ----------
coords <- make_genomiccoord(
  bsgenome    = BSgenome.Hsapiens.UCSC.hg38,
  chromosomes = "chr1",
  window      = 200L,
  slide       = 50L
)
length(coords)      # number of windows on chr1
head(coords, 3)
# region1 "chr1:10001-10200:+:BSgenome.Hsapiens.UCSC.hg38:region1"
# region2 "chr1:10051-10250:+:BSgenome.Hsapiens.UCSC.hg38:region2"
# region3 "chr1:10101-10300:+:BSgenome.Hsapiens.UCSC.hg38:region3"

## -- Non-overlapping tiling (slide == window) -------------------------
coords_nonoverlap <- make_genomiccoord(
  bsgenome    = BSgenome.Hsapiens.UCSC.hg38,
  chromosomes = paste0("chr", 1:22),
  window      = 200L,
  slide       = 200L    # no overlap
)
length(coords_nonoverlap)   # ~15 million windows across autosomes

## -- Custom start/end (e.g. a specific sub-region) --------------------
coords_sub <- make_genomiccoord(
  bsgenome    = BSgenome.Hsapiens.UCSC.hg38,
  chromosomes = "chr1",
  window      = 200L,
  slide       = 50L,
  start       = 1000000L,
  end         = 2000000L
)
length(coords_sub)   # 19,981 windows in 1 Mb region

## -- No N-trimming (use full chromosome length) ------------------------
coords_noN <- make_genomiccoord(
  bsgenome    = BSgenome.Hsapiens.UCSC.hg38,
  chromosomes = "chr1",
  window      = 200L,
  slide       = 50L,
  trim_N      = "none"
)

## -- Per-window N-filtering (removes windows with >10% N) -------------
coords_filt <- make_genomiccoord(
  bsgenome    = BSgenome.Hsapiens.UCSC.hg38,
  chromosomes = "chr1",
  window      = 200L,
  slide       = 50L,
  trim_N      = "filter",
  max_N_frac  = 0.10
)

## -- Get data.frame output for GRanges construction -------------------
df <- make_genomiccoord(
  bsgenome    = BSgenome.Hsapiens.UCSC.hg38,
  chromosomes = "chr1",
  window      = 200L,
  slide       = 50L,
  as_vector   = FALSE
)
gr <- GenomicRanges::GRanges(
  seqnames = df$chr,
  ranges   = IRanges::IRanges(start = df$win_start, end = df$win_end),
  strand   = df$strand
)

## -- Pass directly to tm_nn -------------------------------------------
coords <- make_genomiccoord(
  bsgenome    = BSgenome.Hsapiens.UCSC.hg38,
  chromosomes = "chr1",
  window      = 200L,
  slide       = 200L,
  start       = 1000000L,
  end         = 1010000L
)
tm_results <- tm_nn(coords, Na = 50)

## End(Not run)

Description

Visualise genomic tracks as either a linear karyoploteR plot or a circular plot (base R graphics). Set circular = TRUE to switch to the circular layout; the same track_list works for both views without modification.

Supported track types in each list element:

• type = "rect" : data as GRanges/data.frame, drawn as rectangles

• type = "line" : data as GRanges/data.frame with value_col, drawn as lines

• type = "area" : same as line but filled (linear only)

• type = "region": piled-up regions (linear only)

• type = "coverage": coverage area plot (linear only)

A track with ideogram = TRUE is drawn inside the chromosome bar (linear) or as the outermost ring with a grey background (circular). Only "rect" type is supported for ideogram tracks.

A per-track height field (numeric, default 1) controls relative sizing. In linear mode heights are proportional fractions of the data panel; in circular mode they scale the radial thickness of each ring.

A per-track highlight field (list or list-of-lists with data/col/alpha/border) draws highlight bands within that track only.

Entries with type = "highlight" draw translucent bands spanning all real tracks. These accept: data, col (default "#F1C40F"), alpha (default 0.18), border (default NA), and min.degree (circular only, default 0.4).

Each list element may also contain: name, col, bg.col, border, ylim, lwd, alpha (0–1 transparency), legend_font_col, ideogram, height, highlight.

Usage

plot_genome_track(
  genome_name,
  genome_size = NULL,
  genome = NULL,
  track_list,
  circular = FALSE,
  label = NULL,
  chromosomes = NULL,
  zoom = NULL,
  plot.type = 1,
  plot.params = NULL,
  main = NULL,
  track.gap = 0.01,
  legend.show = TRUE,
  legend.position = NULL,
  legend.cex = 0.6,
  legend.bty = "n",
  legend.border = NA,
  legend.font.col = "black",
  legend.lwd = 1,
  legend.seg.len = -0.5,
  legend.box.lwd = 0.25,
  legend.x.intersp = 1,
  legend.lty = 1,
  axis.cex = NULL,
  title.cex = NULL,
  base.tick.dist = NULL,
  base.tick.units = TRUE,
  start.degree = 34.47,
  track.height = 0.05,
  gap.after = 0,
  cell.padding = c(0, 0, 0, 0),
  track.margin = c(0.005, 0.005),
  circle.margin = c(0.001, 0.001),
  canvas.xlim = c(-1, 1),
  canvas.ylim = c(-1, 1),
  axis.unit = "Mb",
  axis.step = NULL,
  axis.show.unit = FALSE,
  label.column = 4,
  label.niceFacing = TRUE,
  label.cex = 0.6,
  label.side = "inside",
  label.labels_height = 0.01,
  label.connection_height = 0.03,
  label.line_lwd = 0.25
)

Arguments

Value

Invisibly returns the KaryoPlot object (linear) or NULL (circular).

Examples

## Not run: 
data(ecoli_rep_hotspots)

library("BSgenome.Ecoli.NCBI.ASM584v2")
genome_name <- "BSgenome.Ecoli.NCBI.ASM584v2"
chr_name    <- "U00096.3"
genome <- get(genome_name, envir = asNamespace(genome_name))
chr_length  <- length(genome[[chr_name]])
genome_name="BSgenome.Ecoli.NCBI.ASM584v2"
bins_gc <- make_genomiccoord(
  bsgenome    = genome_name,
  chromosomes = chr_name,
  window      = 200L,
  slide       = 200L,
  start       = 1,
  end         = chr_length,
  strand      = "+"
)
input_new <- list(pkg_name = genome_name, seq = bins_gc)
gr_batch <- to_genomic_ranges_fast(input_new)
tm_ASM584v2 <- tm_calculate(
  gr_batch,
  method   = "tm_nn"
)
Tm <- as.data.frame(tm_ASM584v2$gr[, c("Tm", "GC")])

tracks <- list(
  list(type = "rect", data = ecoli_rep_hotspots$all_peaks_IP_mutH,
       col = "#2C3E50", bg.col = "grey", name = "MutL-AR",
       legend_font_col = "#2C3E50", ideogram = TRUE, height = 0.5),
  list(type = "line", data = Tm, value_col = "GC",
       name = "GC content", col = "#4A90E2",
       legend_font_col = "#4A90E2"),
  list(type = "line", data = Tm, value_col = "Tm",
       name = "Melting temp", col = "#E06666",
       legend_font_col = "#E06666", height = 2),
  list(type = "line", data = ecoli_rep_hotspots$bins_rep,
       value_col = "count", name = "Microsatellites", col = "#2ECC71",
       legend_font_col = "#2ECC71"),
  list(type = "line", data = ecoli_rep_hotspots$bins_cru,
       value_col = "count", name = "Cruciform", col = "#3B3E6B",
       legend_font_col = "#3B3E6B"),
  list(data = ecoli_rep_hotspots$ssdna, name = "ssDNA",
       col = "#8E44AD", legend_font_col = "#8E44AD"),
  list(type = "line", data = ecoli_rep_hotspots$bins_gatc,
       value_col = "count", name = "GATC sites", col = "#D35400",
       legend_font_col = "#D35400"),
  list(type = "highlight", data = ecoli_rep_hotspots$all_peaks_IP_mutH,
       col = "#F1C40F", alpha = 0.18)
)

# Circular
plot_genome_track("E. coli", genome_size = 4641652,
                  track_list = tracks, circular = TRUE)

# Linear
plot_genome_track("E. coli", genome_size = 4641652,
                  track_list = tracks)

# Linear zoom (single region)
plot_genome_track("E. coli", genome_size = 4641652,
                  track_list = tracks,
                  zoom = "U00096.3:1000000-2000000")

# Linear zoom (multiple regions — stacked panels)
plot_genome_track("E. coli", genome_size = 4641652,
                  track_list = tracks,
                  zoom = c("U00096.3:1000000-1200000",
                           "U00096.3:3000000-3200000"))

# Circular zoom (multiple regions — concatenated arcs)
plot_genome_track("E. coli", genome_size = 4641652,
                  track_list = tracks, circular = TRUE,
                  zoom = c("U00096.3:1000000-1200000",
                           "U00096.3:3000000-3200000"))

# Circular with canvas panning
plot_genome_track("E. coli", genome_size = 4641652,
                  track_list = tracks, circular = TRUE,
                  canvas.xlim = c(0.5, 1), canvas.ylim = c(0, 1),
                  circle.margin = c(0.05, 0.05))

## End(Not run)

Description

Visualise how melting temperature (and optionally other numeric metadata) differs between region classes such as peak vs. non-peak, mutant vs. wild-type, or any categorical annotation stored in a GRanges object. Designed as the visual companion to compare_groups().

Usage

plot_tm(
  gr,
  group,
  value = "Tm",
  plot_type = c("box", "violin", "rank", "density", "ridgeline", "ecdf", "sina"),
  color_palette = c("viridis", "magma", "plasma", "inferno", "cividis"),
  show_points = FALSE,
  point_size = 1,
  point_alpha = 0.3,
  notch = FALSE,
  add_mean = FALSE,
  facet_by = NULL,
  show_pvalue = FALSE,
  p_adjust_method = "BH",
  title = NULL,
  ylab = "Tm (°C)",
  xlab = NULL,
  show_legend = TRUE
)

Arguments

Details

The function pairs naturally with integrate_granges() and compare_groups():

1. Use integrate_granges() to annotate Tm windows with feature overlaps (e.g. ChIP-seq peaks, repeat classes).

2. Use compare_groups() for the statistical test.

3. Use plot_tm() to visualise the distributional difference.

Value

A ggplot object.

Examples

## Not run: 
library(GenomicRanges)
data(ecoli_rep_hotspots)

library("BSgenome.Ecoli.NCBI.ASM584v2")
genome_name <- "BSgenome.Ecoli.NCBI.ASM584v2"
chr_name    <- "U00096.3"
genome <- get(genome_name, envir = asNamespace(genome_name))
chr_length  <- length(genome[[chr_name]])
genome_name="BSgenome.Ecoli.NCBI.ASM584v2"
bins_gc <- make_genomiccoord(
  bsgenome    = genome_name,
  chromosomes = chr_name,
  window      = 200L,
  slide       = 200L,
  start       = 1,
  end         = chr_length,
  strand      = "+"
)
input_new <- list(pkg_name = genome_name, seq = bins_gc)
gr_batch <- to_genomic_ranges_fast(input_new)
tm_ASM584v2 <- tm_calculate(
  gr_batch,
  method   = "tm_nn"
)
gr_tm <- tm_ASM584v2$gr

# Annotate with MutL-AR peak membership
mutH_peaks <- GRanges(
  seqnames = ecoli_rep_hotspots$all_peaks_IP_mutH$chr,
  ranges   = IRanges(start = ecoli_rep_hotspots$all_peaks_IP_mutH$start,
                     end   = ecoli_rep_hotspots$all_peaks_IP_mutH$end)
)
mutH_peaks$peak_id <- paste0("mutH_", seq_along(mutH_peaks))
gr_annot <- integrate_granges(
  gr_tm = gr_tm, gr_features = mutH_peaks,
  strategy = "overlap", feature_cols = "peak_id",
  keep_unmatched = TRUE
)
gr_annot$in_mutH <- ifelse(is.na(gr_annot$peak_id), "non_peak", "peak")

# Box plot — peak vs non-peak Tm
plot_tm(gr_annot, group = "in_mutH")

# Box plot with Wilcoxon p-value bracket
plot_tm(gr_annot, group = "in_mutH", show_pvalue = TRUE)

# Violin plot with p-value and jittered points
plot_tm(gr_annot, group = "in_mutH", plot_type = "violin",
        show_points = TRUE, show_pvalue = TRUE)

# Rank plot — visual Wilcoxon test with p-value subtitle
plot_tm(gr_annot, group = "in_mutH", plot_type = "rank",
        show_pvalue = TRUE)

# Density overlay with p-value
plot_tm(gr_annot, group = "in_mutH", plot_type = "density",
        show_pvalue = TRUE)

# ECDF — cumulative distribution comparison
plot_tm(gr_annot, group = "in_mutH", plot_type = "ecdf",
        show_pvalue = TRUE)

# Compare GC content instead of Tm
plot_tm(gr_annot, group = "in_mutH", value = "GC", ylab = "GC content",
        show_pvalue = TRUE)

# Notched box plot with mean markers and p-value
plot_tm(gr_annot, group = "in_mutH", notch = TRUE, add_mean = TRUE,
        show_pvalue = TRUE)

## End(Not run)

Description

print.TmCalculator prints to console the melting temperature value from an object of class TmCalculator.

Usage

## S3 method for class 'TmCalculator'
print(x, ...)

Arguments

Value

The melting temperature value.

Description

Apply corrections to melting temperature calculations based on salt concentrations. Different correction methods are available for various experimental conditions.

Usage

salt_correct(
  Na = 0,
  K = 0,
  Tris = 0,
  Mg = 0,
  dNTPs = 0,
  method = c("Schildkraut2010", "Wetmur1991", "SantaLucia1996", "SantaLucia1998-1",
    "SantaLucia1998-2", "Owczarzy2004", "Owczarzy2008"),
  input_seq,
  ambiguous = FALSE
)

Arguments

Details

Different correction methods are available for various experimental conditions:

- Schildkraut2010: Schildkraut & Lifson (1965) monovalent correction; the name is historical - Wetmur1991: Classic salt correction method for monovalent cations - SantaLucia1996: DNA-specific salt correction - SantaLucia1998-1: Improved DNA salt correction - SantaLucia1998-2: Alternative DNA salt correction (requires sequence information) - Owczarzy2004: Sodium correction from Owczarzy et al. (2004) (requires sequence information) - Owczarzy2008: Updated comprehensive salt correction (requires sequence information)

For methods other than Owczarzy2008, the sodium-equivalent magnesium adjustment follows von Ahsen et al. (2001); it is separate from the original monovalent formulas.

Owczarzy2008 is piecewise. Which of its three forms applies is decided by the ratio $R = \sqrt{[Mg^{2+}]_{free}} / [Mon]$, where $[Mon]$ is the total monovalent concentration in molar and the free magnesium is what is left after chelation by dNTPs ($K_a = 3 \times 10^4$ $M^{-1}$). Below 0.22 the monovalent term alone applies; between 0.22 and 6 a competing form whose coefficients are themselves functions of $[Mon]$; at 6 and above the divalent-dominated form, in which $[Mon]$ drops out of the expression entirely. A solution with magnesium and no monovalent cation at all therefore has a defined correction rather than none.

Author(s)

Junhui Li

References

Schildkraut C, Lifson S. Dependence of the melting temperature of DNA on salt concentration. Biopolymers. 1965;3(2):195-208.

Wetmur JG. DNA Probes: Applications of the Principles of Nucleic Acid Hybridization. Critical Reviews in Biochemistry and Molecular Biology. 1991;26(3-4):227-259.

SantaLucia J. A unified view of polymer, dumbbell, and oligonucleotide DNA nearest-neighbor thermodynamics. Proceedings of the National Academy of Sciences. 1998;95(4):1460-1465.

Owczarzy R, Moreira BG, Manthey JA, et al. Predicting stability of DNA duplexes in solutions containing magnesium and monovalent cations. Biochemistry. 2008;47(19):5336-5353.

SantaLucia J Jr, Allawi HT, Seneviratne PA (1996). Improved nearest-neighbor parameters for predicting DNA duplex stability. <doi:10.1021/bi951907q>

Owczarzy R et al. (2004). Effects of sodium ions on DNA duplex oligomers: improved predictions of melting temperatures. <doi:10.1021/bi034621r>

von Ahsen N, Wittwer CT, Schutz E (2001). Oligonucleotide melting temperatures under PCR conditions. <doi:10.1093/clinchem/47.11.1956>

Examples

salt_correct(Na = 50, Mg = 1.5, method = "Owczarzy2008", 
               input_seq = "ATGCGATGCG")

Description

A data frame containing coefficients and parameters for different GC-based Tm calculation methods. Each row represents a different method with its specific coefficients (A, B, C, D) and salt correction method.

Usage

thermodynamic_gc_params

Format

A data frame with 8 rows and 5 columns:

A

Intercept coefficient

B

GC content coefficient

C

Length correction coefficient

D

Mismatch coefficient

salt_correct

Associated salt correction method

Details

The methods included are: - Chester1993: Tm = 69.3 + 0.41(Percentage_GC) - 650/N - QuikChange: Tm = 81.5 + 0.41(Percentage_GC) - 675/N - Percentage_mismatch - Schildkraut1965: Tm = 81.5 + 0.41(Percentage_GC) - 675/N + 16.6 x log10[Na+] - Wetmur1991_MELTING: Tm = 81.5 + 0.41(Percentage_GC) - 500/N + Wetmur salt - - Wetmur1991_RNA: Tm = 78 + 0.7(Percentage_GC) - 500/N + Wetmur salt - - Wetmur1991_RNA/DNA: Tm = 67 + 0.8(Percentage_GC) - 500/N + Wetmur salt - - Primer3Plus: Tm = 81.5 + 0.41(Percentage_GC) - 600/N + 16.6 x log10[Na+] - vonAhsen2001: Tm = 77.1 + 0.41(Percentage_GC) - 528/N + 11.7 x log10[Na+]

Description

A comprehensive collection of thermodynamic parameters used for calculating melting temperatures of nucleic acid duplexes. The dataset includes parameters for DNA/DNA, RNA/RNA, RNA/DNA and 2'-O-methylRNA/RNA hybridizations, plus parameters for mismatches, dangling ends, internal / bulge loops, GU wobble, CNG repeats, inosine, hydroxyadenine, azobenzene and locked nucleic acids (LNA).

Usage

thermodynamic_nn_params

Format

A named list of two-column matrices with colnames = c("left", "right") (left = $\Delta H$, right = $\Delta S$). Rownames are dinucleotide step keys such as "AT/TA", mismatch / dangling-end keys, or feature-specific keys (e.g. "CAG_4" for the (CAG)4 CNG entry).

Populated tables

DNA_NN_Breslauer_1986

DNA/DNA NN, Breslauer et al. (1986)

DNA_NN_Sugimoto_1996

DNA/DNA NN, Sugimoto et al. (1996)

DNA_NN_Allawi_1998

DNA/DNA Watson-Crick NN, Allawi & SantaLucia (1997), Table 1; legacy identifier

DNA_NN_SantaLucia_2004

DNA/DNA NN, SantaLucia & Hicks (2004)

RNA_NN_Freier_1986

RNA/RNA NN, Freier (1986)

RNA_NN_Xia_1998

RNA/RNA NN, Xia (1998)

RNA_NN_Chen_2012

RNA/RNA NN with GU, Chen / Serra (2012)

RNA_DNA_NN_Sugimoto_1995

RNA/DNA hybrid NN, Sugimoto (1995)

DNA_IMM_Peyret_1999

DNA internal mismatches, six sources (1997-2005); see tm_nn

DNA_TMM_Bommarito_2000

DNA terminal mismatch, SantaLucia & Peyret (2001), WO2001094611A2, Tables 2-3; legacy identifier

DNA_DE_Bommarito_2000

DNA single dangling end, Bommarito (2000)

RNA_DE_Turner_2010

RNA single dangling end, NNDB Turner 2004 compilation; database described in 2010

Registered placeholders (values TBD - populate from primary literature)

DNA_NN_AllSan_1997

Allawi & SantaLucia (1997) Biochemistry 36:10581

DNA_NN_SantaLucia_1996

SantaLucia et al. (1996) Biochemistry 35:3555

DNA_NN_Tanaka_2004

Tanaka et al. (2004) Biochemistry 43:7143

MeRNA_RNA_NN_Kierzek_2006

Kierzek et al. (2006) Biochemistry 45:581

DNA_RNA_IMM_Watkins_2011

Watkins et al. (2011) Nucleic Acids Res 39:1894

RNA_IMM_Lu_2006

Lu et al. (2006) Nucleic Acids Res 34:4912

RNA_IMM_Davis_Znosko_2007

Davis & Znosko (2007) Biochemistry 46:13425

RNA_IMM_Davis_Znosko_2008

Davis & Znosko (2008) Biochemistry 47:10178

DNA_TandemMM_AllSanPey

Allawi & SantaLucia (1997/1998); Peyret (1999)

RNA_TandemMM_Mathews_1999

Mathews et al. (1999) JMB 288:911

DNA_DE_Ohmichi_2002

Ohmichi et al. (2002) JACS 124:10367

RNA_DE_Ohmichi_2002

Ohmichi et al. (2002) JACS 124:10367

RNA_DE_Serra_2008

O'Toole / Serra (2006, 2008)

DNA_DDE_Ohmichi_2002

Ohmichi et al. (2002) JACS 124:10367

RNA_DDE_Ohmichi_2002

Ohmichi et al. (2002) JACS 124:10367

RNA_DDE_OToole_2005

O'Toole et al. (2005) Biochemistry 44:14914

RNA_DDE_OToole_2006

O'Toole et al. (2006) NAR 34:3338

DNA_LDE_Ohmichi_2002

Ohmichi et al. (2002) JACS 124:10367

RNA_LDE_Ohmichi_2002

Ohmichi et al. (2002) JACS 124:10367

DNA_IL_SantaLucia_2004

SantaLucia & Hicks (2004)

RNA_IL_Lu_2006

Lu et al. (2006) NAR 34:4912

RNA_IL_Badhwar_2007

Badhwar et al. (2007) Biochemistry 46:14715

DNA_SBL_Tanaka_2007

Tan & Chen (2007) Biophys J 92:3615

DNA_SBL_SantaLucia_2004

SantaLucia & Hicks (2004)

RNA_SBL_Lu_2006

Lu et al. (2006) NAR 34:4912

RNA_SBL_Blose_2007

Blose et al. (2007) Biochemistry 46:15123

DNA_LBL_SantaLucia_2004

SantaLucia & Hicks (2004)

RNA_LBL_Lu_2006

Mathews (1999); Lu et al. (2006)

RNA_GU_Mathews_1999

Mathews & Turner (1999) JMB 288:911

RNA_GU_Chen_2012

Chen / Serra (2012)

DNA_CNG_Broda_2005

Broda et al. (2005) Biochemistry 44:13851 - rownames use paste0("C", N, "G_", n), e.g. "CAG_4".

DNA_INO_SantaLucia_2005

Watkins & SantaLucia (2005) NAR 33:6258

RNA_INO_Wright_2007

Wright et al. (2007) Biochemistry 46:4625

DNA_HA_Kawakami_2001

Kawakami / Sugimoto (2001) Biochemistry 40:14040

DNA_AZB_Asanuma_2005

Asanuma et al. (2005) Angew Chem Int Ed 44:2747

DNA_LNA_Owczarzy_2011

Owczarzy et al. (2011) Biochemistry 50:9352

DNA_LNA_McTigue_2004

McTigue et al. (2004) Biochemistry 43:5388

DNA_cLNA_Owczarzy_2011

Owczarzy et al. (2011) Biochemistry 50:9352

DNA_cLNA_MM_Owczarzy_2011

Owczarzy et al. (2011) Biochemistry 50:9352

To add the numeric values for any placeholder table, construct a matrix with two columns ($\Delta H$ kcal/mol, $\Delta S$ cal/(mol*K)) and the dinucleotide-step rownames described above, then assign it:

thermodynamic_nn_params$DNA_CNG_Broda_2005 <- new_table

Details

Coverage mirrors the rmelting (Bioconductor) vignette sections 4.2.1 - 4.4. Tables whose numeric values are not yet ported from primary literature are included in the list as NULL placeholders so that tm_nn can still dispatch on the table name; the corresponding contribution is silently skipped after a one-time package-startup notice.

Source

Various publications as cited above. See the installed extdata/tm_nn_reference_audit.md for the audit of runtime tables.

References

Breslauer K J (1986) <doi:10.1073/pnas.83.11.3746> Sugimoto N (1996) <doi:10.1093/nar/24.22.4501> Allawi H T & SantaLucia J Jr (1997), Table 1 <doi:10.1021/bi962590c>; also SantaLucia (1998), Table 2 <doi:10.1073/pnas.95.4.1460> SantaLucia J (2004) <doi:10.1146/annurev.biophys.32.110601.141800> Freier S (1986) <doi:10.1073/pnas.83.24.9373> Xia T (1998) <doi:10.1021/bi9809425> Chen JL (2012) <doi:10.1021/bi3002709> Sugimoto N (1995) <doi:10.1021/bi00035a029> Bommarito S (2000) <doi:10.1093/nar/28.9.1929> Peyret N (1999) <doi:10.1021/bi9825091> Allawi H T & SantaLucia J (1997) <doi:10.1021/bi962590c> Watkins N E Jr & SantaLucia J Jr (2005) <doi:10.1093/nar/gki918> Allawi H T & SantaLucia J Jr (1998, G.A) <doi:10.1021/bi9724873> Allawi H T & SantaLucia J Jr (1998, A.C) <doi:10.1021/bi9803729> Allawi H T & SantaLucia J Jr (1998, C.T) <doi:10.1093/nar/26.11.2694> SantaLucia J Jr & Peyret N (2001), WO2001094611A2, Tables 2-3. https://patents.google.com/patent/WO2001094611A2/en Turner D H & Mathews D H (2010) <doi:10.1093/nar/gkp892> Tanaka F (2004) Biochemistry 43:7143 Kierzek E (2006) Biochemistry 45:581 Watkins N E (2011) Nucleic Acids Res 39:1894 Lu Z J (2006) Nucleic Acids Res 34:4912 Davis A R & Znosko B M (2007, 2008) Biochemistry Mathews D H (1999) <doi:10.1006/jmbi.1999.2700> Ohmichi T (2002) J Am Chem Soc 124:10367 O'Toole A S (2005, 2006) Biochemistry / Nucleic Acids Res Badhwar J (2007) Biochemistry 46:14715 Tan Z J & Chen S J (2007) Biophys J 92:3615 Blose J M (2007) Biochemistry 46:15123 Broda M (2005) <doi:10.1021/bi0501447> Wright D J (2007) Biochemistry 46:4625 Kawakami J / Sugimoto N (2001) <doi:10.1021/bi010918b> Asanuma H (2005) Angew Chem Int Ed 44:2747 McTigue P M (2004) Biochemistry 43:5388 Owczarzy R (2011) Biochemistry 50:9352

Examples

# Access DNA/DNA nearest neighbor parameters
thermodynamic_nn_params$DNA_NN_SantaLucia_2004

# Access DNA internal mismatch parameters
thermodynamic_nn_params$DNA_IMM_Peyret_1999

# See which tables are still placeholders (NULL)
names(Filter(is.null, thermodynamic_nn_params))

Description

Calculates nucleic acid melting temperature (Tm) by one of three methods, and returns the result as a GRanges object so that Tm can be used directly as a quantitative genomic feature alongside other assays:

• Nearest neighbor (tm_nn) sums the stacking free energies of adjacent base-pair steps using experimentally derived enthalpy and entropy parameters, and so resolves sequences of identical base composition but different order. It supports salt and chemical corrections, mismatches and dangling ends, and is the default.

• GC content (tm_gc) computes Tm as an empirical function of GC percentage with corrections for length and ionic strength. Cheaper than NN, but blind to sequence order.

• Wallace rule (tm_wallace) assigns a fixed contribution per base. It is calibrated for short oligonucleotides, typically 14 to 20 bp, and ignores sequence context and reaction conditions, so it is not appropriate for long sequences or for genome-wide windows.

Usage

tm_calculate(
  input_seq,
  method = c("tm_nn", "tm_gc", "tm_wallace"),
  complement_seq = NULL,
  ambiguous = FALSE,
  shift = 0,
  nn_table = c("DNA_NN_SantaLucia_2004", "DNA_NN_Ghosh_2020_PEG200",
    "DNA_NN_Breslauer_1986", "DNA_NN_Sugimoto_1996", "DNA_NN_Allawi_1998",
    "RNA_NN_Freier_1986", "RNA_NN_Xia_1998", "RNA_NN_Chen_2012", "RNA_NN_Zuber_2022",
    "RNA_NN_Ghosh_2023_PEG200", "RNA_DNA_NN_Sugimoto_1995", "DNA_NN_Weber_2015",
    "DNA_NN_Weber_OW04_69", "DNA_NN_Weber_OW04_119", "DNA_NN_Weber_OW04_220",
    "DNA_NN_Weber_OW04_621", "DNA_NN_Weber_OW04_1020", "RNA_NN_Weber_VIF_71",
    "RNA_NN_Weber_VIF_121", "RNA_NN_Weber_VIF_221", "RNA_NN_Weber_VIF_621", 
    
    "RNA_NN_Weber_VIF_1021", "RNA_NN_Weber_FIF_71", "RNA_NN_Weber_FIF_121",
    "RNA_NN_Weber_FIF_221", "RNA_NN_Weber_FIF_621", "RNA_NN_Weber_FIF_1021",
    "RNA_DNA_NN_Weber_2019_FT", "RNA_DNA_NN_Weber_2019_VH", "RNA_DNA_NN_Weber_2019_LS",
    "RNA_DNA_NN_Banerjee_2020"),
  tmm_table = "DNA_TMM_Bommarito_2000",
  imm_table = "DNA_IMM_Peyret_1999",
  de_table = c("DNA_DE_Bommarito_2000", "RNA_DE_Turner_2010"),
  dnac_high = 25,
  dnac_low = 25,
  self_comp = FALSE,
  variant = c("Primer3Plus", "Chester1993", "QuikChange", "Schildkraut1965",
    "Wetmur1991_MELTING", "Wetmur1991_RNA", "Wetmur1991_RNA/DNA", "vonAhsen2001"),
  userset = NULL,
  Na = 50,
  K = 0,
  Tris = 0,
  Mg = 0,
  dNTPs = 0,
  salt_method = c("Schildkraut2010", "Wetmur1991", "SantaLucia1996", "SantaLucia1998-1",
    "SantaLucia1998-2", "Owczarzy2004", "Owczarzy2008", "none"),
  DMSO = 0,
  formamide_unit = list(value = 0, unit = "percent"),
  dmso_factor = 0.75,
  formamide_factor = 0.65,
  mismatch = TRUE,
  regions = NULL,
  window = NULL,
  slide = window,
  unit = c("segment", "region"),
  segment_size = 5e+07,
  BPPARAM = NULL,
  keep_sequence = NULL,
  tmpdir = tempdir(),
  verbose = FALSE
)

Arguments

Details

The input sequence is processed once and passed to the selected method, which is faster than calling the individual functions separately.

The three methods differ in resolution and in the range of sequence lengths over which they are calibrated, so they are not interchangeable.

tm_nn is the appropriate default. Because it sums sequence-dependent stacking terms, it distinguishes sequences of identical GC content but different base order, which the other two cannot. Its parameters were derived from short duplexes under a two-state assumption; when applied to long sequences or to fixed-width genomic windows the resulting value is best read as a relative measure of local thermodynamic stability for comparison across windows, rather than as an absolute experimental melting temperature.

tm_gc computes Tm from GC percentage using one of several published empirical formulas selected by variant, with corrections for length and ionic strength. It extends to longer sequences at low computational cost but cannot resolve base order.

tm_wallace applies the 2 + 4 rule, adding 2 $^{\circ}$C per A or T and 4 $^{\circ}$C per G or C. It is calibrated for short oligonucleotides, typically 14 to 20 bp, and takes no account of sequence context, salt or chemical additives. Accuracy degrades quickly with length, so it is retained for compatibility rather than recommended for genome-scale work.

Salt and chemical corrections apply to tm_nn and tm_gc only. The input sequence is parsed and validated once and reused by the selected method, which is faster than calling the individual functions directly.

Value

A TmCalculator list with:

Salt handling

Most nearest-neighbor parameter sets were fitted at a single reference sodium concentration, and other conditions are reached through the salt_method correction formulas. The Weber/VarGibbs sets were instead fitted directly at a stated sodium concentration and are meant to replace salt correction. When such a set is selected and Na matches the concentration it was fitted at, correction is skipped automatically; when it does not, correction is applied with a warning. See tm_nn for details.

Available Options

Method Selection:

• method: c("tm_nn", "tm_gc", "tm_wallace")

Nearest Neighbor (NN) Method Options:

• nn_table:

  • DNA/DNA: "DNA_NN_Breslauer_1986", "DNA_NN_Sugimoto_1996", "DNA_NN_Allawi_1998", "DNA_NN_SantaLucia_2004" (default)

  • DNA/DNA, molecular crowding: "DNA_NN_Ghosh_2020_PEG200" (40 wt

  • DNA/DNA, salt-optimized: "DNA_NN_Weber_2015" (1020 mM), "DNA_NN_Weber_OW04_69", "..._119", "..._220", "..._621", "..._1020" (fitted at 69 to 1020 mM sodium)

  • RNA/RNA: "RNA_NN_Freier_1986", "RNA_NN_Xia_1998", "RNA_NN_Chen_2012", "RNA_NN_Zuber_2022" (improved end effects), "RNA_NN_Ghosh_2023_PEG200" (molecular crowding, cell-like)

  • RNA/RNA, salt-optimized: "RNA_NN_Weber_VIF_71", "..._121", "..._221", "..._621", "..._1021" and the corresponding "RNA_NN_Weber_FIF_*" sets

  • RNA/DNA: "RNA_DNA_NN_Sugimoto_1995", "RNA_DNA_NN_Weber_2019_FT", "RNA_DNA_NN_Weber_2019_VH" (1000 mM), "RNA_DNA_NN_Weber_2019_LS" (100 mM), "RNA_DNA_NN_Banerjee_2020" (100 mM)

• tmm_table (Terminal Mismatches):

  • "DNA_TMM_Bommarito_2000" (default)

• imm_table (Internal Mismatches):

  • "DNA_IMM_Peyret_1999" (default)

• de_table (Dangling Ends):

  • "DNA_DE_Bommarito_2000" (default)

  • "RNA_DE_Turner_2010"

GC Method Options:

• variant:

  • "Primer3Plus" (default)

  • "Chester1993"

  • "QuikChange"

  • "Schildkraut1965"

  • "Wetmur1991_MELTING"

  • "Wetmur1991_RNA"

  • "Wetmur1991_RNA/DNA"

  • "vonAhsen2001"

Salt Correction Options:

• salt_method:

  • "Schildkraut2010" (default)

  • "Wetmur1991"

  • "SantaLucia1996"

  • "SantaLucia1998-1"

  • "SantaLucia1998-2" (method = "tm_nn" only)

  • "Owczarzy2004" (method = "tm_nn" only)

  • "Owczarzy2008" (method = "tm_nn" only)

  • "none" (also selected automatically when nn_table was fitted at the requested Na)

With method = "tm_gc" the salt term belongs to the published formula selected by variant, so naming a different one is ignored there, with a warning, unless userset is supplied. "none" and NA drop the correction and are honoured on either path.

Formamide Unit Options:

• formamide_unit$unit:

  • "percent" (default)

  • "molar"

Other Parameters:

• ambiguous: TRUE/FALSE (default: FALSE)

• shift: Integer value (default: 0)

• dnac_high: Numeric value in nM (default: 25)

• dnac_low: Numeric value in nM (default: 25)

• self_comp: TRUE/FALSE (default: FALSE)

• Na: Millimolar concentration (default: 50)

• K: Millimolar concentration (default: 0)

• Tris: Millimolar concentration (default: 0)

• Mg: Millimolar concentration (default: 0)

• dNTPs: Millimolar concentration (default: 0)

• DMSO: Percent concentration (default: 0)

• dmso_factor: Numeric value (default: 0.75)

• formamide_factor: Numeric value (default: 0.65)

• mismatch: TRUE/FALSE (default: TRUE)

Author(s)

Junhui Li

See Also

tm_nn for the nearest-neighbor method and the full list of thermodynamic parameter sets.

Examples

## Not run: 
input_seq <- c("chr1:1000100-1000150:+:BSgenome.Hsapiens.UCSC.hg38")
result <- tm_calculate(
  input_seq,
  method = "tm_nn",
  nn_table = "DNA_NN_SantaLucia_2004",
  salt_method = "Owczarzy2008"
)

# A hybrid parameter set fitted at 100 mM sodium. Salt correction is
# skipped automatically because Na matches the fitted condition.
result_ls <- tm_calculate(
  input_seq,
  method = "tm_nn",
  nn_table = "RNA_DNA_NN_Weber_2019_LS",
  Na = 100
)

# Genome scale. The source is named rather than loaded, because each
# worker opens it for itself; `regions` says what to cover and `window`
# says at what resolution.
hg38 <- "BSgenome.Hsapiens.UCSC.hg38"
tm_calculate(hg38, regions = "chr21:10e6-20e6", window = 200, slide = 200)

# Five processes. Tasks are divided by region, never by splitting the
# sequences of one region, so only coordinates cross between them.
library(BiocParallel)
whole <- tm_calculate(hg38, window = 200, slide = 200,
                      BPPARAM = SnowParam(workers = 5))
whole$gr

# The same, for a FASTA file and for sequences already in R. With
# window = NULL, the default, each record returns a single Tm, which is
# what a file of probes or primers calls for.
tm_calculate("probes.fa", BPPARAM = SnowParam(workers = 4))
tm_calculate(c("ACGTGCTAGCTAGCTAGC", "GGCCATATATGCGC"))

## End(Not run)

Description

Calculate the melting temperature using empirical formulas based on GC content with different options. The function returns a list of sequences with updated Tm attributes and calculation options.

Usage

tm_gc(
  gr_seq,
  ambiguous = FALSE,
  userset = NULL,
  variant = c("Primer3Plus", "Chester1993", "QuikChange", "Schildkraut1965",
    "Wetmur1991_MELTING", "Wetmur1991_RNA", "Wetmur1991_RNA/DNA", "vonAhsen2001"),
  Na = 50,
  K = 0,
  Tris = 0,
  Mg = 0,
  dNTPs = 0,
  salt_method = NULL,
  mismatch = TRUE,
  DMSO = 0,
  formamide_unit = list(value = 0, unit = "percent"),
  dmso_factor = 0.75,
  formamide_factor = 0.65
)

Arguments

Value

Returns a list with two components: - Tm: A list of sequences with updated Tm attributes - Options: A list containing calculation parameters and method information

Author(s)

Junhui Li

References

Marmur J, Doty P. Determination of the base composition of deoxyribonucleic acid from its thermal denaturation temperature. Journal of Molecular Biology, 1962, 5(1):109-118.

Schildkraut C, Lifson S. Dependence of the melting temperature of DNA on salt concentration. Biopolymers, 1965, 3(2):195-208.

Wetmur JG. DNA Probes: Applications of the Principles of Nucleic Acid Hybridization. CRC Critical Reviews in Biochemistry, 1991, 26(3-4):33.

Untergasser A, Cutcutache I, Koressaar T, et al. Primer3–new capabilities and interfaces. Nucleic Acids Research, 2012, 40(15):e115-e115.

von Ahsen N, Wittwer CT, Schutz E, et al. Oligonucleotide melting temperatures under PCR conditions: deoxynucleotide Triphosphate and Dimethyl sulfoxide concentrations with comparison to alternative empirical formulas. Clin Chem 2001, 47:1956-1961.

Examples

# Example with multiple sequences
input_seq <- c("ATCGTGCGTAGCAGTACGATCAGTAG", "ATCGTGCGTAGCAGTACGATCAGTAG")
gr_seq <- to_genomic_ranges(input_seq)
out <- tm_gc(gr_seq, ambiguous = TRUE, variant = "Primer3Plus", Na = 50, mismatch = TRUE)
out
out$options

Description

Calculates melting temperature (Tm) from nearest-neighbor (NN) thermodynamics, summing the stacking enthalpies and entropies of adjacent base-pair steps and applying initiation, symmetry, salt and chemical corrections. Terminal mismatches, internal mismatches and dangling ends are supported through dedicated parameter tables. The function verifies that every dinucleotide step in the input sequence is present in the selected tables before calculating.

Usage

tm_nn(
  gr_seq,
  ambiguous = FALSE,
  shift = 0,
  nn_table = c("DNA_NN_SantaLucia_2004", "DNA_NN_Ghosh_2020_PEG200",
    "DNA_NN_Breslauer_1986", "DNA_NN_Sugimoto_1996", "DNA_NN_Allawi_1998",
    "RNA_NN_Freier_1986", "RNA_NN_Xia_1998", "RNA_NN_Chen_2012", "RNA_NN_Zuber_2022",
    "RNA_NN_Ghosh_2023_PEG200", "RNA_DNA_NN_Sugimoto_1995", "DNA_NN_Weber_2015",
    "DNA_NN_Weber_OW04_69", "DNA_NN_Weber_OW04_119", "DNA_NN_Weber_OW04_220",
    "DNA_NN_Weber_OW04_621", "DNA_NN_Weber_OW04_1020", "RNA_NN_Weber_VIF_71",
    "RNA_NN_Weber_VIF_121", "RNA_NN_Weber_VIF_221", "RNA_NN_Weber_VIF_621", 
    
    "RNA_NN_Weber_VIF_1021", "RNA_NN_Weber_FIF_71", "RNA_NN_Weber_FIF_121",
    "RNA_NN_Weber_FIF_221", "RNA_NN_Weber_FIF_621", "RNA_NN_Weber_FIF_1021",
    "RNA_DNA_NN_Weber_2019_FT", "RNA_DNA_NN_Weber_2019_VH", "RNA_DNA_NN_Weber_2019_LS",
    "RNA_DNA_NN_Banerjee_2020"),
  tmm_table = "DNA_TMM_Bommarito_2000",
  imm_table = "DNA_IMM_Peyret_1999",
  de_table = c("DNA_DE_Bommarito_2000", "RNA_DE_Turner_2010"),
  dnac_high = 25,
  dnac_low = 25,
  self_comp = FALSE,
  Na = 50,
  K = 0,
  Tris = 0,
  Mg = 0,
  dNTPs = 0,
  salt_method = c("Schildkraut2010", "Wetmur1991", "SantaLucia1996", "SantaLucia1998-1",
    "SantaLucia1998-2", "Owczarzy2004", "Owczarzy2008", "none"),
  DMSO = 0,
  formamide_unit = list(value = 0, unit = "percent"),
  dmso_factor = 0.75,
  formamide_factor = 0.65
)

Arguments

Details

DNA_NN_Breslauer_1986: Breslauer K J (1986) <doi:10.1073/pnas.83.11.3746>

DNA_NN_Sugimoto_1996: Sugimoto N (1996) <doi:10.1093/nar/24.22.4501>

DNA_NN_Allawi_1998: Allawi H T & SantaLucia J Jr (1997), Table 1 <doi:10.1021/bi962590c>; also SantaLucia (1998), Table 2 <doi:10.1073/pnas.95.4.1460>

DNA_NN_SantaLucia_2004: SantaLucia J (2004) <doi:10.1146/annurev.biophys.32.110601.141800>

RNA_NN_Freier_1986: Freier S (1986) <doi:10.1073/pnas.83.24.9373>

RNA_NN_Xia_1998: Xia T (1998) <doi:10.1021/bi9809425>

RNA_NN_Chen_2012: Chen JL (2012) <doi:10.1021/bi3002709>

RNA_DNA_NN_Sugimoto_1995: Sugimoto N (1995)<doi:10.1021/bi00035a029>

The following sets were derived by melting-temperature optimization and are fitted at the sodium concentration given in parentheses. They are not salt-corrected further; see the “Choosing a parameter set” section.

DNA_NN_Weber_2015 (1020 mM): Weber G (2015) <doi:10.1093/bioinformatics/btu751>

DNA_NN_Weber_OW04_69 (69 mM), DNA_NN_Weber_OW04_119 (119 mM), DNA_NN_Weber_OW04_220 (220 mM), DNA_NN_Weber_OW04_621 (621 mM), DNA_NN_Weber_OW04_1020 (1020 mM): Weber G (2015) <doi:10.1093/bioinformatics/btu751>

RNA_NN_Weber_VIF_71 (71 mM), RNA_NN_Weber_VIF_121 (121 mM), RNA_NN_Weber_VIF_221 (221 mM), RNA_NN_Weber_VIF_621 (621 mM), RNA_NN_Weber_VIF_1021 (1021 mM): Ferreira I (2019) <doi:10.1016/j.chemphys.2019.01.016>, variable initiation factors

RNA_NN_Weber_FIF_71 (71 mM), RNA_NN_Weber_FIF_121 (121 mM), RNA_NN_Weber_FIF_221 (221 mM), RNA_NN_Weber_FIF_621 (621 mM), RNA_NN_Weber_FIF_1021 (1021 mM): Ferreira I (2019) <doi:10.1016/j.chemphys.2019.01.016>, fixed initiation factors

RNA_DNA_NN_Weber_2019_FT (1000 mM), RNA_DNA_NN_Weber_2019_VH (1000 mM), RNA_DNA_NN_Weber_2019_LS (100 mM): Basilio Barbosa V (2019) <doi:10.1016/j.bpc.2019.106189>

RNA_DNA_NN_Banerjee_2020 (100 mM): Banerjee D (2020) <doi:10.1093/nar/gkaa572>

RNA_NN_Zuber_2022: Zuber J (2022) <doi:10.1093/nar/gkac261>

RNA_NN_Ghosh_2023_PEG200 (100 mM, 40 wt

DNA_NN_Ghosh_2020_PEG200 (100 mM, 40 wt

DNA_TMM_Bommarito_2000: SantaLucia J Jr & Peyret N (2001), WO2001094611A2, Tables 2-3 (https://patents.google.com/patent/WO2001094611A2/en)

DNA_IMM_Peyret_1999: Allawi H T & SantaLucia J Jr (1997, G.T) <doi:10.1021/bi962590c>; (1998, G.A) <doi:10.1021/bi9724873>; (1998, A.C) <doi:10.1021/bi9803729>; (1998, C.T) <doi:10.1093/nar/26.11.2694>; Peyret N et al. (1999, A.A/C.C/G.G/T.T) <doi:10.1021/bi9825091>; Watkins N E Jr & SantaLucia J Jr (2005, inosine) <doi:10.1093/nar/gki918>

DNA_DE_Bommarito_2000: Bommarito S (2000) <doi:10.1093/nar/28.9.1929>

RNA_DE_Turner_2010: NNDB Turner 2004 dangling-end compilation; Turner D H & Mathews D H (2010, database description) <doi:10.1093/nar/gkp892>

Value

A TmCalculator list with:

Choosing a parameter set

Parameter sets fall into two families that are handled differently.

The two families are distinguished by one mechanical criterion: a set is condition-specific if and only if it carries a salt_mM attribute.

Reference-salt sets (no salt_mM; Breslauer 1986, Sugimoto 1996, Allawi 1997 (legacy name: Allawi_1998), SantaLucia 2004, Freier 1986, Xia 1998, Chen 2012, Zuber 2022, Sugimoto 1995) were fitted at a single reference sodium concentration, and other conditions are reached by applying one of the salt_method correction formulas.

Condition-specific sets (the Weber/VarGibbs series, Banerjee 2020, and the molecular-crowding sets of Ghosh 2020 and Ghosh 2023) were instead fitted directly at a stated sodium concentration and are intended to replace salt correction rather than be corrected. When the requested Na matches the set's salt_mM value, salt correction is skipped automatically; when it does not, the correction is still applied but a warning is issued, since correcting an already condition-specific set double-counts the ionic effect. Whether a correction was applied is recorded in the returned options.

As a rule of thumb, pick the set whose fitted salt is closest to your experimental condition rather than correcting a distant one.

What the initiation terms are charged on

A dangling end or a terminal mismatch is accounted for by its own parameter, from de_table or tmm_table, and the position it covers is then consumed. The initiation terms – init, init_A/T, init_G/C, init_5T/A and the init_allA/T / init_oneG/C pair – are charged on what remains, that is, on the duplex that actually closes.

This matters because the terminal AT penalty is a property of the closing base pair. SantaLucia and Hicks (2004) define it as applied for each end of a duplex that has a terminal AT, so the question it asks is whether the pair holding that end shut is AT or GC. A mismatched terminus is neither an AT nor a GC pair, and a dangling residue is not in a base pair at all, so neither one can carry the penalty. Charging it there would also double-count the end, since the terminal-mismatch and dangling-end parameters were measured as terminal and already carry that end's contribution.

Some implementations index these terms on the sequence as supplied rather than on the trimmed duplex. That makes the melting temperature depend on which of the two strands is handed over as gr_seq: the same duplex, written from the other strand, can then come back with a different answer. The convention here is invariant under that swap by construction. See the 1.1.2 entry in NEWS for a worked example.

Parameter provenance and limitations

Table identifiers are retained for compatibility and are not always literal author-year citations. The installed files extdata/tm_nn_reference_audit.md and extdata/tm_nn_parameter_sources.tsv record the source and verification status of every built-in parameter row.

The audit identified unresolved numerical/model issues: the patent's GG/CG terminal-mismatch enthalpy is inconsistent with its printed free energy; RNA_NN_Chen_2012 combines refitted Watson-Crick terms with GU terms fitted against Xia (1998); and Banerjee (2020) initiation conditions in Table 2 footnotes are not reproduced by the current per-end application. The audit corrects citations without changing these numerical results. Consult the audit before using these affected models.

Author(s)

Junhui Li

References

Breslauer KJ, Frank R, Blocker H, Marky LA (1986). Predicting DNA duplex stability from the base sequence. PNAS 83:3746-3750. <doi:10.1073/pnas.83.11.3746>

Sugimoto N, Nakano S, Yoneyama M, Honda K (1996). Improved thermodynamic parameters and helix initiation factor to predict stability of DNA duplexes. NAR 24:4501-4505. <doi:10.1093/nar/24.22.4501>

Allawi HT, SantaLucia J Jr (1997). Thermodynamics and NMR of internal G.T mismatches in DNA. Biochemistry 36:10581-10594. Table 1: Watson-Crick parameters; Table 5: G.T mismatches. <doi:10.1021/bi962590c>

SantaLucia J Jr (1998). A unified view of polymer, dumbbell, and oligonucleotide DNA nearest-neighbor thermodynamics. PNAS 95:1460-1465. <doi:10.1073/pnas.95.4.1460>

SantaLucia J Jr, Hicks D (2004). The thermodynamics of DNA structural motifs. Annual Review of Biophysics and Biomolecular Structure 33:415-440. <doi:10.1146/annurev.biophys.32.110601.141800>

Freier SM et al. (1986). Improved free-energy parameters for predictions of RNA duplex stability. PNAS 83:9373-9377. <doi:10.1073/pnas.83.24.9373>

Xia T et al. (1998). Thermodynamic parameters for an expanded nearest-neighbor model for formation of RNA duplexes with Watson-Crick base pairs. Biochemistry 37:14719-14735. <doi:10.1021/bi9809425>

Chen JL et al. (2012). Testing the nearest neighbor model for canonical RNA base pairs: revision of GU parameters. Biochemistry 51:3508-3522. <doi:10.1021/bi3002709>

Sugimoto N et al. (1995). Thermodynamic parameters to predict stability of RNA/DNA hybrid duplexes. Biochemistry 34:11211-11216. <doi:10.1021/bi00035a029>

Allawi HT, SantaLucia J Jr (1998). Nearest-neighbor thermodynamics of internal A.C mismatches in DNA: sequence dependence and pH effects. Biochemistry 37:9435-9444. <doi:10.1021/bi9803729>

Allawi HT, SantaLucia J Jr (1998). Nearest neighbor thermodynamic parameters for internal G.A mismatches in DNA. Biochemistry 37:2170-2179. <doi:10.1021/bi9724873>

Allawi HT, SantaLucia J Jr (1998). Thermodynamics of internal C.T mismatches in DNA. NAR 26:2694-2701. <doi:10.1093/nar/26.11.2694>

Peyret N, Seneviratne PA, Allawi HT, SantaLucia J Jr (1999). Nearest-neighbor thermodynamics and NMR of DNA sequences with internal A.A, C.C, G.G, and T.T mismatches. Biochemistry 38:3468-3477. <doi:10.1021/bi9825091>

Watkins NE Jr, SantaLucia J Jr (2005). Nearest-neighbor thermodynamics of deoxyinosine pairs in DNA duplexes. NAR 33:6258-6267. <doi:10.1093/nar/gki918>

SantaLucia J Jr, Peyret N (2001). Method and system for predicting nucleic acid hybridization thermodynamics and computer-readable storage medium for use therein. Patent WO2001094611A2, published December 13, 2001, Tables 2-3. https://patents.google.com/patent/WO2001094611A2/en

Bommarito S, Peyret N, SantaLucia J Jr (2000). Thermodynamic parameters for DNA sequences with dangling ends. NAR 28:1929-1934. <doi:10.1093/nar/28.9.1929>

Turner DH, Mathews DH (2010). NNDB: the nearest neighbor parameter database for predicting stability of nucleic acid secondary structure. NAR 38:D280-D282. <doi:10.1093/nar/gkp892> The RNA dangling-end values are the Turner 2004 compilation: https://rna.urmc.rochester.edu/NNDB/rna_2004/rna_2004_dangling_ends.html

Weber G (2015). Optimization method for obtaining nearest-neighbour DNA entropies and enthalpies directly from melting temperatures. Bioinformatics 31:871-877. <doi:10.1093/bioinformatics/btu751>

Ferreira I, Jolley EA, Znosko BM, Weber G (2019). Replacing salt correction factors with optimized RNA nearest-neighbour enthalpy and entropy parameters. Chemical Physics 521:69-76. <doi:10.1016/j.chemphys.2019.01.016>

Basilio Barbosa V, de Oliveira Martins E, Weber G (2019). Nearest-neighbour parameters optimized for melting temperature prediction of DNA/RNA hybrids at high and low salt concentrations. Biophysical Chemistry 251:106189. <doi:10.1016/j.bpc.2019.106189>

Banerjee D et al. (2020). Improved nearest-neighbor parameters for the stability of RNA/DNA hybrids under a physiological condition. NAR 48:12042-12054. <doi:10.1093/nar/gkaa572>

Zuber J, Schroeder SJ, Sun H, Turner DH, Mathews DH (2022). Nearest neighbor rules for RNA helix folding thermodynamics: improved end effects. NAR 50:5251-5262. <doi:10.1093/nar/gkac261>

Ghosh S et al. (2020). Nearest-neighbor parameters for predicting DNA duplex stability in diverse molecular crowding conditions. PNAS 117:14194-14201. <doi:10.1073/pnas.1920886117>

Ghosh S et al. (2023). Nearest-neighbor parameters for the prediction of RNA duplex stability in diverse in vitro and cellular-like crowding conditions. NAR 51:4101-4111. <doi:10.1093/nar/gkad020>

Blake RD, Delcourt SG (1996). Thermodynamic effects of formamide on DNA stability. NAR 24:2095-2103. <doi:10.1093/nar/24.11.2095>

See Also

tm_calculate for a single entry point to the nearest-neighbor, GC-content and Wallace methods.

Examples

input_seq <- c("AAAATTTTTTTCCCCCCCCCCCCCCGGGGGGGGGGGGTGTGCGCTGC",
"AAAATTTTTTTCCCCCCCCCCCCCCGGGGGGGGGGGGTGTGCGCTGC")
seqs <- to_genomic_ranges(input_seq)
out <- tm_nn(seqs, Na=50)
out

# A parameter set fitted at a stated sodium concentration. Because Na
# matches the concentration the set was fitted at, salt correction is
# skipped automatically rather than applied on top of it.
out_ls <- tm_nn(seqs, nn_table = "RNA_DNA_NN_Weber_2019_LS", Na = 100)
out_ls$options[["Salt correction applied"]]

# -- A parameter table supplied by the user --------------------------------
# Any of nn_table, tmm_table, imm_table and de_table also accepts a matrix,
# which is how a set the package does not ship, most obviously one covering
# modified bases, is used without waiting for a new release.
#
# Start from a built-in set to get the required keys, then add a stack. The
# key naming follows the built-in convention: top strand, "/", bottom
# strand, so "MG/CG" would be a 5-methylcytosine followed by G, paired with
# CG. The values here are illustrative and are NOT measured parameters.
tbl <- TmCalculator:::get_table("DNA_NN_SantaLucia_2004")
tbl <- rbind(tbl, "MG/CG" = c(-9.1, -24.0))

# Optional attributes. "reference" names the built-in whose key set must be
# present, which matters for RNA and hybrid tables because the default
# reference is a DNA/DNA set. "salt_mM" marks the table as fitted at a
# stated sodium concentration, so that it suppresses the salt correction at
# that concentration exactly as the built-in condition-specific sets do.
attr(tbl, "reference") <- "DNA_NN_SantaLucia_2004"

out_user <- tm_nn(seqs, nn_table = tbl, Na = 50)
out_user$options[["Thermodynamic NN values"]]
# "user-supplied (reference: DNA_NN_SantaLucia_2004)"

# Passing a built-in table back in through this route changes nothing: the
# supplied table is reordered to the reference key order before use, so row
# order carries no information.
identical(tm_nn(seqs, nn_table = "DNA_NN_SantaLucia_2004")$gr$Tm,
          tm_nn(seqs,
                nn_table = TmCalculator:::get_table("DNA_NN_SantaLucia_2004")
                )$gr$Tm)

# A table missing a key is refused rather than tolerated. The compiled core
# resolves each stack by name, so an absent key would contribute zero
# enthalpy and entropy to every sequence containing that step, silently.
try(tm_nn(seqs, nn_table = tbl[setdiff(rownames(tbl), "AA/TT"), ]))
out_ls$options[["Parameter set fitted at [Na+] (mM)"]]

Description

The Wallace rule is often used as rule of thumb for approximate melting temperature calculations for primers with 14 to 20 nt length.

Usage

tm_wallace(gr_seq, ambiguous = FALSE)

Arguments

Value

Returns a list of sequences with updated Tm attributes

Length of validity

The 2 + 4 rule was calibrated on 14 to 20 nt hybridisation probes and carries no length-dependent, salt-dependent or concentration-dependent term. Its error therefore grows without bound as sequences lengthen, and on genome-scale windows it returns a number that is not a melting temperature in any useful sense: a 200 bp window is reported at several hundred degrees Celsius.

Sequences longer than 30 nt raise a warning that names tm_nn as the appropriate alternative. A warning rather than an error, because the rule stays a legitimate rule of thumb and users who knowingly apply it outside its calibrated range should not be blocked; wrap the call in suppressWarnings() in that case.

Author(s)

Junhui Li

References

Thein S L , Lynch J R , Weatherall D J , et al. DIRECT DETECTION OF HAEMOGLOBIN E WITH SYNTHETIC OLIGONUCLEOTIDES[J]. The Lancet, 1986, 327(8472):93.

Examples

input_seq = c('acgtTGCAATGCCGTAWSDBSY','acgtTGCCCCGGCCGCGCCGTAWSDBSY') #for wallace rule
gr_seq <- to_genomic_ranges(input_seq)
out <- tm_wallace(gr_seq, ambiguous = TRUE)
out
out$options

Description

Drop-in replacement for to_genomic_ranges that uses the fast coordinate backend in coor_to_genomic_ranges for genomic coordinate input. Accepts the same input_seq and complement_seq arguments as to_genomic_ranges, plus a list input list(pkg_name = ..., seq = ...) for large tiling jobs.

This function processes a vector of sequences string, a FASTA file, or a character vector with genomic coordinates into a GenomicRanges object, optionally including complementary sequences. sequence names are parsed based on their format: - If names have this pattern "chr:start-end:strand:species[:name]" (e.g., "chr1:1-5:+:seq_1"), parse components into seqnames, ranges, strand, and name. - If names have this pattern "chr:start-end:strand" (e.g., "chr1:1-5:+"), parse components into seqnames, ranges, and strand. - If names have this pattern "chr:start-end" (e.g., "chr1:1-5"), parse components into seqnames and ranges. - If no names are provided, use default values: seqnames = "chr1", start = 1, width = sequence length, strand = "*", name = "1", etc. Complementary sequences are either provided or automatically generated.

Usage

to_genomic_ranges_fast(
  input_seq,
  complement_seq = NULL,
  method = c("vectorized", "preload_chr")
)

to_genomic_ranges(input_seq, complement_seq = NULL)

Arguments

  input_seq       5'-G C A T C G-3'
  complement_seq  3'-C G T A G C-5'
  

Value

A GRanges object. See coor_to_genomic_ranges for metadata columns when coordinate input is used.

A GenomicRanges object with seqnames, ranges, strand, name, sequence, Complement, and Tm as metadata.

Author(s)

Junhui Li

Examples

## Not run: 
gr <- to_genomic_ranges_fast(
  list(
    pkg_name = "BSgenome.Hsapiens.UCSC.hg38",
    seq = c("chr1:1000-1199:+:win1", "chr1:1200-1399:+:win2")
  )
)

## End(Not run)

# Using a character vector with auto-generated complementary sequences
seqs <- c("ATGCG", "GCTAG")
names(seqs) <- c("chr1:1-5:+:seq_1", "chr2:1-5:+")
gr <- to_genomic_ranges(seqs)
gr

# Using a character vector with provided complementary sequences
seqs <- c("ATGCG", "GCTAG")
comp_seqs <- c("TACGC", "CGTA")
gr <- to_genomic_ranges(seqs, comp_seqs)
gr

# Using a FASTA file
gr <- to_genomic_ranges(system.file("extdata", "example1.fasta", package = "TmCalculator"))
## Not run: 
# Using a character vector with genomic coordinates
seqs <- c(
  "chr1:1898000-1898050:+:BSgenome.Hsapiens.UCSC.hg38",
  "chr2:2563000-2563050:-:BSgenome.Hsapiens.UCSC.hg38"
)
gr <- to_genomic_ranges(seqs)
gr

## End(Not run)

Description

This function converts sequence strings to a GenomicRanges object, handling both named and unnamed sequences. It can also process complementary sequences if provided. sequence names can be in the format ">chr2:1-10:+:seq2" which will be parsed into chromosome, position, strand, and name components.

Usage

vec_to_genomic_ranges(input_seq)

Arguments

Value

A GenomicRanges object containing: - GRanges information (seqnames, ranges, strand) - sequence data - Complementary sequences - Names from input or auto-generated

Examples

# Example with named sequences in GRanges format
seqs <- c("ATGCG", "GCTAG")
names(seqs) <- c("chr1:1111-1115:+:seq1", "chr2:1221-1225:+")
gr <- vec_to_genomic_ranges(seqs)

# Example with unnamed sequences
seqs <- c("ATGCG", "GCTAG")
gr <- vec_to_genomic_ranges(seqs)