inst/doc/MAGeCKFlute.R

In MAGeCKFlute: Integrative Analysis Pipeline for Pooled CRISPR Functional Genetic Screens

## ----setup, echo=FALSE, fig.height=6, fig.width=9, dpi=300--------------------
knitr::opts_chunk$set(tidy=FALSE, cache=TRUE,
                      dev="png", message=FALSE, error=FALSE, warning=TRUE)

## ----library, eval=TRUE, message=FALSE----------------------------------------
library(MAGeCKFlute)
library(ggplot2)

## ----gene_summary-------------------------------------------------------------
file1 = file.path(system.file("extdata", package = "MAGeCKFlute"),
                  "testdata/rra.gene_summary.txt")
# Read and visualize the file format
gdata = read.delim(file1, check.names = FALSE)
head(gdata)

## ----readrra------------------------------------------------------------------
gdata = ReadRRA(file1)
head(gdata)

## ----sgrna--------------------------------------------------------------------
file2 = file.path(system.file("extdata", package = "MAGeCKFlute"),
                  "testdata/rra.sgrna_summary.txt")
sdata = read.delim(file2)
head(sdata)

## ----readsgrra----------------------------------------------------------------
sdata = ReadsgRRA(file2)
head(sdata)

## ----FluteRRA1, eval=FALSE----------------------------------------------------
#  FluteRRA(file1, file2, proj="Test", organism="hsa", scale_cutoff = 1, outdir = "./")
#  # Or
#  FluteRRA(gdata, sdata, proj="Test", organism="hsa", scale_cutoff = 1, outdir = "./")

## ----FluteRRA2, eval=FALSE----------------------------------------------------
#  FluteRRA(file1, proj="Test", organism="hsa", outdir = "./")
#  # Or
#  FluteRRA(gdata, proj="Test", organism="hsa", outdir = "./")

## ----DepmapRRA, eval=FALSE----------------------------------------------------
#  FluteRRA(gdata, proj="Test", organism="hsa", incorporateDepmap = TRUE,
#           outdir = "./")

## ----DepmapRRA2, eval=FALSE---------------------------------------------------
#  FluteRRA(gdata, proj="Test", organism="hsa", incorporateDepmap = TRUE,
#           omitEssential = TRUE, outdir = "./")

## ----helpRRA, eval=FALSE------------------------------------------------------
#  ?FluteRRA

## ----mledata------------------------------------------------------------------
file3 = file.path(system.file("extdata", package = "MAGeCKFlute"),
                  "testdata/mle.gene_summary.txt")
# Read and visualize the file format
gdata = read.delim(file3, check.names = FALSE)
head(gdata)

## ----readbeta-----------------------------------------------------------------
gdata = ReadBeta(file3)
head(gdata)

## ----FluteMLE, eval=FALSE-----------------------------------------------------
#  FluteMLE(file3, treatname="plx", ctrlname="dmso", proj="Test", organism="hsa")
#  # Or
#  FluteMLE(gdata, treatname="plx", ctrlname="dmso", proj="Test", organism="hsa")

## ---- eval=FALSE--------------------------------------------------------------
#  ## Take Depmap screen as control
#  FluteMLE(gdata, treatname="plx", ctrlname="Depmap", proj="PLX", organism="hsa", incorporateDepmap = TRUE)

## ---- eval=FALSE--------------------------------------------------------------
#  FluteMLE(gdata, treatname="plx", ctrlname="Depmap", proj="PLX", organism="hsa", incorporateDepmap = TRUE, omitEssential = TRUE)

## ----helpMLE, eval=FALSE------------------------------------------------------
#  ?FluteMLE

## ----CheckCountSummary--------------------------------------------------------
file4 = file.path(system.file("extdata", package = "MAGeCKFlute"),
                  "testdata/countsummary.txt")
countsummary = read.delim(file4, check.names = FALSE)
head(countsummary)

## ----CountQC, fig.height=5, fig.width=4.5-------------------------------------
# Gini index
BarView(countsummary, x = "Label", y = "GiniIndex",
        ylab = "Gini index", main = "Evenness of sgRNA reads")

# Missed sgRNAs
countsummary$Missed = log10(countsummary$Zerocounts)
BarView(countsummary, x = "Label", y = "Missed", fill = "#394E80",
        ylab = "Log10 missed gRNAs", main = "Missed sgRNAs")

