inst/doc/NOISeq.R

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

###################################################
### code chunk number 1: options
###################################################
options(digits=3, width=95)


###################################################
### code chunk number 2: data
###################################################
library(NOISeq)
data(Marioni)


###################################################
### code chunk number 3: NOISeq.Rnw:88-89
###################################################
head(mycounts)


###################################################
### code chunk number 4: factors
###################################################
myfactors = data.frame(Tissue=c("Kidney","Liver","Kidney","Liver","Liver","Kidney","Liver",
                                "Kidney","Liver","Kidney"),
                       TissueRun = c("Kidney_1","Liver_1","Kidney_1","Liver_1","Liver_1",
                                     "Kidney_1","Liver_1","Kidney_2","Liver_2","Kidney_2"),
                       Run = c(rep("R1", 7), rep("R2", 3)))
myfactors


###################################################
### code chunk number 5: NOISeq.Rnw:119-123
###################################################
head(mylength)
head(mygc)
head(mybiotypes)
head(mychroms)


###################################################
### code chunk number 6: readData
###################################################
mydata <- readData(data=mycounts,length=mylength, gc=mygc, biotype=mybiotypes,
                   chromosome=mychroms, factors=myfactors)
mydata


###################################################
### code chunk number 7: NOISeq.Rnw:152-156
###################################################
str(mydata)
head(assayData(mydata)$exprs)
head(pData(mydata))
head(featureData(mydata)@data)


###################################################
### code chunk number 8: readData2
###################################################
mydata <- readData(data=mycounts,chromosome=mychroms, factors=myfactors)


###################################################
### code chunk number 9: readData3
###################################################
mydata <- addData(mydata, length=mylength, biotype=mybiotypes, gc = mygc)


###################################################
### code chunk number 10: dat
###################################################
myexplodata <- dat(mydata, type = "biodetection")
explo.plot(myexplodata, plottype = "persample")


###################################################
### code chunk number 11: nicedata
###################################################
mynicedata <- dat2save(myexplodata)


###################################################
### code chunk number 12: fig_biodetection
###################################################
mybiodetection <- dat(mydata, k = 0, type = "biodetection", factor = NULL)
par(mfrow = c(1,2))  # we need this instruction because two plots (one per sample) will be generated
explo.plot(mybiodetection, samples=c(1,2), plottype = "persample")


###################################################
### code chunk number 13: fig_biodetection2
###################################################
par(mfrow = c(1,2))  # we need this instruction because two plots (one per sample) will be generated
explo.plot(mybiodetection, samples=c(1,2), toplot = "protein_coding", plottype = "comparison")


###################################################
### code chunk number 14: fig_boxplot1
###################################################
mycountsbio = dat(mydata, factor = NULL, type = "countsbio")
explo.plot(mycountsbio, toplot = 1, samples = 1, plottype = "boxplot")


###################################################
### code chunk number 15: fig_sat1
###################################################
mysaturation = dat(mydata, k = 0, ndepth = 7, type = "saturation")
explo.plot(mysaturation, toplot = 1, samples = 1:2, yleftlim = NULL, yrightlim = NULL)


###################################################
### code chunk number 16: fig_sat2
###################################################
explo.plot(mysaturation, toplot = "protein_coding", samples = 1:4)


###################################################
### code chunk number 17: fig_boxplot2
###################################################
explo.plot(mycountsbio, toplot = "protein_coding", samples = NULL, plottype = "boxplot")


###################################################
### code chunk number 18: fig_boxplot3
###################################################
explo.plot(mycountsbio, toplot = 1, samples = NULL, plottype = "barplot")


###################################################
### code chunk number 19: fig_length
###################################################
mylengthbias = dat(mydata, factor = "Tissue", type = "lengthbias")
explo.plot(mylengthbias, samples = NULL, toplot = "global")


###################################################
### code chunk number 20: showmodels
###################################################
show(mylengthbias)


###################################################
### code chunk number 21: fig_GC
###################################################
myGCbias = dat(mydata, factor = "Tissue", type = "GCbias")
explo.plot(myGCbias, samples = NULL, toplot = "global")


###################################################
### code chunk number 22: fig_countdistr
###################################################
mycd = dat(mydata, type = "cd", norm = FALSE, refColumn = 1)
explo.plot(mycd)


###################################################
### code chunk number 23: randomBatchEffect
###################################################
set.seed(123)
mycounts2 = mycounts
mycounts2[,1:4] = mycounts2[,1:4] + runif(nrow(mycounts2)*4, 3, 5)
myfactors = data.frame(myfactors, "batch" = c(rep(1,4), rep(2,6)))
mydata2 = readData(mycounts2, factors = myfactors)


