MAGeCKFlute - Integrative analysis pipeline for pooled CRISPR functional genetic screens

Binbin Wang, Wubing Zhang, Feizhen Wu, Wei Li & X. Shirley Liu

27 April 2020

Abstract

CRISPR (clustered regularly interspaced short palindrome repeats) coupled with nuclease Cas9 (CRISPR/Cas9) screens represent a promising technology to systematically evaluate gene functions. Data analysis for CRISPR/Cas9 screens is a critical process that includes identifying screen hits and exploring biological functions for these hits in downstream analysis. We have previously developed two algorithms, MAGeCK (Wei Li and Liu. 2014) and MAGeCK-VISPR (Wei Li and Liu. 2015), to analyze CRISPR/Cas9 screen data in various scenarios. These two algorithms allow users to perform quality control, read count generation and normalization, and calculate beta score to evaluate gene selection performance. In downstream analysis, the biological functional analysis is required for understanding biological functions of these identified genes with different screening purposes. Here, We developed MAGeCKFlute for supporting downstream analysis. MAGeCKFlute provides several strategies to remove potential biases within sgRNA-level read counts and gene-level beta scores. The downstream analysis with the package includes identifying essential, non-essential, and target-associated genes, and performing biological functional category analysis, pathway enrichment analysis and protein complex enrichment analysis of these genes. The package also visualizes genes in multiple ways to benefit users exploring screening data. Collectively, MAGeCKFlute enables accurate identification of essential, non-essential, and targeted genes, as well as their related biological functions. This vignette explains the use of the package and demonstrates typical workflows.

1 Quick start
2 Step by step analysis
3 Session info
References

Note: if you use MAGeCKFlute in published research, please cite: Binbin Wang, Mei Wang, Wubing Zhang. “Integrative analysis of pooled CRISPR genetic screens using MAGeCKFlute.” Nature Protocols (2019), doi: 10.1038/s41596-018-0113-7.

1 Quick start

1.1 Load the required packages

library(MAGeCKFlute)
library(ggplot2)

1.2 Downstream analysis of MAGeCK RRA

The MAGeCK (mageck test) uses Robust Rank Aggregation (RRA) for robust identification of CRISPR-screen hits, and outputs the summary results at both sgRNA and gene level. Before performing the downstream analysis, please make sure you have got the gene summary and sgRNA summary results from mageck test. MAGeCKFlute incorporates an example datasets (Ophir Shalem1 2014) for demonstration, shown as below.

1.2.1 gene_summary file (required)

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)

##         id num neg|score neg|p-value neg|fdr neg|rank neg|goodsgrna neg|lfc
## 1   CREBBP   4   1.00000     1.00000       1     8269             0 0.96608
## 2    EP300   4   1.00000     1.00000       1     8270             0 1.02780
## 3      CHD   4   0.99999     0.99999       1     8268             0 0.59265
## 4 C16orf72   4   0.99998     0.99998       1     8267             0 0.82307
## 5   CACNB2   5   0.99863     0.99864       1     8243             0 0.39268
## 6    FGF12   6   0.54609     0.71111       1     4862             1 0.40822
##   pos|score pos|p-value  pos|fdr pos|rank pos|goodsgrna pos|lfc
## 1  3.17e-12  2.9700e-07 0.001650        1             4 0.96608
## 2  3.32e-10  2.9700e-07 0.001650        2             4 1.02780
## 3  1.22e-05  4.4300e-05 0.092203        3             4 0.59265
## 4  2.06e-05  8.0000e-05 0.147965        4             4 0.82307
## 5  4.45e-05  2.1084e-04 0.319082        5             5 0.39268
## 6  5.46e-05  2.7627e-04 0.336280        6             5 0.40822

You can also read the file using ReadRRA in MAGeCKFlute

gdata = ReadRRA(file1)
head(gdata)

##         id   Score      FDR
## 1   CREBBP 0.96608 0.001650
## 2    EP300 1.02780 0.001650
## 3      CHD 0.59265 0.092203
## 4 C16orf72 0.82307 0.147965
## 5   CACNB2 0.39268 0.319082
## 6    FGF12 0.40822 0.336280

Hints: you can also use a data from other analysis, just make sure the three columns (id, Score, and FDR) are avaible in the data.

