Skip to contents

Quantify windows by estimating nucleosome repeat lengths (NRLs)

Source: R/scanGenome.R
estimateNRLwindows.Rd

Estimate nucleosome repeath length in windows using calcModbaseSpacing and estimateNRL.

Usage

estimateNRLwindows(
  se,
  gr,
  dmax = 1000L,
  assayName = "mod_prob",
  minModProb = 0.5,
  minDist = 140L,
  usePeaks = seq_len(5),
  minperiod1 = 30,
  minperiod2 = 450,
  BPPARAM = BiocParallel::MulticoreParam(4L)
)

Arguments

se

A RangedSummarizedExperiment object, as generated by readModBam(..., level = "quickread") containing a read-level assay called assayName.

gr

A GRanges object defining the windows to quantify.

dmax

Numeric scalar specifying the maximal distance between modified bases on the same read to count.

assayName

A character scalar specifying the assay of se containing the read-level modification probabilities.

minModProb

Numeric scalar giving the minimal modification probability for a modified base.

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.

minperiod1, minperiod2

numeric scalars giving the smoothing parameter for de-trending the signal (minimal periods).

BPPARAM

A BiocParallelParam object that controls the number of parallel CPU threads to use for some of the steps in estimateNRLwindows(). The default value is (MulticoreParam(4L)).

Value

A RangedSummarizedExperiment object with one row per window defined by gr and the NRL estimates. This object is suitable for transformation by getRangesWithAssayValues.

Author

Michael Stadler

Examples

library(GenomicRanges)
library(SummarizedExperiment)
modbamfiles <- system.file("extdata",
                           c("6mA_1_10reads.bam", "6mA_2_10reads.bam"),
                           package = "footprintR")
rng <- GRanges("chr1", IRanges(6930000, 6940000))
windowSize <- 2000
windowStep <- 1000
s <- seq(start(rng), end(rng) - windowSize + 1, by = windowStep)
windowgr <- GRanges(seqnames = seqnames(rng),
                    ranges = IRanges(start = s, width = windowSize))
se <- readModBam(bamfiles = modbamfiles, regions = rng, level = "quickread",
                 modbase = "a", trim = TRUE,
                 BPPARAM = BiocParallel::SerialParam())
seNRL <- estimateNRLwindows(se = se, gr = windowgr,
                            BPPARAM = BiocParallel::SerialParam())
assay(seNRL, "NRL")
#>      s1    s2
#> 1 184.4 188.2
#> 2 184.8 191.7
#> 3 175.5 192.4
#> 4 180.8 181.8
#> 5 180.1 181.2
#> 6 187.6 194.8
#> 7 107.2 188.3
#> 8 146.3 187.2
#> 9 211.4 154.4
assay(seNRL, "NRL.CI95low")
#>          s1       s2
#> 1 176.66550 184.0587
#> 2 172.39798 182.5887
#> 3 154.33015 187.7084
#> 4 176.21021 166.8775
#> 5 177.09769 167.5291
#> 6 180.22295 189.7455
#> 7  71.47524 165.1526
#> 8 100.37600 177.2306
#> 9 178.03848 126.9546
assay(seNRL, "NRL.CI95high")
#>         s1       s2
#> 1 192.1345 192.3413
#> 2 197.2020 200.8113
#> 3 196.6698 197.0916
#> 4 185.3898 196.7225
#> 5 183.1023 194.8709
#> 6 194.9771 199.8545
#> 7 142.9248 211.4474
#> 8 192.2240 197.1694
#> 9 244.7615 181.8454