inst/doc/segmentSeq.R

In segmentSeq: Methods for identifying small RNA loci from high-throughput sequencing data

### R code from vignette source 'segmentSeq.Rnw'

###################################################
### code chunk number 1: <style-Sweave
###################################################
BiocStyle::latex()


###################################################
### code chunk number 2: segmentSeq.Rnw:30-31
###################################################
  library(segmentSeq)


###################################################
### code chunk number 3: segmentSeq.Rnw:36-43
###################################################
if(require("parallel")) 
{
    numCores <- min(8, detectCores())
    cl <- makeCluster(numCores)
} else {
    cl <- NULL
}


###################################################
### code chunk number 4: segmentSeq.Rnw:46-47 (eval = FALSE)
###################################################
## cl <- NULL


###################################################
### code chunk number 5: segmentSeq.Rnw:50-51
###################################################
set.seed(1)


###################################################
### code chunk number 6: segmentSeq.Rnw:59-68
###################################################
datadir <- system.file("extdata", package = "segmentSeq")
libfiles <- c("SL9.txt", "SL10.txt", "SL26.txt", "SL32.txt")
libnames <- c("SL9", "SL10", "SL26", "SL32")
replicates <- c("AGO6", "AGO6", "AGO4", "AGO4")

aD <- readGeneric(files = libfiles, dir = datadir,
                  replicates = replicates, libnames = libnames,                   
                  polyLength = 10, header = TRUE, gap = 200)
aD


###################################################
### code chunk number 7: segmentSeq.Rnw:73-75
###################################################
sD <- processAD(aD, cl = cl)
sD


###################################################
### code chunk number 8: segmentSeq.Rnw:79-80
###################################################
if(nrow(sD) != 1452) stop("sD object is the wrong size (should have 1452 rows). Failure.")


###################################################
### code chunk number 9: segmentSeq.Rnw:89-90
###################################################
hS <- heuristicSeg(sD = sD, aD = aD, RKPM = 1000, largeness = 1e8, getLikes = TRUE, cl = cl)


###################################################
### code chunk number 10: segmentSeq.Rnw:93-95
###################################################
if(nrow(hS) != 507) stop("hS object is the wrong size (should have 507 rows). Failure.")
if(any(abs(colSums(exp(hS@locLikelihoods)) - c(88, 208)) > 2)) stop("hS object contains wrong number of loci. Likely failure.")


###################################################
### code chunk number 11: segmentSeq.Rnw:107-109
###################################################
cS <- classifySeg(sD = sD, aD = aD, cD = hS, cl = cl)      
cS


###################################################
### code chunk number 12: segmentSeq.Rnw:112-114
###################################################
if(abs(nrow(cS) - 64) > 2) stop("cS object is the wrong size (should have ~142 rows). Likely failure.")
if(any(abs(colSums(exp(cS@locLikelihoods)) - c(29,36)) > 2)) stop("cS object contains wrong number of loci. Likely failure.")


###################################################
### code chunk number 13: segPlot
###################################################
par(mfrow = c(2,1), mar = c(2,6,2,2))
plotGenome(aD, hS, chr = ">Chr1", limits = c(1, 1e5),
           showNumber = FALSE, cap = 50)
plotGenome(aD, cS, chr = ">Chr1", limits = c(1, 1e5),
           showNumber = FALSE, cap = 50)


###################################################
### code chunk number 14: figSeg
###################################################
par(mfrow = c(2,1), mar = c(2,6,2,2))
plotGenome(aD, hS, chr = ">Chr1", limits = c(1, 1e5),
           showNumber = FALSE, cap = 50)
plotGenome(aD, cS, chr = ">Chr1", limits = c(1, 1e5),
           showNumber = FALSE, cap = 50)


###################################################
### code chunk number 15: segmentSeq.Rnw:146-148
###################################################
loci <- selectLoci(cS, FDR = 0.05)
loci


###################################################
### code chunk number 16: segmentSeq.Rnw:154-155
###################################################
groups(cS) <- list(NDE = c(1,1,1,1), DE = c("AGO6", "AGO6", "AGO4", "AGO4"))


###################################################
### code chunk number 17: segmentSeq.Rnw:159-160
###################################################
cS <- getPriors(cS, cl = cl)


###################################################
### code chunk number 18: segmentSeq.Rnw:164-165
###################################################
cS <- getLikelihoods(cS, nullData = TRUE, cl = cl)


###################################################
### code chunk number 19: segmentSeq.Rnw:170-171
###################################################
topCounts(cS, NULL, number = 3)


###################################################
### code chunk number 20: segmentSeq.Rnw:176-177
###################################################
topCounts(cS, "NDE", number = 3)


###################################################
### code chunk number 21: segmentSeq.Rnw:182-183
###################################################
topCounts(cS, "DE", number = 3)


###################################################
### code chunk number 22: <stopCluster
###################################################
if(!is.null(cl))
    stopCluster(cl)


###################################################
### code chunk number 23: segmentSeq.Rnw:195-196
###################################################
sessionInfo()