1.2.2 sgrna_summary file (optional)

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

##     sgrna  Gene control_count treatment_count control_mean treat_mean     LFC
## 1 s_36798   NF2    8917/21204   5020.7/5127.9      15061.0    5074.30 -1.5693
## 2 s_45763 RAB6A 3375.8/3667.7   372.88/357.79       3521.8     365.33 -3.2655
## 3 s_50164   SF1 3657.8/3352.6   453.62/628.28       3505.2     540.95 -2.6937
## 4 s_20780  FDXR 3444.8/3191.1   552.51/555.38       3317.9     553.95 -2.5803
## 5 s_36796   NF2  5492.9/11396   3832.2/3794.6       8444.6    3813.40 -1.1467
## 6 s_36799   NF2 1805.2/5541.7   695.86/882.47       3673.5     789.16 -2.2173
##   control_var adj_var  score       p.low p.high  p.twosided         FDR
## 1    75491000   78871 35.559 2.9804e-277      1 5.9609e-277 1.2732e-272
## 2       42617   15519 25.338 6.1638e-142      1 1.2328e-141 1.9748e-137
## 3       46575   15438 23.857 4.2365e-126      1 8.4731e-126 9.0487e-122
## 4       32179   14520 22.938 9.7946e-117      1 1.9589e-116 1.5690e-112
## 5    17425000   41247 22.803 2.1321e-115      1 4.2641e-115 3.0359e-111
## 6     6980700   16267 22.615 1.5592e-113      1 3.1183e-113 1.9981e-109
##   high_in_treatment
## 1             False
## 2             False
## 3             False
## 4             False
## 5             False
## 6             False

You can also read the file using ReadsgRRA in MAGeCKFlute

sdata = ReadsgRRA(file2)
head(sdata)

##     sgrna  Gene     LFC         FDR
## 1 s_36798   NF2 -1.5693 1.2732e-272
## 2 s_45763 RAB6A -3.2655 1.9748e-137
## 3 s_50164   SF1 -2.6937 9.0487e-122
## 4 s_20780  FDXR -2.5803 1.5690e-112
## 5 s_36796   NF2 -1.1467 3.0359e-111
## 6 s_36799   NF2 -2.2173 1.9981e-109

1.2.3 Run the FluteRRA pipeline

Run the downstream analysis pipeline with both gene summary and sgrna summary

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

Run the downstream analysis pipeline with only gene summary file

FluteRRA(file1, proj="Test", organism="hsa",
         scale_cutoff = 1, outdir = "./")
# Or
FluteRRA(gdata, proj="Test", organism="hsa",
         scale_cutoff = 1, outdir = "./")

Incorporate Depmap data into analysis

FluteRRA(gdata, proj="Test", organism="hsa", incorporateDepmap = TRUE,
         outdir = "./")

Omit common essential genes in the analysis

FluteRRA(gdata, proj="Test", organism="hsa", incorporateDepmap = TRUE,
         omitEssential = TRUE, outdir = "./")

Hints: all figures and intermediate data are saved into local directory “./MAGeCKFlute_Test/”, and all figures are integrated into file “FluteRRA_Test.pdf”.

For more available parameters in FluteRRA, please read the documentation

?FluteRRA

1.3 Downstream analysis of MAGeCK MLE

The MAGeCK-VISPR (mageck mle) computes beta scores and the associated statistics for all genes in multiple conditions. The beta score describes how the gene is selected: a positive beta score indicates a positive selection, and a negative beta score indicates a negative selection. Before using FluteMLE, you should first get gene summary result from MAGeCK-VISPR (mageck mle). MAGeCKFlute incorporates an example datasets for demonstration.

1.3.1 gene_summary file (required)

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)

