Integrative analysis pipeline for pooled CRISPR functional genetic screens - MAGeCKFlute
WubingZhang, Feizhen Wu, and Binbin Wang
5 October 2018
Abstract
Genome-wide 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 and MAGeCK-VISPR, 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, 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, utilizing the data provided through MAGeCK and MAGeCK-VISPR. 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 and pathway enrichment analysis of these genes. The package also visualizes genes in the context of pathways 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. MAGeCKFlute package version: 1.0.1- How to get help for MAGeCKFlute
- Input data
- Quick start
- Section I: Run pipeline from gene summary file generated by MAGeCK MLE
- Gene summary file
- Read beta scores
- Batch effect removal
- Normalization of beta scores
- Distribution of all gene beta scores
- Estimate cell cycle time by linear fitting
- Positive selection and negative selection
- Functional analysis of selected genes
- Identify treatment-associated genes using 9-square model
- Functional analysis for treatment-associated genes
- Section II: Run pipeline from gene summary file generated by MAGeCK RRA
- Session info
- References
Note: if you use MAGeCKFlute in published research, please cite:
How to get help for MAGeCKFlute
Any and all MAGeCKFlute questions should be posted to the Bioconductor support site, which serves as a searchable knowledge base of questions and answers:
https://support.bioconductor.org
Posting a question and tagging with “MAGeCKFlute” will automatically send an alert to the package authors to respond on the support site. See the first question in the list of Frequently Asked Questions (FAQ) for information about how to construct an informative post.
You should not email your question to the package authors, as we will just reply that the question should be posted to the Bioconductor support site.
Input data
As input, the MAGeCKFlute package expects gene summary file as obtained by running commands mageck test or mageck mle in MAGeCK (Wei Li and Liu. 2014) and MAGeCK-VISPR (Wei Li and Liu. 2015), which are developed by our lab previously, to analyze CRISPR/Cas9 screen data in different scenarios(Tim Wang 2014, Hiroko Koike-Yusa (2014), Ophir Shalem1 (2014), Luke A.Gilbert (2014), Silvana Konermann (2015)). Both algorithms use negative binomial models to model the variances of sgRNAs, and use Robust Rank Aggregation (for MAGeCK) or maximum likelihood framework (for MAGeCK-VISPR) for a robust identification of selected genes.
MAGeCK-MLE can be used to analyze CRISPR screen data from multi-conditioned experiments. MAGeCK-MLE also normalizes the data across multiple samples, making them comparable to each other. The most important ouput of MAGeCK MLE is “gene_summary” file, which includes the beta scores of multiple conditions and the associated statistics. The ‘beta score’ for each gene describes how the gene is selected: a positive beta score indicates a positive selection, and a negative beta score indicates a negative selection.
MAGeCK RRA allows for the comparison between two experimental conditions. It can identify genes and sgRNAs are significantly selected between the two conditions. The most important output of MAGeCK RRA is the file “gene_summary.txt”. MAGeCK RRA will output both the negative score and positive score for each gene. A smaller score indicates higher gene importance. MAGeCK RRA will also output the statistical value for the scores of each gene. Genes that are significantly positively and negatively selected can be identified based on the p-value or FDR.
Quick start
Here we show the most basic steps for integrative analysis pipeline using gene summary file. Before using MAGeCKFlute, analysing CRISPR/Cas9 screen data using MAGeCK RRA (in MAGeCK (Wei Li and Liu. 2014)) or MAGeCK MLE (in MAGeCK-VISPR (Wei Li and Liu. 2015)) is neccessary, which result in the generation of the gene summary file.
To run MAGeCKFlute pipeline, we need gene summary file generated by running MAGeCK RRA or MAGeCK MLE. MAGeCKFlute package provides two example data, one is MLE_Data and the other is RRA_Data. We will work with them in this document.
Downstream analysis pipeline for MAGeCK MLE
library(MAGeCKFlute)
##Load gene summary data in MAGeCK MLE results
data("MLE_Data")
##Run the downstream analysis pipeline for MAGeCK MLE
FluteMLE(MLE_Data, ctrlname=c("D7_R1", "D7_R2"), treatname=c("PLX7_R1","PLX7_R2"), prefix="BRAF_", organism="hsa")All pipeline results are written into local directory “./BRAF_Flute_Results/”, and all figures are integrated into file “BRAF_Flute.mle_summary.pdf”.
Downstream analysis pipeline for MAGeCK RRA
##Load gene summary data in MAGeCK RRA results
data("RRA_Data")
##Run the downstream analysis pipeline for MAGeCK RRA
FluteRRA(RRA_Data, prefix="BRAF", organism="hsa")All pipeline results are written into local directory “./BRAF_Flute_Results/” too, and all figures are integrated into file “BRAF_Flute.rra_summary.pdf”.
Section I: Run pipeline from gene summary file generated by MAGeCK MLE
Gene summary file
As the input of MAGeCKFlute package, the gene summary file in MAGeCK-MLE results includes beta scores of all genes in multiple condition samples. Use function ‘data’ to load the dataset, and have a look at the file with a text editor to see how it is formatted.
## Gene D7_R1.beta D7_R2.beta PLX7_R1.beta PLX7_R2.beta
## 1 SMPD2 -0.458770 -0.514690 -0.67770 -0.510360
## 2 MIXL1 0.094571 0.122780 -0.14383 -0.119670
## 3 ENPEP -0.082621 -0.082305 -0.11271 -0.061796
## 4 VAV2 -0.024270 -0.034077 0.14537 -0.015095
## 5 C5orf24 0.452310 0.226940 0.15809 0.191940
## 6 TLCD1 -0.883320 -0.779930 -0.80811 -0.854900Read beta scores
Then, extract beta scores of control and treatment samples from the gene summary table(can be a file path of ‘gene_summary’ or data frame), and have a look at the beta score matrix.
data("MLE_Data")
gene_summary = MLE_Data
ctrlname = c("D7_R1", "D7_R2")
treatname = c("PLX7_R1", "PLX7_R2")
#Read beta scores from gene summary table in MAGeCK MLE results
dd=ReadBeta(gene_summary, organism="hsa")
head(dd)## D7_R1 D7_R2 PLX7_R1 PLX7_R2
## SMPD2 -0.458770 -0.514690 -0.67770 -0.510360
## MIXL1 0.094571 0.122780 -0.14383 -0.119670
## ENPEP -0.082621 -0.082305 -0.11271 -0.061796
## VAV2 -0.024270 -0.034077 0.14537 -0.015095
## C5orf24 0.452310 0.226940 0.15809 0.191940
## TLCD1 -0.883320 -0.779930 -0.80811 -0.854900# Transform gene symbol to Entrez gene id
dd$Gene = rownames(dd)
tmp = TransGeneID(dd$Gene, "Symbol", "Entrez")
idx = is.na(tmp) | duplicated(tmp)
dd = dd[!idx,]
rownames(dd) = tmp[!idx]
head(dd)## D7_R1 D7_R2 PLX7_R1 PLX7_R2 Gene
## 6610 -0.458770 -0.514690 -0.67770 -0.510360 SMPD2
## 83881 0.094571 0.122780 -0.14383 -0.119670 MIXL1
## 2028 -0.082621 -0.082305 -0.11271 -0.061796 ENPEP
## 7410 -0.024270 -0.034077 0.14537 -0.015095 VAV2
## 134553 0.452310 0.226940 0.15809 0.191940 C5orf24
## 116238 -0.883320 -0.779930 -0.80811 -0.854900 TLCD1Batch effect removal
Is there batch effect? This is a common asked question before perform later analysis. In our package, we provide HeatmapView to ensure whether batch effect exists in data, and use BatchRemove to remove easily if same batch samples cluster together.
batchMat = data.frame(samples = c(ctrlname, treatname), batch = c(1,2,1,2), cov = c(1,1,2,2))
dd1 = BatchRemove(dd[,c(ctrlname, treatname)], batchMat)$data## Standardizing Data across genes
Normalization of beta scores
It is difficult to control all samples with a consistent cell cycle in an CRISPR screen experiment with multi conditions. Besides, beta score among different conditions with an inconsistent cell cycle are incomparable. So it is necessary to do the normalization when comparing the beta scores in different conditions. Essential genesare thosegenes 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 cell cycle normalization method to shorten the gap of cell cycle in different conditions. In addition, a previous normalization method called loess normalization is available in this package.(Laurent Gautier 2004)
dd_essential = NormalizeBeta(dd, samples=c(ctrlname, treatname), method="cell_cycle")
head(dd_essential)## D7_R1 D7_R2 PLX7_R1 PLX7_R2 Gene
## 6610 -0.50335740 -0.58521416 -0.8882277 -0.66804545 SMPD2
## 83881 0.10376226 0.13960363 -0.1885108 -0.15664433 MIXL1
## 2028 -0.09065085 -0.09358264 -0.1477234 -0.08088905 ENPEP
## 7410 -0.02662878 -0.03874632 0.1905292 -0.01975889 VAV2
## 134553 0.49626956 0.25803591 0.2072007 0.25124352 C5orf24
## 116238 -0.96916899 -0.88679803 -1.0591497 -1.11903765 TLCD1## D7_R1 D7_R2 PLX7_R1 PLX7_R2 Gene
## 6610 -0.43596609 -0.50733383 -0.6981720 -0.52004811 SMPD2
## 83881 0.09677572 0.13618156 -0.1525602 -0.12654604 MIXL1
## 2028 -0.08013862 -0.07305255 -0.1188903 -0.06735051 ENPEP
## 7410 -0.01964304 -0.02103090 0.1336521 -0.02105016 VAV2
## 134553 0.45732275 0.25111644 0.1422490 0.17859184 C5orf24
## 116238 -0.83430196 -0.77372518 -0.8411143 -0.87711860 TLCD1Distribution 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 function ‘ViolinView’, ‘DensityView’, and ‘DensityDiffView’.
#we can also use function 'MAView' to evaluate the data quality of normalized
#beta score profile.
MAView(dd_essential, ctrlname, treatname, cex=1, main="Cell cycle normalized")
Estimate cell cycle time by linear fitting
After normalization, the cell cycle time in different condition should be almost consistent. Here we use linear fitting to estimate the cell cycle time, and use function CellCycleView to view the cell cycle time of all samples.
##Fitting beta score of all genes
CellCycleView(dd_essential, ctrlname, treatname, main="Cell cycle normalized")
Positive selection and negative selection
The function ScatterView can group all genes into three groups, positive selection genes (GroupA), negative selection genes (GroupB), and others, and visualize these three grouped genes in scatter plot. We can also use function RankView to rank the beta score deviation between control and treatment and mark top selected genes in figure.
## Add column of 'diff'
dd_essential$Control = rowMeans(dd_essential[,ctrlname])
dd_essential$Treatment = rowMeans(dd_essential[,treatname])
rankdata = dd_essential$Treatment - dd_essential$Control
names(rankdata) = rownames(dd_essential)
p2 = RankView(rankdata, main="Cell cycle normalized")
print(p2)
Functional analysis of selected genes
For gene set enrichment analysis, we provide five methods in this package, including “ORT”(Over-Representing Test (Guangchuang Yu and He. 2012)), “GSEA”(Gene Set Enrichment Analysis (Aravind Subramanian and Mesirov. 2005)), “DAVID” (Huang da W 2009), “GOstats” (S. Falcon 2007), and “HGT”(HyperGemetric test). And gene ontology(GO) term (Consortium. 2014) and Kyoto encyclopedia of genes and genomes (KEGG) pathway (Minoru Kanehisa 2014) enrichment analysis are available in this package. The enrichment analysis can be done easily using function enrichment_analysis, which return a list containing gridPlot (ggplot object) and enrichRes (enrichResult instance). Alternatively, you can do enrichment analysis using function enrich.ORT for “ORT”, enrich.GSE for GSEA, enrich.DAVID for DAVID, enrich.GOstats for GOstats, and enrich.HGT for “HGT”, which return an enrichResult instance. Function EnrichedView and EnrichedGSEView (for enrich.GSE) can be used to generate gridPlot from enrichReseasily, as shown below.
## Get information of positive and negative selection genes
groupAB = p1$data
## select positive selection genes
idx1=groupAB$group=="up"
genes=rownames(groupAB)[idx1]
geneList=groupAB$diff[idx1]
names(geneList)=genes
universe=rownames(groupAB)
## Do enrichment analysis using HGT method
keggA = enrichment_analysis(geneList = genes, universe = universe, method = "HGT",
type = "KEGG", organism = "human", plotTitle = "Positive selection")
#same as
kegg_A = enrich.HGT(genes, universe, type = "KEGG", organism = "human")
keggA_grid = EnrichedView(kegg_A@result, plotTitle = "Positive selection")
## look at the results
head(keggA$enrichRes@result)## ID Description pvalue p.adjust
## 17 hsa00270 Cysteine and methionine metabolism 0.002413521 0.1675247
## 31 hsa00480 Glutathione metabolism 0.004713313 0.1675247
## 68 hsa00970 Aminoacyl-tRNA biosynthesis 0.005318244 0.1675247
## 85 hsa04122 Sulfur relay system 0.004230677 0.1675247
## nLogpvalue GeneRatio BgRatio
## 17 0.7759212 14/43 43/29909
## 31 0.7759212 15/50 50/29909
## 68 0.7759212 13/42 42/29909
## 85 0.7759212 4/8 8/29909
## geneID
## 17 191/58478/10993/875/4357/23382/10768/3948/6723/2805/4190/1788/2937/635
## 31 2941/2944/2539/257202/2876/6723/79017/2879/5226/119391/2936/2938/2882/2937/9027
## 68 10352/3735/25973/10667/51067/57176/92935/57505/123283/23395/5859/124454/123263
## 85 27304/9054/348180/4357
## geneName
## 17 AHCY/ENOPH1/SDS/CBS/MPST/AHCYL2/AHCYL1/LDHC/SRM/GOT1/MDH1/DNMT3A/GSS/BHMT
## 31 GSTA4/GSTM1/G6PD/GPX6/GPX1/SRM/GGCT/GPX4/PGD/GSTO2/GSR/GSTA1/GPX7/GSS/NAT8
## 68 WARS2/KARS/PARS2/FARS2/YARS2//MARS2/AARS2/TARSL2/LARS2/QARS/EARS2/MTFMT
## 85 MOCS3/NFS1/CTU2/MPST
## Count
## 17 14
## 31 15
## 68 13
## 85 4
## ID Description pvalue
## 14 hsa00190 Oxidative phosphorylation 5.302278e-09
## 191 hsa04714 Thermogenesis 2.243007e-08
## 90 hsa03010 Ribosome 2.901736e-07
## 249 hsa05012 Parkinson disease 6.529432e-07
## 251 hsa05016 Huntington disease 3.817661e-06
## 228 hsa04932 Non-alcoholic fatty liver disease (NAFLD) 1.373651e-04
## p.adjust nLogpvalue GeneRatio BgRatio
## 14 1.707333e-06 5.767682 40/101 101/29909
## 191 3.611241e-06 5.442344 61/190 190/29909
## 90 3.114530e-05 4.506607 40/114 114/29909
## 249 5.256193e-05 4.279329 38/109 109/29909
## 251 2.458574e-04 3.609317 49/163 163/29909
## 228 6.318793e-03 2.199366 39/138 138/29909
## geneID
## 14 9377/1327/1337/155066/4716/4726/1329/54539/6389/4708/23545/1340/4719/4702/6392/4695/525/4717/4715/55967/4713/51079/4709/4701/7385/29796/534/4720/10063/7386/1355/9114/1339/4728/4724/4718/526/1353/7381/4722
## 191 9377/6195/60/84699/8110/1327/1337/4716/4726/1329/57521/196883/54539/51103/6389/4708/1340/4719/4702/6602/6392/4695/4717/5606/4715/5568/55967/8289/114/4713/51079/6598/4709/57104/107/4701/6199/79133/7385/29796/493753/7248/64223/2778/4720/10063/7386/1355/6605/1339/4728/4724/4718/90639/1353/90993/7381/1268/7249/4722/112
## 90 51021/6231/51116/65003/6230/6208/64979/6135/63931/124995/51318/29093/64968/51264/29088/65005/28998/6183/6191/55052/64981/6193/64969/79590/6234/6181/63875/55168/6182/6150/51187/6146/6168/11222/6209/54460/6169/64963/54948/51073
## 249 9377/1327/1337/4716/7317/4726/1329/7332/54539/6389/4708/1340/4719/4702/6392/4695/4717/4715/5568/55967/6531/4713/51079/7054/4709/135/4701/7385/29796/4720/7386/1339/4728/4724/4718/7381/7345/4722
## 251 9377/7019/5330/84699/581/160/1327/7157/1337/4716/4726/1329/6648/54539/6389/4708/1340/4719/4702/6392/4695/2876/4717/4715/55967/6667/4713/5434/51079/1173/4709/4701/1387/7385/29796/6647/4720/7386/4899/5436/1175/1339/4728/4724/4718/90993/7381/1767/4722
## 228 9377/581/1327/1337/4716/4726/1329/54539/6389/4708/1340/4719/4702/7494/6392/4695/5465/4717/4715/55967/4713/51079/1050/4709/4701/9451/356/6256/7385/29796/7186/4720/7386/1339/4728/4724/4718/7381/4722
## geneName
## 14 COX5A/COX4I1/COX6A1/ATP6V0E2/NDUFB10/NDUFS6/COX5B/NDUFB11/SDHA/NDUFB2/ATP6V0A2/COX6B1/NDUFS1/NDUFA8/SDHD/NDUFA2/ATP6V1B1/NDUFC1/NDUFB9/NDUFA12/NDUFB7/NDUFA13/NDUFB3/NDUFA7/UQCRC2/UQCR10/ATP6V1G2/NDUFS2/COX17/UQCRFS1/COX15/ATP6V0D1/COX6A2/NDUFS8/NDUFS4/NDUFC2/ATP6V1B2/COX11/UQCRB/NDUFS3
## 191 COX5A/RPS6KA1/ACTB/CREB3L3/DPF3/COX4I1/COX6A1/NDUFB10/NDUFS6/COX5B/RPTOR/ADCY4/NDUFB11/NDUFAF1/SDHA/NDUFB2/COX6B1/NDUFS1/NDUFA8/SMARCD1/SDHD/NDUFA2/NDUFC1/MAP2K3/NDUFB9/PRKACG/NDUFA12/ARID1A/ADCY8/NDUFB7/NDUFA13/SMARCB1/NDUFB3/PNPLA2/ADCY1/NDUFA7/RPS6KB2/NDUFAF5/UQCRC2/UQCR10/COA5/TSC1/MLST8/GNAS/NDUFS2/COX17/UQCRFS1/COX15/SMARCE1/COX6A2/NDUFS8/NDUFS4/NDUFC2/COX19/COX11/CREB3L1/UQCRB/CNR1/TSC2/NDUFS3/ADCY6
## 90 MRPS16/RPS26/MRPS2/MRPL11/RPS25/RPS14/MRPL36/RPL11/MRPS14/MRPL10/MRPL35/MRPL22/MRPS6/MRPL27/MRPL15/MRPL9/MRPL13/MRPS12/RPS4X/MRPL20/MRPL34/RPS5/MRPS5/MRPL24/RPS28/RPLP2/MRPL17/MRPS18A/MRPL12/MRPL23/RSL24D1/RPL22/RPL37A/MRPL3/RPS15/MRPS21/RPL38/MRPS11/MRPL16/MRPL4
## 249 COX5A/COX4I1/COX6A1/NDUFB10/UBA1/NDUFS6/COX5B/UBE2L3/NDUFB11/SDHA/NDUFB2/COX6B1/NDUFS1/NDUFA8/SDHD/NDUFA2/NDUFC1/NDUFB9/PRKACG/NDUFA12/SLC6A3/NDUFB7/NDUFA13/TH/NDUFB3/ADORA2A/NDUFA7/UQCRC2/UQCR10/NDUFS2/UQCRFS1/COX6A2/NDUFS8/NDUFS4/NDUFC2/UQCRB/UCHL1/NDUFS3
## 251 COX5A/TFAM/PLCB2/CREB3L3/BAX/AP2A1/COX4I1/TP53/COX6A1/NDUFB10/NDUFS6/COX5B/SOD2/NDUFB11/SDHA/NDUFB2/COX6B1/NDUFS1/NDUFA8/SDHD/NDUFA2/GPX1/NDUFC1/NDUFB9/NDUFA12/SP1/NDUFB7/POLR2E/NDUFA13/AP2M1/NDUFB3/NDUFA7/CREBBP/UQCRC2/UQCR10/SOD1/NDUFS2/UQCRFS1/NRF1/POLR2G/AP2S1/COX6A2/NDUFS8/NDUFS4/NDUFC2/CREB3L1/UQCRB/DNAH5/NDUFS3
## 228 COX5A/BAX/COX4I1/COX6A1/NDUFB10/NDUFS6/COX5B/NDUFB11/SDHA/NDUFB2/COX6B1/NDUFS1/NDUFA8/XBP1/SDHD/NDUFA2/PPARA/NDUFC1/NDUFB9/NDUFA12/NDUFB7/NDUFA13/CEBPA/NDUFB3/NDUFA7/EIF2AK3/FASLG/RXRA/UQCRC2/UQCR10/TRAF2/NDUFS2/UQCRFS1/COX6A2/NDUFS8/NDUFS4/NDUFC2/UQCRB/NDUFS3
## Count
## 14 40
## 191 61
## 90 40
## 249 38
## 251 49
## 228 39
## Do enrichment analysis using GSEA method
gseA = enrichment_analysis(geneList = geneList, method = "GSEA", type = "KEGG",
organism="human", plotTitle="Positive selection")
#same as
gse_A = enrich.GSE(geneList, type = "KEGG", organism = "human")
gseA_grid = EnrichedGSEView(gse_A@result, plotTitle = "Positive selection")
head(gseA$enrichRes@result)## ID Description setSize
## hsa05016 hsa05016 Huntington disease 49
## hsa05010 hsa05010 Alzheimer disease 41
## hsa00190 hsa00190 Oxidative phosphorylation 40
## hsa04932 hsa04932 Non-alcoholic fatty liver disease (NAFLD) 39
## hsa05012 hsa05012 Parkinson disease 38
## hsa05169 hsa05169 Epstein-Barr virus infection 27
## enrichmentScore NES pvalue p.adjust qvalues rank
## hsa05016 0.4514230 1.942223 0.001019368 0.05356516 0.05102324 615
## hsa05010 0.4923908 2.050555 0.001029866 0.05356516 0.05102324 615
## hsa00190 0.4758918 1.969399 0.001031992 0.05356516 0.05102324 615
## hsa04932 0.4931197 2.028598 0.001033058 0.05356516 0.05102324 615
## hsa05012 0.5027018 2.053564 0.001036269 0.05356516 0.05102324 615
## hsa05169 0.5280043 2.006460 0.001053741 0.05356516 0.05102324 630
## leading_edge
## hsa05016 tags=49%, list=23%, signal=38%
## hsa05010 tags=54%, list=23%, signal=42%
## hsa00190 tags=52%, list=23%, signal=41%
## hsa04932 tags=54%, list=23%, signal=42%
## hsa05012 tags=55%, list=23%, signal=43%
## hsa05169 tags=59%, list=24%, signal=46%
## core_enrichment
## hsa05016 7019/1327/1387/9377/7381/160/5434/6392/4715/4708/4702/4701/1329/7386/29796/7385/6389/1337/1340/4716/4719/6647/4709/55967
## hsa05010 3028/1327/8883/9377/7381/6392/4715/4708/4702/4701/1329/7386/29796/7385/6389/1337/1340/4716/2906/4719/4709/55967
## hsa00190 10063/1327/9377/7381/6392/4715/4708/4702/1355/4701/1329/7386/29796/7385/6389/1337/1340/4716/4719/4709/55967
## hsa04932 1050/1327/9377/7381/6392/4715/4708/4702/4701/1329/7186/7386/29796/7385/6389/1337/1340/4716/4719/4709/55967
## hsa05012 1327/9377/7381/6392/4715/4708/4702/7332/4701/1329/7386/29796/7385/6389/1337/1340/4716/4719/4709/135/55967
## hsa05169 1460/1387/5608/5434/3312/3305/171568/5704/51728/7186/3106/6693/10621/890/9612/3133
## geneName
## hsa05016 TFAM/COX4I1/CREBBP/COX5A/UQCRB/AP2A1/POLR2E/SDHD/NDUFB9/NDUFB2/NDUFA8/NDUFA7/COX5B/UQCRFS1/UQCR10/UQCRC2/SDHA/COX6A1/COX6B1/NDUFB10/NDUFS1/SOD1/NDUFB3/NDUFA12
## hsa05010 HSD17B10/COX4I1/NAE1/COX5A/UQCRB/SDHD/NDUFB9/NDUFB2/NDUFA8/NDUFA7/COX5B/UQCRFS1/UQCR10/UQCRC2/SDHA/COX6A1/COX6B1/NDUFB10/GRIN2D/NDUFS1/NDUFB3/NDUFA12
## hsa00190 COX17/COX4I1/COX5A/UQCRB/SDHD/NDUFB9/NDUFB2/NDUFA8/COX15/NDUFA7/COX5B/UQCRFS1/UQCR10/UQCRC2/SDHA/COX6A1/COX6B1/NDUFB10/NDUFS1/NDUFB3/NDUFA12
## hsa04932 CEBPA/COX4I1/COX5A/UQCRB/SDHD/NDUFB9/NDUFB2/NDUFA8/NDUFA7/COX5B/TRAF2/UQCRFS1/UQCR10/UQCRC2/SDHA/COX6A1/COX6B1/NDUFB10/NDUFS1/NDUFB3/NDUFA12
## hsa05012 COX4I1/COX5A/UQCRB/SDHD/NDUFB9/NDUFB2/NDUFA8/UBE2L3/NDUFA7/COX5B/UQCRFS1/UQCR10/UQCRC2/SDHA/COX6A1/COX6B1/NDUFB10/NDUFS1/NDUFB3/ADORA2A/NDUFA12
## hsa05169 CSNK2B/CREBBP/MAP2K6/POLR2E/HSPA8/HSPA1L/POLR3H/PSMC4/POLR3K/TRAF2/HLA-B/SPN/POLR3F/CCNA2/NCOR2/HLA-E
## ID Description setSize
## hsa04714 hsa04714 Thermogenesis 61
## hsa05016 hsa05016 Huntington disease 49
## hsa04932 hsa04932 Non-alcoholic fatty liver disease (NAFLD) 39
## hsa00190 hsa00190 Oxidative phosphorylation 40
## hsa05010 hsa05010 Alzheimer disease 41
## hsa05012 hsa05012 Parkinson disease 38
## enrichmentScore NES pvalue p.adjust qvalues rank
## hsa04714 0.4228475 1.893086 0.001005025 0.02911956 0.0261445 615
## hsa05016 0.4514230 1.946578 0.001014199 0.02911956 0.0261445 615
## hsa04932 0.4931197 2.049043 0.001023541 0.02911956 0.0261445 615
## hsa00190 0.4758918 1.981589 0.001024590 0.02911956 0.0261445 615
## hsa05010 0.4923908 2.052409 0.001024590 0.02911956 0.0261445 615
## hsa05012 0.5027018 2.069244 0.001027749 0.02911956 0.0261445 615
## leading_edge
## hsa04714 tags=46%, list=23%, signal=36%
## hsa05016 tags=49%, list=23%, signal=38%
## hsa04932 tags=54%, list=23%, signal=42%
## hsa00190 tags=52%, list=23%, signal=41%
## hsa05010 tags=54%, list=23%, signal=42%
## hsa05012 tags=55%, list=23%, signal=43%
## core_enrichment
## hsa04714 6598/10063/1327/9377/7381/6392/4715/79133/4708/4702/1355/4701/1329/8110/7386/29796/7385/6389/1337/1340/6199/60/4716/4719/4709/8289/51103/55967
## hsa05016 7019/1327/1387/9377/7381/160/5434/6392/4715/4708/4702/4701/1329/7386/29796/7385/6389/1337/1340/4716/4719/6647/4709/55967
## hsa04932 1050/1327/9377/7381/6392/4715/4708/4702/4701/1329/7186/7386/29796/7385/6389/1337/1340/4716/4719/4709/55967
## hsa00190 10063/1327/9377/7381/6392/4715/4708/4702/1355/4701/1329/7386/29796/7385/6389/1337/1340/4716/4719/4709/55967
## hsa05010 3028/1327/8883/9377/7381/6392/4715/4708/4702/4701/1329/7386/29796/7385/6389/1337/1340/4716/2906/4719/4709/55967
## hsa05012 1327/9377/7381/6392/4715/4708/4702/7332/4701/1329/7386/29796/7385/6389/1337/1340/4716/4719/4709/135/55967
## geneName
## hsa04714 SMARCB1/COX17/COX4I1/COX5A/UQCRB/SDHD/NDUFB9/NDUFAF5/NDUFB2/NDUFA8/COX15/NDUFA7/COX5B/DPF3/UQCRFS1/UQCR10/UQCRC2/SDHA/COX6A1/COX6B1/RPS6KB2/ACTB/NDUFB10/NDUFS1/NDUFB3/ARID1A/NDUFAF1/NDUFA12
## hsa05016 TFAM/COX4I1/CREBBP/COX5A/UQCRB/AP2A1/POLR2E/SDHD/NDUFB9/NDUFB2/NDUFA8/NDUFA7/COX5B/UQCRFS1/UQCR10/UQCRC2/SDHA/COX6A1/COX6B1/NDUFB10/NDUFS1/SOD1/NDUFB3/NDUFA12
## hsa04932 CEBPA/COX4I1/COX5A/UQCRB/SDHD/NDUFB9/NDUFB2/NDUFA8/NDUFA7/COX5B/TRAF2/UQCRFS1/UQCR10/UQCRC2/SDHA/COX6A1/COX6B1/NDUFB10/NDUFS1/NDUFB3/NDUFA12
## hsa00190 COX17/COX4I1/COX5A/UQCRB/SDHD/NDUFB9/NDUFB2/NDUFA8/COX15/NDUFA7/COX5B/UQCRFS1/UQCR10/UQCRC2/SDHA/COX6A1/COX6B1/NDUFB10/NDUFS1/NDUFB3/NDUFA12
## hsa05010 HSD17B10/COX4I1/NAE1/COX5A/UQCRB/SDHD/NDUFB9/NDUFB2/NDUFA8/NDUFA7/COX5B/UQCRFS1/UQCR10/UQCRC2/SDHA/COX6A1/COX6B1/NDUFB10/GRIN2D/NDUFS1/NDUFB3/NDUFA12
## hsa05012 COX4I1/COX5A/UQCRB/SDHD/NDUFB9/NDUFB2/NDUFA8/UBE2L3/NDUFA7/COX5B/UQCRFS1/UQCR10/UQCRC2/SDHA/COX6A1/COX6B1/NDUFB10/NDUFS1/NDUFB3/ADORA2A/NDUFA12
For enriched pathways, we can use function KeggPathwayView to visualize the beta score level in control and treatment on pathway map.(Weijun Luo 2013)
genedata = dd_essential[,c("Control","Treatment")]
keggID = keggA$enrichRes@result$ID[1]
#The pathway map will be located on current workspace
KeggPathwayView(gene.data = genedata, pathway.id = keggID, species="hsa")## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.
## Warning in structure(x$children, class = "XMLNodeList"): Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
## Consider 'structure(list(), *)' instead.##Read the figure into R
pngname=paste0(keggID, ".pathview.multi.png")
grid.arrange(grid::rasterGrob(png::readPNG(pngname)))
## [1] TRUE TRUE TRUEIdentify treatment-associated genes using 9-square model
Considering the difference of beta scores in control and treatment sample, we developed a 9-square model, which group all genes into several subgroups. Among these subgroups, four subgroup genes are treatment-associated, which correspond to specific functions. Group1 and Group3 genes are not selected in control sample, while they are significantly selected in treatment sample, so they may related to drug resistance. Group2 and Group4 genes are selected in control, but they are not selected in treatment, so maybe these genes are associated with drug targets.
In this package, use function SquareView can select these four treatment-associated subgroup genes and view them in 9-Square scatter plot.
Functional analysis for treatment-associated genes
Same as section above. We can do enrichment analysis for treatment-associated genes.
##Get information of treatment-associated genes
Square9 = p3$data
##==select group1 genes in 9-Square
idx=Square9$group=="Group1"
genes=rownames(Square9)[idx]
universe=rownames(Square9)
#====KEGG_enrichment=====
kegg1=enrichment_analysis(geneList = genes, universe = universe,
type = "KEGG", plotTitle = "KEGG: Group1")
## look at the results
head(kegg1$enrichRes@result)## ID Description
## hsa00130 hsa00130 Ubiquinone and other terpenoid-quinone biosynthesis
## hsa00030 hsa00030 Pentose phosphate pathway
## hsa00360 hsa00360 Phenylalanine metabolism
## hsa00100 hsa00100 Steroid biosynthesis
## hsa00051 hsa00051 Fructose and mannose metabolism
## hsa03460 hsa03460 Fanconi anemia pathway
## GeneRatio BgRatio pvalue p.adjust qvalue
## hsa00130 3/325 10/6513 0.01136634 0.4740760 0.4740760
## hsa00030 5/325 28/6513 0.01142352 0.4740760 0.4740760
## hsa00360 3/325 16/6513 0.04251604 0.9193596 0.9193596
## hsa00100 3/325 18/6513 0.05759572 0.9193596 0.9193596
## hsa00051 4/325 31/6513 0.06627724 0.9193596 0.9193596
## hsa03460 5/325 46/6513 0.07738612 0.9193596 0.9193596
## geneID Count geneName
## hsa00130 84274/27235/51805 3 COQ5/COQ2/COQ3
## hsa00030 226/7086/5226/5213/51071 5 /TKT/PGD/PFKM/DERA
## hsa00360 314/4128/5053 3 AOC2/MAOA/PAH
## hsa00100 120227/8435/1595 3 CYP2R1/SOAT2/CYP51A1
## hsa00051 226/3099/5372/5213 4 /HK2/PMM1/PFKM
## hsa03460 9894/4292/8940/3280/5429 5 TELO2/MLH1/TOP3B/HES1/POLH
Also, pathway visualization can be done using function KeggPathwayView, the same as section above.
genedata = dd_essential[, c("Control","Treatment")]
keggID = kegg1$enrichRes@result$ID[1]
KeggPathwayView(gene.data = genedata, pathway.id = keggID, species="hsa")
##Read the figure into R
pngname=paste0(keggID, ".pathview.multi.png")
grid.arrange(grid::rasterGrob(png::readPNG(pngname)))
file.remove(paste0(keggID, c(".pathview.multi.png", ".png", ".xml")))Section II: Run pipeline from gene summary file generated by MAGeCK RRA
Gene summary file
For experiments with two experimental conditions, we recommend to use MAGeCK-RRA to identify essential genes from CRISPR/Cas9 knockout screens and tests the statistical significance of each observed change between two conditions. Gene summary file in MAGeCK-RRA results summarize the statistical significance of positive selection and negative selection. Use function ‘data’ to load the dataset, and have a look at the file with a text editor to see how it is formatted.
## id neg.fdr pos.fdr
## 1 CMIP 0.001650 0.977480
## 2 MCL1 0.001650 1.000000
## 3 ITGB2 0.001650 0.999967
## 4 SCN2A 0.003713 0.886655
## 5 GRB2 0.048680 0.984303
## 6 EGFR 0.048680 0.993803Read FDR of negative selection and positive selection
Then, extract “neg.fdr” and “pos.fdr” from the gene summary table.
## Official neg.fdr pos.fdr ENTREZID
## 1 CMIP 0.001650 0.977480 80790
## 2 MCL1 0.001650 1.000000 4170
## 3 ITGB2 0.001650 0.999967 3689
## 4 SCN2A 0.003713 0.886655 6326
## 5 GRB2 0.048680 0.984303 2885
## 6 EGFR 0.048680 0.993803 1956Negative selection and positive selection
Take 0.05 as cutoff, get negative selection and positive selection genes and do enrichment analysis on KEGG pathway and GO BP terms.
##Negative selection
universe=dd.rra$ENTREZID
genes = dd.rra[dd.rra$neg.fdr<0.05, "ENTREZID"]
kegg.neg=enrichment_analysis(geneList = genes, universe=universe,
type = "KEGG", plotTitle="KEGG: neg")
bp.neg=enrichment_analysis(geneList = genes, universe=universe,
type = "BP", plotTitle="BP: neg")
print(kegg.neg$gridPlot)
##Positive selection
universe=dd.rra$ENTREZID
genes = dd.rra[dd.rra$pos.fdr<0.05, "ENTREZID"]
kegg.pos=enrichment_analysis(geneList = genes, universe=universe,
type = "KEGG", plotTitle="KEGG: pos")
bp.pos=enrichment_analysis(geneList = genes, universe=universe,
type = "BP", plotTitle="BP: pos")
print(kegg.pos$gridPlot)
Session info
## R version 3.5.1 Patched (2018-07-12 r74967)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.5 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.7-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.7-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] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] MAGeCKFlute_1.0.1 gridExtra_2.3 pathview_1.20.0
## [4] org.Hs.eg.db_3.6.0 AnnotationDbi_1.42.1 IRanges_2.14.12
## [7] S4Vectors_0.18.3 Biobase_2.40.0 BiocGenerics_0.26.0
## [10] ggplot2_3.0.0
##
## loaded via a namespace (and not attached):
## [1] fgsea_1.6.0 colorspace_1.3-2 ggridges_0.5.1
## [4] rprojroot_1.3-2 qvalue_2.12.0 XVector_0.20.0
## [7] farver_1.0 ggrepel_0.8.0 bit64_0.9-7
## [10] codetools_0.2-15 splines_3.5.1 GOSemSim_2.6.2
## [13] knitr_1.20 annotate_1.58.0 GO.db_3.6.0
## [16] png_0.1-7 pheatmap_1.0.10 graph_1.58.1
## [19] shiny_1.1.0 ggforce_0.1.3 compiler_3.5.1
## [22] httr_1.3.1 rvcheck_0.1.1 GOstats_2.46.0
## [25] backports_1.1.2 assertthat_0.2.0 Matrix_1.2-14
## [28] lazyeval_0.2.1 limma_3.36.5 later_0.7.5
## [31] tweenr_1.0.0 htmltools_0.3.6 prettyunits_1.0.2
## [34] tools_3.5.1 bindrcpp_0.2.2 igraph_1.2.2
## [37] gtable_0.2.0 glue_1.3.0 Category_2.46.0
## [40] reshape2_1.4.3 DO.db_2.9 dplyr_0.7.6
## [43] fastmatch_1.1-0 Rcpp_0.12.19 enrichplot_1.0.2
## [46] Biostrings_2.48.0 nlme_3.1-137 ggraph_1.0.2
## [49] stringr_1.3.1 mime_0.6 miniUI_0.1.1.1
## [52] clusterProfiler_3.8.1 XML_3.98-1.16 DOSE_3.6.1
## [55] zlibbioc_1.26.0 MASS_7.3-50 scales_1.0.0
## [58] promises_1.0.1 hms_0.4.2 RBGL_1.56.0
## [61] KEGGgraph_1.40.0 RColorBrewer_1.1-2 curl_3.2
## [64] yaml_2.2.0 memoise_1.1.0 UpSetR_1.3.3
## [67] biomaRt_2.36.1 ggExtra_0.8 stringi_1.2.4
## [70] RSQLite_2.1.1 genefilter_1.62.0 BiocParallel_1.14.2
## [73] matrixStats_0.54.0 rlang_0.2.2 pkgconfig_2.0.2
## [76] bitops_1.0-6 evaluate_0.11 lattice_0.20-35
## [79] purrr_0.2.5 bindr_0.1.1 labeling_0.3
## [82] cowplot_0.9.3 bit_1.1-14 tidyselect_0.2.4
## [85] GSEABase_1.42.0 AnnotationForge_1.22.2 ggsci_2.9
## [88] plyr_1.8.4 magrittr_1.5 R6_2.3.0
## [91] DBI_1.0.0 mgcv_1.8-24 pillar_1.3.0
## [94] withr_2.1.2 units_0.6-1 survival_2.42-6
## [97] KEGGREST_1.20.2 RCurl_1.95-4.11 tibble_1.4.2
## [100] crayon_1.3.4 rmarkdown_1.10 viridis_0.5.1
## [103] progress_1.2.0 grid_3.5.1 sva_3.28.0
## [106] data.table_1.11.8 blob_1.1.1 Rgraphviz_2.24.0
## [109] digest_0.6.17 xtable_1.8-3 httpuv_1.4.5
## [112] tidyr_0.8.1 munsell_0.5.0 viridisLite_0.3.0References
Aravind Subramanian, Vamsi K. Moothaa, Pablo Tamayo, and Jill P. Mesirov. 2005. “Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles.” http://www.pnas.org/content/102/43/15545.full.
Consortium., The Gene Ontology. 2014. “Gene Ontology Consortium: going forward.” https://academic.oup.com/nar/article/43/D1/D1049/2439067.
Guangchuang Yu, Yanyan Han, Li-Gen Wang, and Qing-Yu He. 2012. “clusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters.” http://online.liebertpub.com/doi/abs/10.1089/omi.2011.0118.
Hiroko Koike-Yusa, E-Pien Tan, Yilong Li. 2014. “Genome-wide recessive genetic screening in mammalian cells with a lentiviral CRISPR-guide RNA library.” http://science.sciencemag.org/content/343/6166/80.long.
Huang da W, Lempicki RA., Sherman BT. 2009. “Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists.” https://academic.oup.com/nar/article-lookup/doi/10.1093/nar/gkn923.
Laurent Gautier, Benjamin M. Bolstad, Leslie Cope. 2004. “affy—analysis of Affymetrix GeneChip data at the probe level.” https://academic.oup.com/bioinformatics/article/20/3/307/185980.
Luke A.Gilbert, BrittAdamson, Max A.Horlbeck. 2014. “Genome-Scale CRISPR-Mediated Control of Gene Repression and Activation.” https://linkinghub.elsevier.com/retrieve/pii/S0092-8674(14)01178-7.
Minoru Kanehisa, Yoko Sato, Susumu Goto. 2014. “Data, information, knowledge and principle: back to metabolism in KEGG.” https://academic.oup.com/nar/article-lookup/doi/10.1093/nar/gkt1076.
Ophir Shalem1, *, 2. 2014. “Genome-scale CRISPR-Cas9 knockout screening in human cells.” http://science.sciencemag.org/content/343/6166/84.long.
S. Falcon, R. Gentleman. 2007. “Using GOstats to test gene lists for GO term association.” https://academic.oup.com/bioinformatics/article/23/2/257/204776.
Silvana Konermann, Alexandro E. Trevino, Mark D. Brigham. 2015. “Genome-scale transcriptional activation by an engineered CRISPR-Cas9 complex.” https://www.nature.com/nature/journal/vnfv/ncurrent/full/nature14136.html.
Tim Wang, David M. Sabatini, Jenny J. Wei1. 2014. “Genetic Screens in Human Cells Using the CRISPR-Cas9 System.” http://science.sciencemag.org/content/343/6166/80.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.
Weijun Luo, Cory Brouwer. 2013. “Pathview: an R/Bioconductor package for pathway-based data integration and visualization.” https://academic.oup.com/bioinformatics/article-lookup/doi/10.1093/bioinformatics/btt285.