vignette.R
In MPRAnalyze: Statistical Analysis of MPRA data

## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)

## ----load and examine data, include=TRUE--------------------------------------
library(MPRAnalyze)
data("ChrEpi")
summary(ce.colAnnot)
head(ce.colAnnot)

## ----init object, include=TRUE------------------------------------------------
obj <- MpraObject(dnaCounts = ce.dnaCounts, rnaCounts = ce.rnaCounts, 
                  dnaAnnot = ce.colAnnot, rnaAnnot = ce.colAnnot, 
                  controls = ce.control)

## ----library size estimation--------------------------------------------------
## If the library factors are different for the DNA and RNA data, separate 
## estimation of these factors is needed. We can also change the estimation 
## method (Upper quartile by default)
obj <- estimateDepthFactors(obj, lib.factor = c("batch", "condition"),
                            which.lib = "dna", 
                            depth.estimator = "uq")
obj <- estimateDepthFactors(obj, lib.factor = c("condition"),
                            which.lib = "rna", 
                            depth.estimator = "uq")

## In this case, the factors are the same - each combination of batch and 
## condition is a single library, and we'll use the default estimation
obj <- estimateDepthFactors(obj, lib.factor = c("batch", "condition"),
                            which.lib = "both")

## ----quant model fit----------------------------------------------------------
obj <- analyzeQuantification(obj = obj, 
                              dnaDesign = ~ barcode + batch + condition,
                              rnaDesign = ~ condition)

## ----quant extract and viz----------------------------------------------------
##extract alpha values from the fitted model
alpha <- getAlpha(obj, by.factor = "condition")

##visualize the estimates
boxplot(alpha)

## ----quant test and viz-------------------------------------------------------
## test 
res.epi <- testEmpirical(obj = obj, 
                               statistic = alpha$MT)
summary(res.epi)

res.chr <- testEmpirical(obj = obj,
                               statistic = alpha$WT)
par(mfrow=c(2,2))

hist(res.epi$pval.mad, main="episomal, all")
hist(res.epi$pval.mad[ce.control], main="episomal, controls")
hist(res.chr$pval.mad, main="chromosomal, all")
hist(res.chr$pval.mad[ce.control], main="chromosomal, controls")

par(mfrow=c(1,1))

## ----comp fit-----------------------------------------------------------------
obj <- analyzeComparative(obj = obj, 
                           dnaDesign = ~ barcode + batch + condition, 
                           rnaDesign = ~ condition, 
                           reducedDesign = ~ 1)

## ----comp lrt-----------------------------------------------------------------
res <- testLrt(obj)

head(res)
summary(res)

## ----foldchange---------------------------------------------------------------

## plot log Fold-Change
plot(density(res$logFC))
## plot volcano
plot(res$logFC, -log10(res$pval))