##       Gene sgRNA dmso|beta   dmso|z dmso|p-value dmso|fdr dmso|wald-p-value
## 1     FEZ1     6 -0.045088 -0.66798      0.79649  0.97939        5.0415e-01
## 2      TNN     6  0.094325  1.36120      0.34176  0.89452        1.7344e-01
## 3     OAS2     8 -0.271210 -4.76860      0.46995  0.93572        1.8555e-06
## 4    JOSD1     4  0.092888  1.08610      0.34538  0.89543        2.7745e-01
## 5      CFH     6 -0.113230 -1.71070      0.95069  0.99513        8.7140e-02
## 6 C11orf68     4 -0.163690 -2.26570      0.77504  0.97694        2.3472e-02
##   dmso|wald-fdr  plx|beta    plx|z plx|p-value plx|fdr plx|wald-p-value
## 1    6.3060e-01 -0.036721 -0.54346     0.81604 0.98345       5.8681e-01
## 2    2.8578e-01  0.065533  0.94344     0.47309 0.93207       3.4546e-01
## 3    1.4126e-05 -0.289010 -5.07170     0.40411 0.90933       3.9431e-07
## 4    4.0754e-01  0.094913  1.10840     0.39207 0.90504       2.6769e-01
## 5    1.6606e-01 -0.060018 -0.90751     0.90090 0.98876       3.6414e-01
## 6    5.7351e-02 -0.094403 -1.30760     0.97492 0.99833       1.9103e-01
##   plx|wald-fdr
## 1   6.9940e-01
## 2   4.7400e-01
## 3   3.5296e-06
## 4   3.9013e-01
## 5   4.9438e-01
## 6   2.9934e-01

You can also read beta scores from the data using ReadBeta in MAGeCKFlute

gdata = ReadBeta(file3)
head(gdata)

##       Gene      dmso       plx
## 1     FEZ1 -0.045088 -0.036721
## 2      TNN  0.094325  0.065533
## 3     OAS2 -0.271210 -0.289010
## 4    JOSD1  0.092888  0.094913
## 5      CFH -0.113230 -0.060018
## 6 C11orf68 -0.163690 -0.094403

Hints: you can also run FluteMLE using other data, in which the first column is “Gene”, and other columns represent samples.

1.3.2 Run the FluteMLE pipeline

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

1.3.2.1 Incorporate Depmap data into analysis

If your data only include one condition, you can take Depmap screens as control.

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

If you are not interested in common essential genes, you can omit them in the analysis by setting a parameter “omitEssential”

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

Hint: All pipeline results are written into local directory “./MAGeCKFlute_Test/”, and all figures are integrated into file “FluteMLE_Test.pdf”.

For more available parameters in FluteMLE, please read the documentation

?FluteMLE

2 Step by step analysis

2.1 Section I: Quality control

2.1.1 Input data

MAGeCK/MAGeCK-VISPR outputs a count summary file, which summarizes some basic QC scores at raw count level, including map ratio, Gini index, and NegSelQC. MAGeCKFlute incorporates an example datasets (Ophir Shalem1 2014) for demonstration.

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

##                                   File    Label    Reads   Mapped Percentage
## 1 ../data/GSC_0131_Day23_Rep1.fastq.gz day23_r1 62818064 39992777     0.6366
## 2  ../data/GSC_0131_Day0_Rep2.fastq.gz  day0_r2 47289074 31709075     0.6705
## 3  ../data/GSC_0131_Day0_Rep1.fastq.gz  day0_r1 51190401 34729858     0.6784
## 4 ../data/GSC_0131_Day23_Rep2.fastq.gz day23_r2 58686580 37836392     0.6447
##   TotalsgRNAs Zerocounts GiniIndex NegSelQC NegSelQCPval
## 1       64076         57   0.08510        0            1
## 2       64076         17   0.07496        0            1
## 3       64076         14   0.07335        0            1
## 4       64076         51   0.08587        0            1
##   NegSelQCPvalPermutation NegSelQCPvalPermutationFDR NegSelQCGene
## 1                       1                          1            0
## 2                       1                          1            0
## 3                       1                          1            0
## 4                       1                          1            0

2.1.2 Visualize the QC results

# 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 = data.table::melt(countsummary[, c("Label", "Mapped", "Unmapped")], id.vars = "Label")

## Warning in data.table::melt(countsummary[, c("Label", "Mapped", "Unmapped")], :
## The melt generic in data.table has been passed a data.frame and will attempt
## to redirect to the relevant reshape2 method; please note that reshape2 is
## deprecated, and this redirection is now deprecated as well. To continue using
## melt methods from reshape2 while both libraries are attached, e.g. melt.list,
## you can prepend the namespace like reshape2::melt(countsummary[, c("Label",
## "Mapped", "Unmapped")]). In the next version, this warning will become an error.

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

