## ----style-knitr, eval=TRUE, echo=FALSE, results="asis"-------------------- BiocStyle::latex() ## -------------------------------------------------------------------------- suppressPackageStartupMessages({ library(BUMHMM) library(Biostrings) library(SummarizedExperiment) }) se ## -------------------------------------------------------------------------- controls <- se[, se$replicate == "control"] head(assay(controls, 'coverage')) ## -------------------------------------------------------------------------- rowData(controls)[1:4,] ## -------------------------------------------------------------------------- colData(controls) ## -------------------------------------------------------------------------- pos <- 300 assay(controls, 'coverage')[pos, 1] assay(controls, 'dropoff_count')[pos, 1] assay(controls, 'dropoff_rate')[pos, 1] ## -------------------------------------------------------------------------- treatments <- se[, se$replicate == "treatment"] assay(treatments, 'coverage')[pos, 1] assay(treatments, 'dropoff_count')[pos, 1] assay(treatments, 'dropoff_rate')[pos, 1] ## -------------------------------------------------------------------------- Nc <- Nt <- 3 t <- 1 nuclSelection <- selectNuclPos(se, Nc, Nt, t) List(nuclSelection) ## -------------------------------------------------------------------------- t(combn(Nc, 2)) ## -------------------------------------------------------------------------- length(nuclSelection$analysedC[[1]]) length(nuclSelection$analysedCT[[1]]) ## -------------------------------------------------------------------------- ## Medians of original drop-off rates in each replicate apply(assay(se, 'dropoff_rate'), 2, median) ## Scale drop-off rates assay(se, "dropoff_rate") <- scaleDOR(se, nuclSelection, Nc, Nt) ## Medians of scaled drop-off rates in each replicate apply(assay(se, 'dropoff_rate'), 2, median) ## -------------------------------------------------------------------------- stretches <- computeStretches(se, t) ## -------------------------------------------------------------------------- head(stretches) assay(se, 'dropoff_count')[1748,] ## -------------------------------------------------------------------------- varStab <- stabiliseVariance(se, nuclSelection, Nc, Nt) LDR_C <- varStab$LDR_C LDR_CT <- varStab$LDR_CT hist(LDR_C, breaks = 30, main = 'Null distribution of LDRs') ## -------------------------------------------------------------------------- nuclNum <- 3 patterns <- nuclPerm(nuclNum) patterns ## -------------------------------------------------------------------------- ## Extract the DNA sequence sequence <- subject(rowData(se)$nucl) sequence nuclPosition <- findPatternPos(patterns, sequence, '+') patterns[[1]] head(nuclPosition[[1]]) ## -------------------------------------------------------------------------- nuclPosition <- list() nuclPosition[[1]] <- 1:nchar(sequence) ## Start of the stretch nuclPosition[[1]][1] ## End of the stretch nuclPosition[[1]][length(nuclPosition[[1]])] ## -------------------------------------------------------------------------- posteriors <- computeProbs(LDR_C, LDR_CT, Nc, Nt, '+', nuclPosition, nuclSelection$analysedC, nuclSelection$analysedCT, stretches) ## -------------------------------------------------------------------------- head(posteriors) ## -------------------------------------------------------------------------- shifted_posteriors <- matrix(, nrow=dim(posteriors)[1], ncol=1) shifted_posteriors[1:(length(shifted_posteriors) - 1)] <- posteriors[2:dim(posteriors)[1], 2] ## -------------------------------------------------------------------------- plot(shifted_posteriors, xlab = 'Nucleotide position', ylab = 'Probability of modification', main = 'BUMHMM output for 18S DMS data set') ## ----eval=FALSE------------------------------------------------------------ # ## Call the function with the additonal tolerance parameter # posteriors <- computeProbs(LDR_C, LDR_CT, Nc, Nt, '+', nuclPosition, # nuclSelection$analysedC, nuclSelection$analysedCT, # stretches, 0.001) ## -------------------------------------------------------------------------- sessionInfo()