###################################################
### code chunk number 24: fig_PCA
###################################################
myPCA = dat(mydata2, type = "PCA")
par(mfrow = c(1,2))
explo.plot(myPCA, factor = "Tissue")
explo.plot(myPCA, factor = "batch")


###################################################
### code chunk number 25: QCreportExample
###################################################
QCreport(mydata, samples = NULL, factor = "Tissue", norm = FALSE)


###################################################
### code chunk number 26: normalization
###################################################
myRPKM = rpkm(assayData(mydata)$exprs, long = mylength, k = 0, lc = 1)
myUQUA = uqua(assayData(mydata)$exprs, long = mylength, lc = 0.5, k = 0)
myTMM = tmm(assayData(mydata)$exprs, long = 1000, lc = 0)
head(myRPKM[,1:4])


###################################################
### code chunk number 27: filtering
###################################################
myfilt = filtered.data(mycounts, factor = myfactors$Tissue, norm = FALSE, depth = NULL, method = 1, cv.cutoff = 100, cpm = 1, p.adj = "fdr")


###################################################
### code chunk number 28: fig_knownBatch
###################################################
mydata2corr1 = ARSyNseq(mydata2, factor = "batch", batch = TRUE, norm = "rpkm",  logtransf = FALSE)
myPCA = dat(mydata2corr1, type = "PCA")
par(mfrow = c(1,2))
explo.plot(myPCA, factor = "Tissue")
explo.plot(myPCA, factor = "batch")


###################################################
### code chunk number 29: fig_unknownBatch
###################################################
mydata2corr2 = ARSyNseq(mydata2, factor = "Tissue", batch = FALSE, norm = "rpkm",  logtransf = FALSE)
myPCA = dat(mydata2corr2, type = "PCA")
par(mfrow = c(1,2))
explo.plot(myPCA, factor = "Tissue")
explo.plot(myPCA, factor = "batch")


###################################################
### code chunk number 30: results
###################################################
mynoiseq = noiseq(mydata, k = 0.5, norm = "rpkm", factor="Tissue", pnr = 0.2, 
                  nss = 5, v = 0.02, lc = 1, replicates = "technical")
head(mynoiseq@results[[1]])


###################################################
### code chunk number 31: NOISeq.Rnw:801-803
###################################################
mynoiseq.tmm = noiseq(mydata, k = 0.5, norm = "tmm", factor="TissueRun", 
                      conditions = c("Kidney_1","Liver_1"), lc = 0, replicates = "technical")


###################################################
### code chunk number 32: NOISeq.Rnw:825-827
###################################################
myresults <- noiseq(mydata, factor = "Tissue", k = NULL, norm="n", pnr = 0.2, 
                    nss = 5, v = 0.02, lc = 1, replicates = "no")


###################################################
### code chunk number 33: NOISeq.Rnw:879-881
###################################################
mynoiseqbio = noiseqbio(mydata, k = 0.5, norm = "rpkm", factor="Tissue", lc = 1, r = 20, adj = 1.5, plot = FALSE,
                        a0per = 0.9, random.seed = 12345, filter = 2)


###################################################
### code chunk number 34: NOISeq.Rnw:926-927
###################################################
head(mynoiseq@results[[1]])


###################################################
### code chunk number 35: NOISeq.Rnw:947-950
###################################################
mynoiseq.deg = degenes(mynoiseq, q = 0.8, M = NULL)
mynoiseq.deg1 = degenes(mynoiseq, q = 0.8, M = "up")
mynoiseq.deg2 = degenes(mynoiseq, q = 0.8, M = "down")


###################################################
### code chunk number 36: fig_summ_expr
###################################################
DE.plot(mynoiseq, q = 0.9, graphic = "expr", log.scale = TRUE)


###################################################
### code chunk number 37: fig_summ_MD
###################################################
DE.plot(mynoiseq, q = 0.8, graphic = "MD")


###################################################
### code chunk number 38: fig_manhattan
###################################################
DE.plot(mynoiseq, chromosomes = c(1,2), log.scale = TRUE,
        join = FALSE, q = 0.8, graphic = "chrom")


###################################################
### code chunk number 39: fig_distrDEG
###################################################
DE.plot(mynoiseq, chromosomes = NULL, q = 0.8, graphic = "distr")


###################################################
### code chunk number 40: session
###################################################
sessionInfo()

Try the NOISeq package in your browser

Any scripts or data that you put into this service are public.

NOISeq documentation built on Nov. 8, 2020, 5:10 p.m.