Score DNA sequences with a PWM over a region of interest

Scores full DNA sequences using a Position Weight Matrix (PWM) over a specified region of interest (ROI). The ROI is defined by start_pos and end_pos (1-based, inclusive), with optional extension controlled by extend. All reported positions are on the full input sequence.

Usage

gseq.pwm(
  seqs,
  pssm,
  mode = c("lse", "max", "pos", "count"),
  bidirect = TRUE,
  strand = 1L,
  score.thresh = NULL,
  start_pos = NULL,
  end_pos = NULL,
  extend = FALSE,
  spat.factor = NULL,
  spat.bin = 1L,
  spat.min = NULL,
  spat.max = NULL,
  return_strand = FALSE,
  skip_gaps = TRUE,
  gap_chars = c("-", "."),
  neutral_chars = character(0),
  neutral_chars_policy = c("average", "log_quarter", "na"),
  prior = 0.01
)

Arguments

seqs: character vector of DNA sequences (A/C/G/T/N; case-insensitive)
pssm: numeric matrix or data frame with columns named A, C, G, T (additional columns are allowed and will be ignored). Each row holds the probabilities or counts of the four bases at one motif position; log-odds matrices, NA values and rows summing to 0 are rejected.
mode: character; one of "lse", "max", "pos", or "count"
bidirect: logical; if TRUE, scans both strands (default: TRUE)
strand: integer; 1 for the forward strand, -1 for the reverse (default: 1). Used only when bidirect = FALSE, and required there - 0 is refused, since exactly one strand is read and 0 cannot name it. Under bidirect = TRUE both strands are read and this is ignored.
score.thresh: single number; windows scoring at or above this value are counted. Required when mode="count" and ignored otherwise. PWM scores are log-likelihoods, so the usable range depends on the PSSM, the prior and any spatial weights - calibrate with mode="max" before choosing one.
start_pos: integer or NULL; 1-based inclusive start of ROI (default: 1)
end_pos: integer or NULL; 1-based inclusive end of ROI (default: sequence length)
extend: logical or integer; if TRUE, a motif anchored in the ROI may run up to w - 1 bases past its end (w being the motif length), so boundary-anchored motifs keep their full context; anchors must still start inside the ROI. If FALSE (default), only windows fully contained in the ROI are scored. An integer is the number of bases a window may run past end_pos: TRUE is w - 1 and FALSE is 0. An integer above w - 1 also scores windows that start after end_pos.
spat.factor: numeric vector; spatial weighting factors (optional)
spat.bin: integer; bin size for spatial weighting
spat.min: numeric; 1-based position, counted from start_pos, where the scanning window starts; spatial bins are counted from it (optional)
spat.max: numeric; 1-based position, counted from start_pos, where the scanning window ends: only windows ending at or before it are scored (optional)
return_strand: logical; if TRUE and mode="pos", returns data.frame with pos and strand columns
skip_gaps: logical; if TRUE, treat gap characters as holes and skip them while scanning. Windows are w consecutive non-gap bases (default: TRUE)
gap_chars: character vector; which characters count as gaps (default: c("-", "."))
neutral_chars: character vector; characters to score as unknown bases, by neutral_chars_policy, instead of leaving their windows out. Default: none, so a window holding N or any other non-ACGT character is not scored, as in the PWM virtual tracks. Pass e.g. c("N", "n", "*") to score them.
neutral_chars_policy: character string; how to score the characters listed in neutral_chars. One of "average" (the column's mean log-probability), "log_quarter" (always log(1/4)), or "na" (return NA when one is in the scanning window). Has no effect when neutral_chars is empty.
prior: numeric; pseudocount added to frequencies (default: 0.01). Set to 0 for no pseudocounts.

Value

Numeric vector (for "lse"/"max"/"count" modes), integer vector (for "pos" mode), or data.frame with pos and strand columns (for "pos" mode with return_strand=TRUE). A sequence whose windows are all unscorable (for example all N) gives 0 for "count", as pwm.count does. A sequence or ROI with no window at all (shorter than the motif), or an NA sequence, gives NA in every mode, "count" included, as the PWM virtual tracks give for an interval without an anchor.

Details

This function scores DNA sequences directly without requiring a genomics database. For detailed documentation on PWM scoring modes, parameters, and spatial weighting, see gvtrack.create (functions "pwm", "pwm.max", "pwm.max.pos", "pwm.count").

The ROI (region of interest) is defined by start_pos and end_pos. As in the PWM virtual tracks, a window on either strand is anchored at its first (leftmost) base, which must not lie before start_pos. extend sets how many bases past end_pos a window may run (w - 1 with TRUE, so every anchor in the ROI is scored), never past the end of the sequence.

When skip_gaps=TRUE, characters specified in gap_chars are treated as gaps. Windows are defined as w consecutive non-gap bases, and extend counts non-gap bases. All positions (pos) are reported as 1-based indices on the original full sequence (including gaps). start_pos and end_pos are interpreted as physical coordinates on the full sequence.

By default a window holding N or any other character that is not A, C, G or T (in either case) is not scored on either strand: it is left out, exactly as the pwm, pwm.max, pwm.max.pos and pwm.count virtual tracks leave it out, and a sequence with no scorable window gives -Inf for "lse" and "max", NA for "pos" and 0 for "count". Characters listed in neutral_chars are instead scored as unknown bases, the same on both strands, as neutral_chars_policy says.

Under bidirect = TRUE the forward and reverse-complement matches are two readings of the same window, and every mode but "pos" combines them at each window by log-sum-exp, after any spatial weights are applied: "lse" then log-sum-exps those unions across the windows, "max" takes the largest of them, and "count" thresholds each one once, so a window is at most a single hit however well its reverse complement matches it. That is the same union the pwm, pwm.max and pwm.count virtual tracks take, and over the same windows, with the default neutral_chars, on a sequence read from the genome (which holds no gap characters), gseq.pwm() agrees with them to the precision the two share - the virtual tracks accumulate in single precision and gseq.pwm() in double, so expect agreement to about 1e-5 relative, not to the last bit. mode = "pos" instead takes the better of the two strands, because it has to name one - what the pwm.max.pos virtual track does - so "max" and "pos" can select different windows, and the value from "max" is not necessarily the score at the position from "pos". The asymmetry is the pwm family's own (pwm.max combines the strands, pwm.max.pos does not) and is mirrored here rather than fixed on one side. With bidirect = FALSE only the strand named by strand is read and there is nothing to combine.

Examples