## ---- 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)

timepoints <- gsub('^(Control|Death)', '', pData(es)$TREATMENT_NAME)
names(timepoints) <- colnames(es)
es.norm <- ct.normalizeGuides(es, 'FQ', annotation = ann, sampleKey = timepoints, 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)

## ---- fig.width=6, fig.height = 4---------------------------------------------
ct.signalSummary(resultsDF,
                 targets = list('TargetSetA' = c(sample(unique(resultsDF$geneSymbol), 3)),
                                'TargetSetB' = c(sample(unique(resultsDF$geneSymbol), 2))))


## ---- 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:
#  genesetdb <- sparrow::getMSigGeneSetDb(collection = 'h', species = 'human', id.type = 'entrez')
#  
#  # If you have a library that targets elements unevenly (e.g., variable numbers of
#  # elements/promoters per genes), you can conform it via GREAT
#  genesetdb.GREAT <- ct.GREATdb(ann, gsdb = genesetdb)
#  
#  ct.seas(resultsDF, gdb = genesetdb)
#  #ct.seas(resultsDF, gdb = genesetdb.GREAT)

## ---- 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, 'enrich')
str(ROC)

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

## ---- fig.width=6, fig.height = 6, warning=FALSE------------------------------
##' tar <-  sample(unique(resultsDF$geneSymbol), 20)
##' res <- ct.targetSetEnrichment(resultsDF, tar)

targetsTest <- ct.targetSetEnrichment(resultsDF, essential.genes, enrich = FALSE)
str(targetsTest)

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