2.2 Section II: Downstream analysis of MAGeCK RRA

For CRISPR/Cas9 screens with two experimental conditions, MAGeCK-RRA is available for identification of essential genes. In MAGeCK-RRA results, the sgRNA summary and gene summary file summarizes the statistical significance of positive selections and negative selections at sgRNA level and gene level.

2.2.1 Read the required data

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

##         id   Score      FDR
## 1   CREBBP 0.96608 0.001650
## 2    EP300 1.02780 0.001650
## 3      CHD 0.59265 0.092203
## 4 C16orf72 0.82307 0.147965
## 5   CACNB2 0.39268 0.319082
## 6    FGF12 0.40822 0.336280

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

##     sgrna  Gene     LFC         FDR
## 1 s_36798   NF2 -1.5693 1.2732e-272
## 2 s_45763 RAB6A -3.2655 1.9748e-137
## 3 s_50164   SF1 -2.6937 9.0487e-122
## 4 s_20780  FDXR -2.5803 1.5690e-112
## 5 s_36796   NF2 -1.1467 3.0359e-111
## 6 s_36799   NF2 -2.2173 1.9981e-109

2.2.2 Compute the similarity between the CRISPR screen with Depmap screens

depmap_similarity = ResembleDepmap(gdata, symbol = "id", score = "Score")
head(depmap_similarity)

##         estimate p.value
## OVCAR8 0.6207002       0
## RKO    0.6178466       0
## A2780  0.6176061       0
## HS578T 0.6163724       0
## MEL202 0.6129839       0
## SF295  0.6124781       0

2.2.3 Omit common essential genes from the data

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")
head(depmap_similarity)

##            estimate      p.value
## TE11      0.2538626 1.153518e-85
## HCC827GR5 0.2511691 2.810572e-86
## ESS1      0.2481296 3.503664e-84
## OVCAR8    0.2475420 8.837071e-84
## NCIH2452  0.2469750 2.152940e-83
## HSC1      0.2466611 3.521411e-83

2.2.4 Visualization of negative selections and positive selections

2.2.4.1 Volcano plot

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)

2.2.4.2 Rank plot

Rank all the genes based on their scores and label genes in the rank plot.

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)

Label interested hits using parameter toplabels (in ScatterView) and genelist (in RankView).

ScatterView(gdata, x = "Rank", y = "Score", label = "id", 
            auto_cut_y = TRUE, groups = c("top", "bottom"), 
            ylab = "Log2FC", toplabels = c("EP300", "NF2"))

Plot Log2FC at x-axis

ScatterView(gdata, x = "Score", y = "Rank", label = "id", 
            auto_cut_x = TRUE, groups = c("left", "right"), 
            xlab = "Log2FC", top = 3)

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

Only plot positive selection

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)

2.2.4.3 Dot plot

Visualize negative and positive selected genes separately.

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

2.2.4.4 sgRankView - visualize the rank of sgRNAs targeting top selected genes.

p2 = sgRankView(sdata, top = 4, bottom = 4)
print(p2)

2.2.5 Enrichment analysis

For more information about functional enrichment analysis in MAGeCKFlute, please read the MAGeCKFlute_enrichment document, in which we introduce all the available options and methods.

geneList= gdata$Score
names(geneList) = gdata$id
enrich = EnrichAnalyzer(geneList = geneList[geneList>0.5], 
                        method = "HGT", type = "KEGG")

2.2.5.1 Visualization of enrichment results

EnrichedView(enrich, mode = 1, top = 5)

EnrichedView(enrich, mode = 2, top = 5)

2.3 Section III: Downstream analysis of MAGeCK MLE

2.3.1 read required data

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

##       Gene      dmso       plx
## 1     FEZ1 -0.045088 -0.036721
## 2      TNN  0.094325  0.065533
## 3     OAS2 -0.271210 -0.289010
## 4    JOSD1  0.092888  0.094913
## 5      CFH -0.113230 -0.060018
## 6 C11orf68 -0.163690 -0.094403

2.3.2 Batch effect removal (not recommended)

Is there batch effects? This is a commonly asked question before perform later analysis. In our package, we provide HeatmapView to ensure whether the batch effect exists in data and use BatchRemove to remove easily if same batch samples cluster together.

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

