transit: TRANSIT algorithm for detecting jumps

In stepR: Multiscale Change-Point Inference

TRANSIT algorithm for detecting jumps

Description

Reimplementation of VanDongen's algorithm for detecting jumps in ion channel recordings.

Deprecation warning: This function is mainly used for patchlamp recordings and may be transferred to a specialised package.

Usage

transit(y, x = 1:length(y), x0 = 2 * x[1] - x[2], sigma.amp = NA, sigma.slope = NA,
  amp.thresh = 3, slope.thresh = 2, rel.amp.n = 3, rel.amp.thresh = 4,
  family = c("gauss", "gaussKern"), param = NULL, refit = FALSE)

Arguments

y	a numeric vector containing the serial data
sigma.amp	amplitude (i.e. raw data within block) standard deviation; estimated using sdrobnorm if omitted
sigma.slope	slope (i.e. central difference within block) standard deviation; estimated using sdrobnorm if omitted
amp.thresh	amplitude threshold
slope.thresh	slope threshold
rel.amp.n	relative amplitude threshold will be used for blocks with no more datapoints than this
rel.amp.thresh	relative amplitude threshold
x	a numeric vector of the same length as y containing the corresponding sample points
x0	a single numeric giving the last unobserved sample point directly before sampling started
family, param	specifies distribution of data, see family
refit	should the values for family = "gaussKern" be obtained by fitting in the end (otherwise they are meaningless)

Value

Returns an object of class stepfit which encodes the jumps and corresponding mean values.

Note

Only central, no forward differences have been used in this implementation. Moreover, the standard deviations will be estimated by sdrobnorm if omitted (respecting the filter's effect if applicable).

References

VanDongen, A. M. J. (1996) A new algorithm for idealizing single ion channel data containing multiple unknown conductance levels. Biophysical Journal 70(3), 1303–1315.

See Also

stepfit, sdrobnorm, jsmurf, stepbound, steppath

Examples

# estimating step-functions with Gaussian white noise added
# simulate a Gaussian hidden Markov model of length 1000 with 2 states
# with identical transition rates 0.01, and signal-to-noise ratio 2
sim <- contMC(1e3, 0:1, matrix(c(0, 0.01, 0.01, 0), 2), param=1/2)
plot(sim$data, cex = 0.1)
lines(sim$cont, col="red")
# maximum-likelihood estimation under multiresolution constraints
fit.MRC <- smuceR(sim$data$y, sim$data$x)
lines(fit.MRC, col="blue")
# choose number of jumps using BIC
path <- steppath(sim$data$y, sim$data$x, max.blocks=1e2)
fit.BIC <- path[[stepsel.BIC(path)]]
lines(fit.BIC, col="green3", lty = 2)

# estimate after filtering
# simulate filtered ion channel recording with two states
set.seed(9)
# sampling rate 10 kHz
sampling <- 1e4
# tenfold oversampling
over <- 10
# 1 kHz 4-pole Bessel-filter, adjusted for oversampling
cutoff <- 1e3
df.over <- dfilter("bessel", list(pole=4, cutoff=cutoff / sampling / over))
# two states, leaving state 1 at 10 Hz, state 2 at 20 Hz
rates <- rbind(c(0, 10), c(20, 0))
# simulate 0.5 s, level 0 corresponds to state 1, level 1 to state 2
# noise level is 0.3 after filtering
Sim <- contMC(0.5 * sampling, 0:1, rates, sampling=sampling, family="gaussKern",
  param = list(df=df.over, over=over, sd=0.3))
plot(Sim$data, pch = ".")
lines(Sim$discr, col = "red")
# fit under multiresolution constraints using filter corresponding to sample rate
df <- dfilter("bessel", list(pole=4, cutoff=cutoff / sampling))
Fit.MRC <- jsmurf(Sim$data$y, Sim$data$x, param=df, r=1e2)
lines(Fit.MRC, col = "blue")
# fit using TRANSIT
Fit.trans <- transit(Sim$data$y, Sim$data$x)
lines(Fit.trans, col = "green3", lty=2)

stepR documentation built on Nov. 14, 2023, 1:09 a.m.