transit {stepR}

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)