##            s1       s2       s3       s4
## [1,] 6.154871 6.044254 6.190671 6.117837
## [2,] 6.147682 5.925400 5.569169 6.451310
## [3,] 6.846471 6.739190 7.228395 6.264156
## [4,] 7.745407 6.144626 7.051242 6.580360
## [5,] 6.281997 5.586695 6.332966 5.462174
## [6,] 5.860856 5.429238 4.689217 6.406219

print(edata1$p)

2.3.3 Normalization of beta scores

It is difficult to control all samples with a consistent cell cycle in a CRISPR screen experiment with multi conditions. Besides, beta score among different states with an inconsistent cell cycle is incomparable. So it is necessary to do the normalization when comparing the beta scores in different conditions. Essential genes are those genes that are indispensable for its survival. The effect generated by knocking out these genes in different cell types is consistent. Based on this, we developed the cell cycle normalization method to shorten the gap of the cell cycle in different conditions.

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

##       Gene        dmso         plx
## 1     FEZ1 -0.05913258 -0.05702550
## 2      TNN  0.12370654  0.10176880
## 3     OAS2 -0.35568991 -0.44881511
## 4    JOSD1  0.12182193  0.14739417
## 5      CFH -0.14850031 -0.09320434
## 6 C11orf68 -0.21467823 -0.14660217

2.3.4 Distribution of all gene beta scores

After normalization, the distribution of beta scores in different conditions should be similar. We can evaluate the distribution of beta scores using the function ‘DensityView’, and ‘ConsistencyView’.

DensityView(gdata_cc, samples=c(ctrlname, treatname))

## Warning in data.table::melt(beta, id = NULL): The melt generic in data.table
## has been passed a data.frame and will attempt to redirect to the relevant
## reshape2 method; please note that reshape2 is deprecated, and this redirection
## is now deprecated as well. To continue using melt methods from reshape2 while
## both libraries are attached, e.g. melt.list, you can prepend the namespace like
## reshape2::melt(beta). In the next version, this warning will become an error.

ConsistencyView(gdata_cc, ctrlname, treatname)

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

2.3.5 Positive selection and negative selection

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)

2.3.5.1 Rank plot - label top hits

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)

2.3.5.2 Nine-square scatter plot - identify treatment-associated genes

p1 = ScatterView(gdata_cc, x = "dmso", y = "plx", label = "Gene", 
                 model = "ninesquare", top = 5, display_cut = TRUE, force = 2)
print(p1)

Customize the cutoff

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)

p2 = SquareView(gdata_cc, label = "Gene", 
                x_cutoff = CutoffCalling(gdata_cc$Control, 2), 
                y_cutoff = CutoffCalling(gdata_cc$Treatment, 2))
print(p2)

2.3.6 Functional analysis for treatment-associated genes

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

Also, pathway visualization can be done using function KeggPathwayView (Luo et al. 2013).

genedata = gdata_cc[, c("Control","Treatment")]
arrangePathview(genedata, pathways = "hsa01521", organism = "hsa", sub = NULL)

3 Session info

sessionInfo()

