inst/doc/gCrisprTools_Vignette.R

In gCrisprTools: Suite of Functions for Pooled Crispr Screen QC and Analysis

## ---- eval = FALSE------------------------------------------------------------
#  if (!requireNamespace("BiocManager", quietly=TRUE))
#      install.packages("BiocManager")
#  BiocManager::install("gCrisprTools")

## ---- message=FALSE, warning=FALSE--------------------------------------------
library(Biobase)
library(limma)
library(gCrisprTools)

## -----------------------------------------------------------------------------
data("es", package = "gCrisprTools")
es
head(exprs(es))

## -----------------------------------------------------------------------------
data("ann", package = "gCrisprTools")
head(ann)

## -----------------------------------------------------------------------------
sk <- relevel(as.factor(pData(es)$TREATMENT_NAME), "ControlReference")
names(sk) <- row.names(pData(es))
sk

## ---- fig.width = 7, fig.height = 5-------------------------------------------
data("aln", package = "gCrisprTools")
head(aln)
ct.alignmentChart(aln, sk)

## ---- fig.width=6, fig.height = 8---------------------------------------------
es.floor <- ct.filterReads(es, read.floor = 30, sampleKey = sk)
es <- ct.filterReads(es, trim = 1000, log2.ratio = 4, sampleKey = sk)

##Convenience function for conforming the annotation object to exclude the trimmed gRNAs
ann <- ct.prepareAnnotation(ann, es, controls = "NoTarget")

## ---- fig.width=6, fig.height = 8---------------------------------------------

es <- ct.normalizeGuides(es, 'scale', annotation = ann, sampleKey = sk, plot.it = TRUE)
es.norm <- ct.normalizeGuides(es, 'slope', annotation = ann, sampleKey = sk, plot.it = TRUE)
es.norm <- ct.normalizeGuides(es, 'controlScale', annotation = ann, sampleKey = sk, plot.it = TRUE, geneSymb = 'NoTarget')
es.norm <- ct.normalizeGuides(es, 'controlSpline', annotation = ann, sampleKey = sk, plot.it = TRUE, geneSymb = 'NoTarget')

## ---- eval= FALSE-------------------------------------------------------------
#  #Not run:
#  path2QC <- ct.makeQCReport(es,
#                   trim = 1000,
#                   log2.ratio = 0.05,
#                   sampleKey = sk,
#                   annotation = ann,
#                   aln = aln,
#                   identifier = 'Crispr_QC_Report',
#                   lib.size = NULL
#                   )

## ---- fig.width=6, fig.height=6-----------------------------------------------
ct.rawCountDensities(es, sk)

## ---- fig.width=6, fig.height = 6---------------------------------------------
ct.gRNARankByReplicate(es, sk)  #Visualization of gRNA abundance distribution

## ---- fig.width=6, fig.height = 6---------------------------------------------
ct.gRNARankByReplicate(es, sk, annotation = ann, geneSymb = "Target1633")

## ---- fig.width=6, fig.height = 6---------------------------------------------
ct.viewControls(es, ann, sk, normalize = FALSE, geneSymb = 'NoTarget')

## ---- fig.width=6, fig.height = 4---------------------------------------------
ct.guideCDF(es, sk, plotType = "gRNA")

## -----------------------------------------------------------------------------
design <- model.matrix(~ 0 + REPLICATE_POOL + TREATMENT_NAME, pData(es))
colnames(design) <- gsub('TREATMENT_NAME', '', colnames(design))
contrasts <-makeContrasts(DeathExpansion - ControlExpansion, levels = design)

vm <- voom(exprs(es), design)

fit <- lmFit(vm, design)
fit <- contrasts.fit(fit, contrasts)
fit <- eBayes(fit)

## ---- message=FALSE, warning=FALSE--------------------------------------------
resultsDF <-
  ct.generateResults(
    fit,
    annotation = ann,
    RRAalphaCutoff = 0.1,
    permutations = 1000,
    scoring = "combined"
  )

## ---- fig.width=6, fig.height = 6---------------------------------------------
ct.topTargets(fit,
              resultsDF,
              ann,
              targets = 10,
              enrich = TRUE)

## ---- fig.width=6, fig.height = 8---------------------------------------------
ct.stackGuides(
  es,
  sk,
  plotType = "Target",
  annotation = ann,
  subset = names(sk)[grep('Expansion', sk)]
)

## ---- fig.width=6, fig.height = 4---------------------------------------------
ct.viewGuides("Target1633", fit, ann)

## ---- eval= FALSE-------------------------------------------------------------
#  #Not run:
#  path2Contrast <-
#    ct.makeContrastReport(eset = es,
#                          fit = fit,
#                          sampleKey = sk,
#                          results = resultsDF,
#                          annotation = ann,
#                          comparison.id = NULL,
#                          identifier = 'Crispr_Contrast_Report')

## ---- eval=FALSE--------------------------------------------------------------
#  #Not run:
#  path2report <-
#    ct.makeReport(fit = fit,
#                  eset = es,
#                  sampleKey = sk,
#                  annotation = ann,
#                  results = resultsDF,
#                  aln = aln,
#                  outdir = ".")

## ---- eval= FALSE-------------------------------------------------------------
#  #Not run:
#  enrichmentResults <-
#    ct.PantherPathwayEnrichment(
#      resultsDF,
#      pvalue.cutoff = 0.01,
#      enrich = TRUE,
#      organism = 'mouse'
#    )
#  
#  > head(enrichmentResults)   #Note: Pathway names have been edited for display purposes.
#                           PATHWAY nGenes sigGenes expected     odds            p        FDR
#  1 EGF receptor signaling pathway    200       14 5.240550 3.332647 0.0004498023 0.03958260
#  2          FGF signaling pathway    230       14 5.949958 2.869779 0.0016284304 0.07165094
#  3     Insulin/MAP kinase cascade    138        9 3.714916 2.785331 0.0101632822 0.20272465
#  4          CCKR signaling map ST    331       15 8.211061 2.148459 0.0126368744 0.20272465
#  5               p38 MAPK pathway    145        9 3.891333 2.641968 0.0135928688 0.20272465
#  6              B cell activation    204       11 5.336192 2.368705 0.0146442541 0.20272465

## ---- fig.width=6, fig.height = 6, warning=FALSE------------------------------
data("essential.genes", package = "gCrisprTools")  #Artificial list created for demonstration
data("resultsDF", package = "gCrisprTools")
ROC <- ct.ROC(resultsDF, essential.genes, stat = "deplete.p")
str(ROC)

## ---- fig.width=6, fig.height = 6, warning=FALSE------------------------------
PRC <- ct.PRC(resultsDF, essential.genes, stat = "deplete.p")
str(PRC)

## ---- fig.width=6, fig.height = 6, warning=FALSE------------------------------
targetsTest <- ct.targetSetEnrichment(resultsDF, essential.genes, enrich = FALSE)
str(targetsTest)

## -----------------------------------------------------------------------------
sessionInfo()