# Read mapping
MapRatesView(countsummary)
# Or
countsummary$Unmapped = countsummary$Reads - countsummary$Mapped
gg = reshape2::melt(countsummary[, c("Label", "Mapped", "Unmapped")], id.vars = "Label")
gg$variable = factor(gg$variable, levels = c("Unmapped", "Mapped"))
p = BarView(gg, x = "Label", y = "value", fill = "variable", 
            position = "stack", xlab = NULL, ylab = "Reads", main = "Map ratio")
p + scale_fill_manual(values = c("#9BC7E9", "#1C6DAB"))


## ----CheckRRARes--------------------------------------------------------------
file1 = file.path(system.file("extdata", package = "MAGeCKFlute"),
                  "testdata/rra.gene_summary.txt")
gdata = ReadRRA(file1)
head(gdata)

file2 = file.path(system.file("extdata", package = "MAGeCKFlute"),
                  "testdata/rra.sgrna_summary.txt")
sdata = ReadsgRRA(file2)
head(sdata)

## ----Depmap, eval=FALSE-------------------------------------------------------
#  ## The first run must be time-consuming for downloading Depmap data automatically.
#  depmap_similarity = ResembleDepmap(gdata, symbol = "id", score = "Score")

## ----omitessential, eval=FALSE------------------------------------------------
#  gdata = OmitCommonEssential(gdata)
#  sdata = OmitCommonEssential(sdata, symbol = "Gene")
#  # Compute the similarity with Depmap screens based on subset genes
#  depmap_similarity = ResembleDepmap(gdata, symbol = "id", score = "Score")

## ----selection1, fig.height=4, fig.width=7------------------------------------
gdata$LogFDR = -log10(gdata$FDR)
p1 = ScatterView(gdata, x = "Score", y = "LogFDR", label = "id", model = "volcano", top = 5)
print(p1)

# Or
p2 = VolcanoView(gdata, x = "Score", y = "FDR", Label = "id")
print(p2)


## ----rankrra, fig.height=6, fig.width=4---------------------------------------
gdata$Rank = rank(gdata$Score)
p1 = ScatterView(gdata, x = "Rank", y = "Score", label = "id", 
                 top = 5, auto_cut_y = TRUE, ylab = "Log2FC", 
                 groups = c("top", "bottom"))
print(p1)

## ---- fig.height=4, fig.width=2.5---------------------------------------------
ScatterView(gdata, x = "Rank", y = "Score", label = "id", 
            auto_cut_y = TRUE, groups = c("top", "bottom"), 
            ylab = "Log2FC", toplabels = c("EP300", "NF2"))

## ---- warning=FALSE, fig.height=4, fig.width=6--------------------------------
ScatterView(gdata, x = "Score", y = "Rank", label = "id", 
            auto_cut_x = TRUE, groups = c("left", "right"), 
            xlab = "Log2FC", top = 3)

## ----rankrra2, fig.height=5, fig.width=6--------------------------------------
geneList= gdata$Score
names(geneList) = gdata$id
p2 = RankView(geneList, top = 5, bottom = 10) + xlab("Log2FC")
print(p2)
RankView(geneList, top = 0, bottom = 0, genelist = c("EP300", "NF2")) + xlab("Log2FC")

## ---- warning=FALSE, fig.height=5.5, fig.width=4------------------------------
gdata$Rank = rank(-gdata$Score)
ScatterView(gdata[gdata$Score>0,], x = "Rank", y = "Score", label = "id", 
            auto_cut_y = TRUE, groups = c("top", "bottom"), 
            ylab = "Log2FC", top = 5)

## ----scatter, fig.height=4, fig.width=6---------------------------------------
gdata$RandomIndex = sample(1:nrow(gdata), nrow(gdata))
gdata = gdata[order(-gdata$Score), ]
gg = gdata[gdata$Score>0, ]
p1 = ScatterView(gg, x = "RandomIndex", y = "Score", label = "id",
                 y_cut = CutoffCalling(gdata$Score,2), 
                 groups = "top", top = 5, ylab = "Log2FC")
p1
gg = gdata[gdata$Score<0, ]
p2 = ScatterView(gg, x = "RandomIndex", y = "Score", label = "id",
                 y_cut = CutoffCalling(gdata$Score,2), 
                 groups = "bottom", top = 5, ylab = "Log2FC")
p2