## R version 4.0.0 (2020-04-24)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.4 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.11-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.11-bioc/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] ggplot2_3.3.0     MAGeCKFlute_1.8.0 BiocStyle_2.16.0 
## 
## loaded via a namespace (and not attached):
##   [1] fastmatch_1.1-0        BiocFileCache_1.12.0   plyr_1.8.6            
##   [4] igraph_1.2.5           splines_4.0.0          BiocParallel_1.22.0   
##   [7] pathview_1.28.0        sva_3.36.0             urltools_1.7.3        
##  [10] digest_0.6.25          htmltools_0.4.0        GOSemSim_2.14.0       
##  [13] viridis_0.5.1          magick_2.3             GO.db_3.10.0          
##  [16] magrittr_1.5           memoise_1.1.0          limma_3.44.0          
##  [19] Biostrings_2.56.0      annotate_1.66.0        graphlayouts_0.7.0    
##  [22] matrixStats_0.56.0     askpass_1.1            enrichplot_1.8.0      
##  [25] prettyunits_1.1.1      colorspace_1.4-1       blob_1.2.1            
##  [28] rappdirs_0.3.1         ggrepel_0.8.2          xfun_0.13             
##  [31] dplyr_0.8.5            crayon_1.3.4           RCurl_1.98-1.2        
##  [34] jsonlite_1.6.1         graph_1.66.0           scatterpie_0.1.4      
##  [37] genefilter_1.70.0      survival_3.1-12        glue_1.4.0            
##  [40] polyclip_1.10-0        gtable_0.3.0           zlibbioc_1.34.0       
##  [43] XVector_0.28.0         Rgraphviz_2.32.0       BiocGenerics_0.34.0   
##  [46] scales_1.1.0           DOSE_3.14.0            msigdbr_7.0.1         
##  [49] pheatmap_1.0.12        DBI_1.1.0              edgeR_3.30.0          
##  [52] Rcpp_1.0.4.6           viridisLite_0.3.0      xtable_1.8-4          
##  [55] progress_1.2.2         gridGraphics_0.5-0     bit_1.1-15.2          
##  [58] europepmc_0.3          stats4_4.0.0           httr_1.4.1            
##  [61] fgsea_1.14.0           RColorBrewer_1.1-2     ellipsis_0.3.0        
##  [64] pkgconfig_2.0.3        XML_3.99-0.3           farver_2.0.3          
##  [67] dbplyr_1.4.3           locfit_1.5-9.4         ggplotify_0.0.5       
##  [70] tidyselect_1.0.0       labeling_0.3           rlang_0.4.5           
##  [73] reshape2_1.4.4         AnnotationDbi_1.50.0   munsell_0.5.0         
##  [76] tools_4.0.0            downloader_0.4         RSQLite_2.2.0         
##  [79] ggridges_0.5.2         evaluate_0.14          stringr_1.4.0         
##  [82] yaml_2.2.1             knitr_1.28             bit64_0.9-7           
##  [85] tidygraph_1.1.2        purrr_0.3.4            KEGGREST_1.28.0       
##  [88] dendextend_1.13.4      ggraph_2.0.2           nlme_3.1-147          
##  [91] KEGGgraph_1.48.0       DO.db_2.9              xml2_1.3.2            
##  [94] biomaRt_2.44.0         compiler_4.0.0         curl_4.3              
##  [97] png_0.1-7              ggsignif_0.6.0         tibble_3.0.1          
## [100] tweenr_1.0.1           stringi_1.4.6          lattice_0.20-41       
## [103] Matrix_1.2-18          ggsci_2.9              vctrs_0.2.4           
## [106] pillar_1.4.3           lifecycle_0.2.0        BiocManager_1.30.10   
## [109] triebeard_0.3.0        data.table_1.12.8      cowplot_1.0.0         
## [112] bitops_1.0-6           qvalue_2.20.0          R6_2.4.1              
## [115] bookdown_0.18          gridExtra_2.3          IRanges_2.22.0        
## [118] codetools_0.2-16       MASS_7.3-51.6          assertthat_0.2.1      
## [121] openssl_1.4.1          withr_2.2.0            S4Vectors_0.26.0      
## [124] mgcv_1.8-31            parallel_4.0.0         hms_0.5.3             
## [127] clusterProfiler_3.16.0 grid_4.0.0             tidyr_1.0.2           
## [130] rmarkdown_2.1          rvcheck_0.1.8          ggpubr_0.2.5          
## [133] ggforce_0.3.1          Biobase_2.48.0

References

Luo, Weijun, Brouwer, and Cory. 2013. “Pathview: An R/Bioconductor Package for Pathway-Based Data Integration and Visualization.” Bioinformatics 29 (14):1830–1. https://doi.org/10.1093/bioinformatics/btt285.

Ophir Shalem1, *, 2. 2014. “Genome-scale CRISPR-Cas9 knockout screening in human cells.” http://science.sciencemag.org/content/343/6166/84.long.

Wei Li, Han Xu, Johannes Köster, and X. Shirley Liu. 2015. “Quality control, modeling, and visualization of CRISPR screens with MAGeCK-VISPR.” https://genomebiology.biomedcentral.com/articles/10.1186/s13059-015-0843-6.

Wei Li, Tengfei Xiao, Han Xu, and X Shirley Liu. 2014. “MAGeCK enables robust identification of essential genes from genome-scale CRISPR/Cas9 knockout screens.” https://genomebiology.biomedcentral.com/articles/10.1186/s13059-014-0554-4.