Estimate the nucleosome repeat length (NRL) from modified-base distances.

Estimate the nucleosome repeat length (NRL) from the frequencies of same-read modified base distances, e.g. generated by calcModbaseSpacing. The NRL is obtained from the slope of a linear fit to the modes in the distance distribution.

Usage

estimateNRL(
  x,
  minDist = 140L,
  usePeaks = seq_len(5),
  span1 = 100/length(x),
  span2 = 1500/length(x)
)

Arguments

x

numeric vector giving the counts of distances (typically the calculated with calcModbaseSpacing.

minDist

integer(1) specifying the minimal distance to be used for NRL estimation. The default value (140) ignores any distance too short to span at least a single nucleosome.

usePeaks

integer vector selecting the modes (peaks) in the phasogram used in NRL estimation.

span1

numeric(1) giving the smoothing parameter for de-trending loess fit (high pass filter).

span2

numeric(1) giving the smoothing parameter for de-noising loess fit (low pass filter).

Value

A list with elements:

nrl

the estimated nucleosome repeat length

nrl.CI95

the 95-percent confidence interval

xs

smoothed (de-trended) phasogram

loessfit

the de-noising fit to the de-trended phasogram

lmfit

the linear fit to the phasogram peaks

minDist

minimal distance included in the fit

span1

smoothing parameter for de-trending loess fit

span2

smoothing parameter for de-noising loess fit

usePeaks

the peaks used in the fit

See also

calcModbaseSpacing to calculate the distances from base modification data, plotModbaseSpacing to visualize annotated distance frequencies between modified bases.

Author

Michael Stadler

Examples

# read base modifications
modbamfiles <- system.file("extdata",
                           c("6mA_1_10reads.bam", "6mA_2_10reads.bam"),
                           package = "footprintR")
se <- readModBam(modbamfiles, "chr1:6940000-6955000", "a")

# get distances for each sample
moddist <- calcModbaseSpacing(se)

# analyze NRL for each sample
print(estimateNRL(moddist$s1)[1:2])
#> $nrl
#> [1] 182.6
#> 
#> $nrl.CI95
#>    2.5 %   97.5 % 
#> 179.5475 185.6525 
#> 
print(estimateNRL(moddist$s2)[1:2])
#> $nrl
#> [1] 113.4
#> 
#> $nrl.CI95
#>     2.5 %    97.5 % 
#>  77.74145 149.05855 
#> 

# combine samples
moddistComb <- Reduce("+", moddist)
print(estimateNRL(moddistComb)[1:2])
#> $nrl
#> [1] 184.3
#> 
#> $nrl.CI95
#>    2.5 %   97.5 % 
#> 180.7993 187.8007 
#>