tm_calculate: Calculate melting temperature using multiple methods

In TmCalculator: Genome-Wide Nucleic Acid Melting Temperature Profiling and Multi-Omics Integration

Calculate melting temperature using multiple methods

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

input_seq	Where the sequence comes from. One of: Sequences, as a character vector in 5' to 3' direction, e.g. c("ATGCG", "GGCCA"). Names, when present, become the seqnames of the result; unnamed sequences are keyed by their position in the vector. An installed BSgenome package, by name, e.g. "BSgenome.Hsapiens.UCSC.hg38". It is named rather than passed as a loaded object because each worker opens the genome for itself, so no sequence crosses between processes. See BSgenome::available.genomes() for what exists, and install the package before calling. A FASTA file, by path; gzipped files are read directly. A GRanges carrying a sequence metadata column. The complement is derived from it when absent. regions selects from any of the four, and means the same thing in each: the identifier before the colon is resolved against whatever names the source itself offers, and falls back to position. A BSgenome or a FASTA file with no regions is taken whole, which for a genome means its standard chromosomes. Also accepted, and unchanged from earlier versions: a character vector of coordinate strings "chr:start-end:strand:species", for example "chr1:100-200:+:BSgenome.Hsapiens.UCSC.hg38", where strand defaults to "+". This form carries its own genome in every element, so it is read directly and ignores regions; to profile coordinates against a genome, prefer passing the genome here and the coordinates as regions.
method	Method(s) to use for Tm calculation. Can be one or more of: - "tm_nn": Nearest Neighbor thermodynamics (default) - "tm_gc": GC content-based method - "tm_wallace": Wallace rule Default: c("tm_nn", "tm_gc", "tm_wallace")
complement_seq	Complementary sequence(s) in 3' to 5' direction. If not provided, the function will automatically generate it from input_seq. This is the template/target sequence that the input sequence will hybridize with. Can be provided as input_seq format besides A NULL value(default)
ambiguous	Logical. If TRUE, ambiguous bases are taken into account when computing the G and C content. The function handles various ambiguous bases (S, W, M, K, R, Y, V, H, D, B) by proportionally distributing their contribution to GC content based on their possible nucleotide compositions. Default: FALSE
shift	Integer value controlling the alignment offset between primer and template sequences. Only applicable for the NN method. Default: 0
nn_table	Thermodynamic nearest-neighbor parameters for different nucleic acid hybridizations. Only applicable for the NN method. Sets whose name encodes a sodium concentration were fitted at that condition and are not salt-corrected again. See tm_nn for the full list and guidance on choosing between them. Default: "DNA_NN_SantaLucia_2004" Alternatively, supply a matrix or data.frame of parameters directly. This is the route for parameter sets the package does not ship, in particular sets covering modified bases such as 5-methylcytosine. Requirements: numeric, with columns 1 and 2 read as delta H (kcal/mol) and delta S (cal/mol/K); further columns are ignored; row names giving the parameter keys, e.g. "AA/TT", "init", "init_A/T", "sym"; every key of a built-in reference set must be present. The reference is named by attr(x, "reference"), or defaults to the first built-in listed for the argument, which is a DNA/DNA set; RNA and hybrid tables should therefore set the attribute. Extra keys beyond the reference are kept, which is how a modified-base set adds stacks rather than replacing them. The supplied table is reordered to the reference key order before use, so that two tables differing only in row order give identical results. A missing key would otherwise contribute zero to the calculation instead of raising an error, which is why the full key set is required. Keys that disagree with their reverse complement produce a warning: expected for modified bases, a transposition error otherwise. Two optional attributes are honoured. attr(x, "salt_mM") marks a set as fitted at a stated sodium concentration, which suppresses the salt correction at that concentration exactly as for the built-in sets fitted this way; without it the table is treated as a reference-condition set and salt_method is applied. attr(x, "end_table") supplies a companion penultimate-pair end-effect table.
tmm_table	Thermodynamic parameters for terminal mismatches. Only applicable for the NN method. Default: "DNA_TMM_Bommarito_2000"
imm_table	Thermodynamic parameters for internal mismatches. Only applicable for the NN method. Default: "DNA_IMM_Peyret_1999"
de_table	Thermodynamic parameters for dangling ends. Only applicable for the NN method. Default: "DNA_DE_Bommarito_2000"
dnac_high	Concentration of the higher concentrated strand in nM. Only applicable for the NN method. Default: 25
dnac_low	Concentration of the lower concentrated strand in nM. Only applicable for the NN method. Default: 25
self_comp	Logical value indicating if the sequence is self-complementary. Only applicable for the NN method. Default: FALSE
variant	Empirical constants coefficient for GC method. Only applicable for the GC method. Default: "Primer3Plus"
userset	A vector of four coefficient values for GC method. Only applicable for the GC method. Usersets override value sets. Default: NULL
Na	Millimolar concentration of sodium ions. Default: 50
K	Millimolar concentration of potassium ions. Default: 0
Tris	Millimolar concentration of Tris buffer. Default: 0
Mg	Millimolar concentration of magnesium ions. Default: 0
dNTPs	Millimolar concentration of deoxynucleotide triphosphates. Default: 0
salt_method	Salt correction method for Tm. Default: "Schildkraut2010" Available options: - "none": Disables salt correction. Also selected automatically when the chosen nn_table was fitted at the requested Na. - "Schildkraut2010": Schildkraut & Lifson (1965); historical identifier - "Wetmur1991": Classic salt correction method - "SantaLucia1996": DNA-specific salt correction - "SantaLucia1998-1": Improved DNA salt correction, applied to Tm - "SantaLucia1998-2": the same correction applied to the entropy of the nearest-neighbor model rather than to Tm - "Owczarzy2004": Comprehensive salt correction - "Owczarzy2008": Updated comprehensive salt correction Default: "Schildkraut2010" With method = "tm_gc" the salt term is part of the published formula selected by variant, so naming a different one is ignored there, with a warning, unless userset is supplied. "none" (or NA) is not a substitution but a request to drop the correction, and is honoured on either path without a warning. The two Owczarzy corrections are not available to "tm_gc" at all: they apply to the reciprocal of the melting temperature in kelvin, referenced to the same duplex in 1 M Na+, and carry a duplex-length term the GC-content formulas already have.
DMSO	Percent DMSO concentration in the reaction mixture. Default: 0
formamide_unit	Formamide concentration as 'list(value, unit)'. Default: list(value = 0, unit = "percent") - value: Numeric value of formamide concentration - unit: Either "percent" or "molar"
dmso_factor	Coefficient of Tm decreases per percent DMSO. Default: 0.75 Other published values are 0.5, 0.6 and 0.675.
formamide_factor	Tm decrease per percent formamide. Default: 0.65 Several papers report factors between 0.6 and 0.72.
mismatch	Logical. If TRUE, every '.' in the sequence is counted as a mismatch. Only applicable for the GC method. Default: TRUE
regions	What to take from input_seq. NULL, the default, means all of it, except for a BSgenome, where it means GenomeInfoDb::standardChromosomes() of that genome, which for GRCh38 includes chrM. Otherwise: names or numbers (c("chr1", "chr2"), 1:2), coordinate strings "name:start-end" with commas and scientific notation accepted ("chr1:5,000,000-6e6"), a mixture of the two, or a GRanges. The identifier before the colon is resolved against whatever names the source itself offers, and falls back to position. So "chr1" is a chromosome in a BSgenome, a record in a FASTA file and a seqname in a GRanges, and "1:1-200" is the first 200 bases of the first sequence in an unnamed character vector. On a BSgenome the chr prefix is added or removed as the genome requires, since that is a convention rather than information; elsewhere names are matched exactly.
window	Window width in base pairs. NULL, the default, gives one melting temperature per region, which is what short records such as probes, primers and oligonucleotides call for; a region longer than 1 Mb with window = NULL is an error rather than one meaningless Tm. Regions shorter than window are returned whole.
slide	Step between window starts, defaulting to window for a non-overlapping tiling. Ignored when window is NULL.
unit	How regions become tasks: "segment" cuts them into pieces of about segment_size bp, "region" makes one task per region. Segments are faster and need less memory per worker, because no worker then holds a whole large chromosome.
segment_size	Task size in base pairs when unit = "segment", rounded down to a multiple of slide so that the window grid is the one an unsegmented run would produce. Default 50 Mb.
BPPARAM	A BiocParallelParam from BiocParallel, e.g. SnowParam(workers = 5), to spread the tasks over processes. NULL, the default, runs them here, with no dependency on BiocParallel. Parallelism divides the work by region, never the sequences of one region: each task opens the source itself, so only a name and a coordinate pair cross between processes.
keep_sequence	Keep the sequence and complement columns. The default keeps them for sequences the caller supplied and drops them for a genome or a file, where they run to roughly 500 MB per large chromosome.
tmpdir	Directory for the temporary FASTA file written when sequences are supplied directly and there is tiling or parallelism to do. Worth setting on a cluster, where tempdir() is often a small partition.
verbose	Report the task and window counts.

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:

gr	The input GRanges with metadata columns Tm and GC (melting temperature in ^{\circ}C and GC percent).
options	Calculation parameters and method information. For the nearest-neighbor method this includes Salt correction applied (logical) and Parameter set fitted at [Na+] (mM), which record whether a salt correction was actually performed.

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)

TmCalculator documentation built on Oct. 5, 2026, 5:08 p.m.