## ----sgRNARank, fig.height=4, fig.width=7-------------------------------------
p2 = sgRankView(sdata, top = 4, bottom = 4)
print(p2)

## ----enrich_rra, fig.height=4, fig.width=9------------------------------------
geneList= gdata$Score
names(geneList) = gdata$id
enrich = EnrichAnalyzer(geneList = geneList[geneList>0.5], 
                        method = "HGT", type = "KEGG")

## ----enrichview---------------------------------------------------------------
EnrichedView(enrich, mode = 1, top = 5)
EnrichedView(enrich, mode = 2, top = 5)

## ----ReadBeta-----------------------------------------------------------------
file3 = file.path(system.file("extdata", package = "MAGeCKFlute"),
                  "testdata/mle.gene_summary.txt")
# Read and visualize the file format
gdata = ReadBeta(file3)
head(gdata)

## ----BatchRemove, fig.height=6, fig.width=9-----------------------------------
##Before batch removal
edata = matrix(c(rnorm(2000, 5), rnorm(2000, 8)), 1000)
colnames(edata) = paste0("s", 1:4)
HeatmapView(cor(edata))

## After batch removal
batchMat = data.frame(sample = colnames(edata), batch = rep(1:2, each = 2))
edata1 = BatchRemove(edata, batchMat)
head(edata1$data)
print(edata1$p)

## ----NormalizeBeta------------------------------------------------------------
ctrlname = "dmso"
treatname = "plx"
gdata_cc = NormalizeBeta(gdata, samples=c(ctrlname, treatname), method="cell_cycle")
head(gdata_cc)

## ----DistributeBeta, fig.height=5, fig.width=8--------------------------------
DensityView(gdata_cc, samples=c(ctrlname, treatname))
ConsistencyView(gdata_cc, ctrlname, treatname)

# Another option MAView
MAView(gdata_cc, ctrlname, treatname)

## ----selection2, fig.height=5, fig.width=7------------------------------------
gdata_cc$Control = rowMeans(gdata_cc[,ctrlname, drop = FALSE])
gdata_cc$Treatment = rowMeans(gdata_cc[,treatname, drop = FALSE])

p1 = ScatterView(gdata_cc, "Control", "Treatment", groups = c("top", "bottom"), auto_cut_diag = TRUE, display_cut = TRUE, toplabels = c("NF1", "NF2", "EP300"))
print(p1)

## ----rank, fig.height=5, fig.width=7------------------------------------------
gdata_cc$Diff = gdata_cc$Treatment - gdata_cc$Control
gdata_cc$Rank = rank(gdata_cc$Diff)
p1 = ScatterView(gdata_cc, x = "Diff", y = "Rank", label = "Gene", 
                 top = 5, model = "rank")
print(p1)

# Or
rankdata = gdata_cc$Treatment - gdata_cc$Control
names(rankdata) = gdata_cc$Gene
RankView(rankdata)

## ----Square, fig.height=6, fig.width=8----------------------------------------
p1 = ScatterView(gdata_cc, x = "dmso", y = "plx", label = "Gene", 
                 model = "ninesquare", top = 5, display_cut = TRUE, force = 2)
print(p1)

## ----Square2, fig.height=6, fig.width=8---------------------------------------
p1 = ScatterView(gdata_cc, x = "dmso", y = "plx", label = "Gene", 
                 model = "ninesquare", top = 5, display_cut = TRUE, 
                 x_cut = c(-1,1), y_cut = c(-1,1))
print(p1)

## ---- fig.height=6, fig.width=8-----------------------------------------------
p2 = SquareView(gdata_cc, label = "Gene", 
                x_cutoff = CutoffCalling(gdata_cc$Control, 2), 
                y_cutoff = CutoffCalling(gdata_cc$Treatment, 2))
print(p2)

## ----EnrichSquare, fig.height=4, fig.width=9----------------------------------
# 9-square groups
Square9 = p1$data
idx=Square9$group=="topcenter"
geneList = Square9$Diff
names(geneList) = Square9$Gene[idx]
universe = Square9$Gene

# Enrichment analysis
kegg1 = EnrichAnalyzer(geneList = geneList, universe = universe)
EnrichedView(kegg1, top = 6, bottom = 0)

## ----pathview2, eval=FALSE----------------------------------------------------
#  genedata = gdata_cc[, c("Control","Treatment")]
#  arrangePathview(genedata, pathways = "hsa01521", organism = "hsa", sub = NULL)

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