## ----setup-------------------------------------------------------------------- knitr::opts_chunk$set(tidy=FALSE, cache=TRUE, dev="png", message=FALSE, error=FALSE, warning=TRUE) ## ----library, eval=TRUE------------------------------------------------------- if(!"MAGeCKFlute" %in% installed.packages()) BiocManager::install("MAGeCKFlute") if(!"clusterProfiler" %in% installed.packages()) BiocManager::install("clusterProfiler") if(!"ggplot2" %in% installed.packages()) BiocManager::install("ggplot2") library(MAGeCKFlute) library(clusterProfiler) library(ggplot2) ## ----fluterra, eval=FALSE----------------------------------------------------- # ## path to the gene summary file (required) # file1 = file.path(system.file("extdata", package = "MAGeCKFlute"), # "testdata/rra.gene_summary.txt") # ## path to the sgRNA summary file (optional) # file2 = file.path(system.file("extdata", package = "MAGeCKFlute"), # "testdata/rra.sgrna_summary.txt") # # Run FluteRRA with only gene summary file # FluteRRA(file1, proj="Pmel1", organism="mmu", outdir = "./") # # # Run FluteRRA with both gene summary file and sgRNA summary file # FluteRRA(file1, file2, proj="Pmel1", organism="mmu", outdir = "./") ## ----incorporateDepmap, eval=FALSE-------------------------------------------- # FluteRRA(file1, proj="Pmel1", organism="mmu", incorporateDepmap = TRUE, outdir = "./") ## ----omitEssential, eval=FALSE------------------------------------------------ # FluteRRA(gdata, proj="Pmel1", organism="mmu", omitEssential = TRUE, outdir = "./") ## ----flutemle, eval=FALSE----------------------------------------------------- # file3 = file.path(system.file("extdata", package = "MAGeCKFlute"), # "testdata/mle.gene_summary.txt") # FluteMLE(file3, treatname="Pmel1", ctrlname="Pmel1_Ctrl", proj="Pmel1", organism="mmu") ## ----incorporateDepmap2, eval=FALSE------------------------------------------- # ## Take Depmap screen as control # FluteMLE(gdata, treatname="Pmel1_Ctrl", ctrlname="Depmap", proj="Pmel1", organism="mmu", incorporateDepmap = TRUE) ## ----omitEssential2, eval=FALSE----------------------------------------------- # FluteMLE(gdata, treatname="Pmel1", ctrlname="Pmel1_Ctrl", proj="Pmel1", organism="mmu", omitEssential = TRUE) ## ----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")) gg = gg[order(gg$Label, gg$variable), ] 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) ## ----mmu2hsa------------------------------------------------------------------ gdata$HumanGene = TransGeneID(gdata$id, fromType = "symbol", toType = "symbol", fromOrg = "mmu", toOrg = "hsa") sdata$HumanGene = TransGeneID(sdata$Gene, fromType = "symbol", toType = "symbol", fromOrg = "mmu", toOrg = "hsa") ## ----Depmap, eval=FALSE------------------------------------------------------- # ## Remove missing or duplicate human genes # idx = duplicated(gdata$HumanGene)|is.na(gdata$HumanGene) # gdata = gdata[!idx, ] # depmap_similarity = ResembleDepmap(gdata, symbol = "HumanGene", score = "Score") # head(depmap_similarity) ## ----omitessential, eval=FALSE------------------------------------------------ # gdata = OmitCommonEssential(gdata, symbol = "HumanGene") # sdata = OmitCommonEssential(sdata, symbol = "HumanGene") ## ----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", top = 5, auto_cut_y = TRUE, groups = c("top", "bottom"), ylab = "Log2FC", toplabels = c("Pbrm1", "Arid2", "Brd7")) ## ----rankrra2, fig.height=4, fig.width=3-------------------------------------- geneList= gdata$Score names(geneList) = gdata$id p2 = RankView(geneList, top = 5, bottom = 10) print(p2) RankView(geneList, top = 0, bottom = 0, genelist = c("Pbrm1", "Arid2", "Brd7")) ## ----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_pos = EnrichAnalyzer(geneList = geneList[geneList>0.5], method = "HGT", type = "KEGG") enrich_neg = EnrichAnalyzer(geneList = geneList[geneList< -0.5], method = "HGT", type = "KEGG") ## ----enrichview--------------------------------------------------------------- EnrichedView(enrich_pos, mode = 1, top = 5, bottom = 0) EnrichedView(enrich_pos, mode = 2, top = 5, bottom = 0) EnrichedView(enrich_neg, mode = 2, top = 0, bottom = 5) ## ----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) ## ----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) ## ----NormalizeBeta------------------------------------------------------------ ctrlname = "Pmel1_Ctrl" treatname = "Pmel1" gdata$HumanGene = TransGeneID(gdata$Gene, fromType = "symbol", toType = "symbol", fromOrg = "mmu", toOrg = "hsa") gdata_cc = NormalizeBeta(gdata, id = "HumanGene", 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=4, fig.width=5------------------------------------ gdata_cc$Control = rowMeans(gdata_cc[,ctrlname, drop = FALSE]) gdata_cc$Treatment = rowMeans(gdata_cc[,treatname, drop = FALSE]) p1 = ScatterView(gdata_cc, "Control", "Treatment", label = "Gene", auto_cut_diag = TRUE, display_cut = TRUE, groups = c("top", "bottom"), toplabels = c("Pbrm1", "Brd7", "Arid2", "Jak1", "Stat1", "B2m")) print(p1) ## ----rank, fig.height=4, fig.width=3------------------------------------------ gdata_cc$Diff = gdata_cc$Treatment - gdata_cc$Control gdata_cc$Rank = rank(gdata_cc$Diff) p1 = ScatterView(gdata_cc, x = "Rank", y = "Diff", label = "Gene", top = 5, auto_cut_y = TRUE, groups = c("top", "bottom")) print(p1) # Or rankdata = gdata_cc$Treatment - gdata_cc$Control names(rankdata) = gdata_cc$Gene RankView(rankdata) ## ----Square, fig.height=4, fig.width=5---------------------------------------- p1 = ScatterView(gdata_cc, x="Pmel1_Ctrl", y="Pmel1", label = "Gene", model = "ninesquare", top = 5, display_cut = TRUE, force = 2) print(p1) ## ----Square2, fig.height=4, fig.width=5--------------------------------------- p1 = ScatterView(gdata_cc, x="Pmel1_Ctrl", y="Pmel1", label = "Gene", model = "ninesquare", top = 5, display_cut = TRUE, x_cut = c(-1,1), y_cut = c(-1,1)) print(p1) ## ----fig.height=4, fig.width=5------------------------------------------------ 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=6---------------------------------- # 9-square groups Square9 = p1$data idx=Square9$group=="topcenter" geneList = Square9$Diff[idx] names(geneList) = Square9$Gene[idx] universe = Square9$Gene # Enrichment analysis kegg1 = EnrichAnalyzer(geneList = geneList, universe = universe) EnrichedView(kegg1, top = 6, bottom = 0) ## ----pathview2, eval=TRUE----------------------------------------------------- gdata_cc$Entrez = TransGeneID(gdata_cc$Gene, "symbol", "entrez", organism = "mmu") genedata = gdata_cc[, c("Entrez", "Control","Treatment")] arrangePathview(genedata, pathways = "mmu04630", organism = "mmu", sub = NULL) ## ----sessionInfo-------------------------------------------------------------- sessionInfo()