---
title: "Example Workflow For Processing a Single Pooled Screen"
author: "Russell Bainer"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{Example_Workflow_gCrisprTools}
%\VignetteEngine{knitr::knitr}
%\VignetteEncoding{UTF-8}
---
### Example Workflow For Processing a Single Screen
This is an example workflow for processing a pooled screening eperiment using the provided sample data. See the various manpages for additional visualization options and algorithmic details.
Note that what follows describes a very basic analysis. If you are considering integrating the results of many different screen contrasts or even different experiments and/or technologies, refer to the `Advanced Screen Analysis: Contrast Comparisons` vignette.
Load dependencies and data
```{r setup, eval = TRUE}
suppressPackageStartupMessages(library(Biobase))
suppressPackageStartupMessages(library(limma))
suppressPackageStartupMessages(library(gCrisprTools))
data("es", package = "gCrisprTools")
data("ann", package = "gCrisprTools")
data("aln", package = "gCrisprTools")
knitr::opts_chunk$set(message = FALSE, fig.width = 8, fig.height = 8, warning = FALSE)
```
Make a sample key, structured as a factor with control samples in the first level
```{r, eval = TRUE}
sk <- relevel(as.factor(pData(es)$TREATMENT_NAME), "ControlReference")
names(sk) <- row.names(pData(es))
```
Generate a contrast of interest using voom/limma; pairing replicates is a good idea if that information is available.
```{r, eval = TRUE}
design <- model.matrix(~ 0 + REPLICATE_POOL + TREATMENT_NAME, pData(es))
colnames(design) <- gsub('TREATMENT_NAME', '', colnames(design))
contrasts <-makeContrasts(DeathExpansion - ControlExpansion, levels = design)
```
Optionally, trim of trace reads from the unnormalized object (see man page for details)
```{r, eval = TRUE, fig.width = 8, fig.height = 10}
es <- ct.filterReads(es, trim = 1000, sampleKey = sk)
```
Normalize, convert to a voom object, and generate a contrast
```{r, eval = TRUE, fig.width = 6, fig.height = 12}
es <- ct.normalizeGuides(es, method = "scale", plot.it = TRUE) #See man page for other options
vm <- voom(exprs(es), design)
fit <- lmFit(vm, design)
fit <- contrasts.fit(fit, contrasts)
fit <- eBayes(fit)
```
Edit the annotation file if you used `ct.filterReads` above
```{r, eval = TRUE}
ann <- ct.prepareAnnotation(ann, fit, controls = "NoTarget")
```
Summarize gRNA signals to identify target genes of interest
```{r, eval = TRUE}
resultsDF <-
ct.generateResults(
fit,
annotation = ann,
RRAalphaCutoff = 0.1,
permutations = 1000,
scoring = "combined",
permutation.seed = 2
)
```
#### Alternative Annotations
In some cases, reagents might target multiple known elements (e.g., gRNAs in a CRISPRi library that target multiple promoters of the same gene). In such cases, you can specify this via the `alt.annotation` argument to `ct.generateResults()`. Alternative annotations are supplied as a list of character vectors named for the reagents.
```{r, eval = TRUE}
# Create random alternative target associations
altann <- sapply(ann$ID,
function(x){
out <- as.character(ann$geneSymbol)[ann$ID %in% x]
if(runif(1) < 0.01){out <- c(out, sample(as.character(ann$geneSymbol), size = 1))}
return(out)
}, simplify = FALSE)
resultsDF <-
ct.generateResults(
fit,
annotation = ann,
RRAalphaCutoff = 0.1,
permutations = 1000,
scoring = "combined",
alt.annotation = altann,
permutation.seed = 2
)
```
Optionally, just load an example results object for testing purposes (trimming out reads as necessary)
```{r, eval = TRUE}
data("fit", package = "gCrisprTools")
data("resultsDF", package = "gCrisprTools")
fit <- fit[(row.names(fit) %in% row.names(ann)),]
resultsDF <- resultsDF[(row.names(resultsDF) %in% row.names(ann)),]
targetResultDF <- ct.simpleResult(resultsDF) #For a simpler target-level result object
```
### Quality Control
gCrisprTools contains a variety of pooled screen-specific quality control and visualization tools (see man pages for details):
```{r, eval = TRUE}
ct.alignmentChart(aln, sk)
ct.rawCountDensities(es, sk)
```
Visualize gRNA abundance distributions
```{r, eval = TRUE}
ct.gRNARankByReplicate(es, sk)
ct.gRNARankByReplicate(es, sk, annotation = ann, geneSymb = "NoTarget") #Show locations of NTC gRNAs
```
Visualize control guide behavior across conditions
```{r, eval = TRUE}
ct.viewControls(es, ann, sk, normalize = FALSE)
ct.viewControls(es, ann, sk, normalize = TRUE)
```
Visualize GC bias across samples, or within an experimental contrast
```{r, eval = TRUE, fig.width = 8, fig.height = 12}
ct.GCbias(es, ann, sk)
ct.GCbias(fit, ann, sk)
```
View most variable gRNAs/Genes (as % of sequencing library)
```{r, eval = TRUE}
ct.stackGuides(es,
sk,
plotType = "gRNA",
annotation = ann,
nguides = 40)
```
```{r, eval = TRUE}
ct.stackGuides(es,
sk,
plotType = "Target",
annotation = ann)
```
```{r, eval = TRUE}
ct.stackGuides(es,
sk,
plotType = "Target",
annotation = ann,
subset = names(sk)[grep('Expansion', sk)])
```
View a CDF of genes/guides
```{r, eval = TRUE}
ct.guideCDF(es, sk, plotType = "gRNA")
ct.guideCDF(es, sk, plotType = "Target", annotation = ann)
```
#### Target-Level Visualization and Analysis
View the overall enrichment and depletion trends identified in the screen:
```{r, eval = TRUE}
ct.contrastBarchart(resultsDF)
```
View top enriched/depleted candidates
```{r, eval = TRUE}
ct.topTargets(fit,
resultsDF,
ann,
targets = 10,
enrich = TRUE)
ct.topTargets(fit,
resultsDF,
ann,
targets = 10,
enrich = FALSE)
```
View the behavior of reagents targeting a particular gene of interest
```{r, eval = TRUE, fig.width = 8, fig.height = 10}
ct.viewGuides("Target1633", fit, ann)
ct.gRNARankByReplicate(es, sk, annotation = ann, geneSymb = "Target1633")
```
Observe the effects detected for sets of targets within a screen contrast
```{r, eval = TRUE}
ct.signalSummary(resultsDF,
targets = list('TargetSetA' = c(sample(unique(resultsDF$geneSymbol), 3)),
'TargetSetB' = c(sample(unique(resultsDF$geneSymbol), 2))))
```
You could test a known gene set for enrichment within target candidates:
```{r, eval = TRUE}
data("essential.genes", package = "gCrisprTools")
ct.targetSetEnrichment(resultsDF, essential.genes)
```
Or optionally add a visualization:
```{r rocprc, eval = TRUE}
ROC <- ct.ROC(resultsDF, essential.genes, direction = "deplete")
PRC <- ct.PRC(resultsDF, essential.genes, direction = "deplete")
show(ROC) # show(PRC) is equivalent for the PRC analysis
```
Or alternatively you could test for ontological enrichment within the depleted/enriched targets via the `sparrow` package:
```{r, eval = TRUE, warning=FALSE, message = FALSE}
#Create a geneset database using one of the many helper functions
genesetdb <- sparrow::getMSigGeneSetDb(collection = 'h', species = 'human', id.type = 'entrez')
ct.seas(resultsDF, gdb = genesetdb)
# If you have a library that targets elements unevenly (e.g., variable numbers of
# elements/promoters per genes), you can conform it via `sparrow::convertIdentifiers()`
genesetdb.GREAT <- sparrow::convertIdentifiers(genesetdb,
from = 'geneID',
to = 'geneSymbol',
xref = ann)
ct.seas(resultsDF, gdb = genesetdb.GREAT)
```
See the `Contrast_Comparisons` vignette for more advanced use cases of gCrisprTools and extension to complex experiments and study designs.
Finally, you can make reports in a directory of interest:
```{r, eval = FALSE}
path2report <- #Make a report of the whole experiment
ct.makeReport(fit = fit,
eset = es,
sampleKey = sk,
annotation = ann,
results = resultsDF,
aln = aln,
outdir = ".")
path2QC <- #Or one focusing only on experiment QC
ct.makeQCReport(es,
trim = 1000,
log2.ratio = 0.05,
sampleKey = sk,
annotation = ann,
aln = aln,
identifier = 'Crispr_QC_Report',
lib.size = NULL
)
path2Contrast <- #Or Contrast-specific one
ct.makeContrastReport(eset = es,
fit = fit,
sampleKey = sk,
results = resultsDF,
annotation = ann,
comparison.id = NULL,
identifier = 'Crispr_Contrast_Report')
```
```{r}
sessionInfo()
```