In bioinfo-pf-curie/cookieCrispR: Downstream analysis of Crispr screening counts data

opts_chunk$set(
  echo=input$report_echo,
  error=TRUE,
  dpi = 100
)

#Files
evalFiles <- FALSE
evalJoined <- FALSE
evalSampleplan <- FALSE
evalCounts <- FALSE
evalGenesList <- FALSE
if("DataOverview" %in% input$include){
  evalFiles <- TRUE
  if (!is.null(reactives$sampleplan) && !is.null(reactives$countsRaw)){evalJoined <- TRUE}
  if (!is.null(reactives$sampleplan)){evalSampleplan <- TRUE}
  if (!is.null(reactives$countsRaw)){evalCounts <- TRUE}
  if (!is.null(input$essential) && !is.null(input$nonessential)){
    evalGenesList <- TRUE
  }
}
evalJoined2 <- FALSE
evalSampleplan2 <- FALSE
evalCounts2 <- FALSE
evalGenesList2 <- FALSE
if (!is.null(reactives$sampleplan) && !is.null(reactives$countsRaw)){evalJoined2 <- TRUE}
if (!is.null(reactives$sampleplan)){evalSampleplan2 <- TRUE}
if (!is.null(reactives$countsRaw)){evalCounts2 <- TRUE}
if (!is.null(input$essential) && !is.null(input$nonessential)){
  evalGenesList2 <- TRUE
}

evalJoineddistribs <- TRUE
evaldistribs <- FALSE
if (TRUE %in% c(c("read_numbers","density_ridges") %in% input$include)){
  evaldistribs <- TRUE
  if(!is.null(reactives$sampleplan) && !is.null(reactives$countsRaw)){
  evalJoineddistribs <- FALSE
  }
}
## Distribs
evalreadnumbers <- FALSE
if ("read_numbers" %in% input$include){
    evalreadnumbers <- TRUE
}

evaldensity <- FALSE
evaldensityGenesList <- FALSE
evaldensityGenesList2 <- FALSE
if ("density_ridges" %in% input$include){
      evaldensity <- TRUE
      if (!is.null(reactives$sampleplan) && !is.null(reactives$countsRaw)){
      evaldensityGenesList <- TRUE
      if (!is.null(input$essential) && !is.null(input$nonessential)){
        evaldensityGenesList2 <- TRUE
      }
}
}

## Descriptive data analysis
evalJoineddescript <- TRUE
evalddescript <- FALSE
evalddescriptGenesList <- FALSE
if (TRUE %in% c(c("temporal_evolution") %in% input$include)){
  evalddescript <- TRUE
  if(!is.null(reactives$sampleplan) && !is.null(reactives$countsRaw)){
   evalJoineddescript <- FALSE
  }
  if (!is.null(input$essential) && !is.null(input$nonessential)){
  evalddescriptGenesList <- TRUE
}
}
evaltemporal_evolution <- FALSE
evaltemporal_evolutionGenesList <- FALSE
if ("temporal_evolution" %in% input$include){
    evaltemporal_evolution <- TRUE
    if (!is.null(input$essential) && !is.null(input$nonessential)){
  evaltemporal_evolutionGenesList <- TRUE
}
}

evalstatanal <- FALSE
evalJoinedtatanal <- FALSE
## Statistical Analysis
if (TRUE %in% c(c("Volcano plots") %in% input$include)){
  evalstatanal <- TRUE
  if(!is.null(reactives$sampleplan) && !is.null(reactives$countsRaw)){
   evalJoinedtatanal <- TRUE
  }
}

## SessionInfo
evalSessionInfo <- FALSE
if("SessionInfo" %in% input$include){
  evalSessionInfo <- TRUE
}

options(kableExtra.auto_format = F)

```{asis, echo = evalFiles, warning=FALSE}

Data overview

Files

```r
br()
write.table(reactives$countsRaw,file = "/tmp/Rawcounts.csv",sep = ",",quote = FALSE)
br()
xfun::embed_files(path = "/tmp/Rawcounts.csv", text = "Counts Matrix File")

br()
write.table(reactives$sampleplan,file = "/tmp/sampleplan.csv",sep = ",",quote = FALSE)
br()
xfun::embed_files(path = "/tmp/sampleplan.csv", text = "SamplePlan file")

print("<h1 style='color:red;'> Please upload counts matrix first </h1>")

print("<h1 style='color:red;'> Please upload SamplePlan first </h1>")

```{asis, echo = evalFiles, warning=FALSE}

Samples Selected

The following samples were selected in this analysis :

```{asis, echo = evaldistribs, warning=FALSE}

Distributions

<!-- ```r -->
<!-- if(evalJoineddistribs == TRUE){ -->
<!-- print("Please upload both SamplePlan and counts matrix first") -->
<!-- } -->
<!-- ``` -->


```{asis, echo = evalreadnumbers}
### Read Numbers

if(!(is.null(input$timepoints_order))){
        counts <- full_join(reactives$counts, reactives$sampleplan) %>%
          mutate(Timepoint = factor(.data$Timepoint, levels = input$timepoints_order))
      } else {
        print(reactives$sampleplan$Timepoint)
        counts <- full_join(reactives$counts, reactives$sampleplan) %>%
          mutate(Timepoint = factor(.data$Timepoint, levels = unique(.data$Timepoint)))
      }
      counts <- as.data.frame(counts) %>%
        group_by(.data$Sample_ID, .data$Replicate, .data$Cell_line, .data$Timepoint) %>%
        summarise(total = sum(.data$count)) %>% 
        as.data.frame() %>%
        ggplot(aes(x = .data$Timepoint, y = .data$total)) +
        geom_col(position = position_dodge()) + facet_wrap(vars(.data$Replicate, .data$Cell_line), nrow = 1) +
        labs(title = "Number of reads per sample", xlab ="Timepoint", ylab ="Total counts")

plot(counts)

```{asis, echo = evaldensity}

Density ridges

All guides

```r
if(evaldensityGenesList == TRUE){
 counts <- reactives$joined
     plot(counts %>%
        ggplot(aes(x = .data$log_cpm, y = .data$Timepoint)) +
        geom_density_ridges(alpha = 0.6, show.legend = FALSE, fill = "gray50") +
        facet_wrap(vars(.data$Cell_line, .data$Replicate), ncol = 1, strip.position = "right") +
        labs(title = "Distribution of normalzed log-cpm by sample", subtitle = "All guides"))
}

```{asis, echo = evaldensity}

Per genes type density ridges

```r
if(evaldensityGenesList2 == TRUE){
ess_genes <- ess_genes()
non_ess_genes <- non_ess_genes()

#counts <- filter(counts, .data$Gene %in% ess_genes$X)
counts <- filter(reactives$joined, .data$gene %in% ess_genes$X)
plot1 <- ggplot(counts, aes(x = .data$log_cpm, y = .data$Timepoint)) + 
  geom_density_ridges(alpha = 0.6, show.legend = FALSE, fill = "gray50") +
  facet_wrap(vars(.data$Cell_line, .data$Replicate), ncol = 1, strip.position = "right") +
  scale_fill_viridis_c() +
  labs(title = "Distribution of normalised log-cpms", subtitle = "Essential genes")

counts <- reactives$joined %>%
        filter(.data$gene %in% c(as.character(non_ess_genes$X),"Non-Targeting","negative_control"))

  plot2 <- ggplot(counts, aes(x = .data$log_cpm, y = .data$Timepoint)) +
  geom_density_ridges(alpha = 0.6, show.legend = FALSE, fill = "gray50") +
  facet_wrap(vars(.data$Cell_line, .data$Replicate), ncol = 1, strip.position = "right") +
  scale_fill_viridis_c() +
  labs(title = "Distribution of normalised log-cpms", subtitle = "Non Essential genes")

grid.arrange(plot1,plot2,ncol = 2)

}

if (evaldensityGenesList2 == FALSE){
 print("Provide both essential and non essential genes lists first") 
}

```{asis,echo = evalddescript}

Descriptive data analysis

```{asis, echo = evaltemporal_evolution, warning=FALSE}
### Temporal evolution
##### All guides

if(is.null(diff_boxes$diff_box_all) | is.null(diff_boxes$diff_box_ess)){
  print("Please compute diff boxes ont the temporal evolution tab first")
} else {
plot(diff_boxes$diff_box_all)
}

```{asis, echo = evaltemporal_evolution, warning=FALSE}

Essential genes

```r
if(evaltemporal_evolutionGenesList == FALSE){
print("Please provide essential and non essential genes lists first")
}

if(evaltemporal_evolutionGenesList == TRUE){
if(is.null(diff_boxes$diff_box_all) | is.null(diff_boxes$diff_box_ess)){
  print("Please compute diff boxes ont the temporal evolution tab first")
} else {
plot(diff_boxes$diff_box_ess)
}
}

```{asis,echo = evalstatanal}

Statistical data analysis

```{asis, echo = evalJoinedtatanal}
# Volcano plot

if(!is.null(DEA$concatenated$results) & !is.null(DEA$reactives$selectedcomp) & !is.null(DEA$reactives$selectedFC) & !is.null(DEA$reactives$selectedPvalT)) {

volcanolist <- list()

FCT <- as.numeric(DEA$reactives$selectedFC)
PVT <- as.numeric(DEA$reactives$selectedPvalT)
res <- as.data.frame(DEA$concatenated$results[[DEA$reactives$selectedcomp]])

ggplot <- ggplot(res, aes(x = estimate, y = -log10(p.value))) + 
           ggtitle(DEA$reactives$selectedcomp) + 
           scale_color_gradient(low = "lightgray", high = "navy") + 
           geom_point(data = res, 
                  color = "grey", alpha = 0.5) + 
           geom_point(data = subset(res, estimate > FCT), 
                  color = "red", alpha = 0.5) + 
           geom_point(data = subset(res, estimate < -FCT), 
                  color = "blue", alpha = 0.5) + 
           geom_point(data = subset(res, adj_p.value < PVT),
                  color = "green", alpha = 0.5) +
           geom_hline(yintercept = -log10(max(subset(res, adj_p.value < PVT)$p.value)), linetype = "dashed") +
          geom_vline(xintercept = c(-FCT,FCT), linetype = "dashed") + 
           theme_linedraw() + 
           theme(panel.grid = element_blank()) + 
           xlab("Fold change (log2)") +
           ylab("-log10(P-Value)") +
  geom_point(data = subset(res,sgRNA %in% DEA$reactives$GeneVolcano),
                 color = "purple", alpha = 0.6) +
      ggrepel::geom_text_repel(
        data = subset(res,sgRNA %in% DEA$reactives$GeneVolcano),
        aes(label = sgRNA),
        size = 5,
        force = 2,
        box.padding = unit(0.35, "lines"),
        point.padding = unit(0.3, "lines")
      )
plot(ggplot)
}

if(!is.null(DEA$concatenated$results) & !is.null(DEA$reactives$selectedcomp) & !is.null(DEA$reactives$selectedFC) & !is.null(DEA$reactives$selectedPvalT)) {
} else {
  print("First draw a volcano plot in the Explore DEA results tab")
}

```{asis,echo=evalSessionInfo}

Session Information

```r
sessionInfo()

library(shiny)
footertemplate <- function(){
  tags$div(
    class = "footer",
    style = "text-align:center",
    tags$div(
      class = "foot-inner",
      list(
        hr(),
        "This report was generated with", tags$a(href="http://bioconductor.org/packages/COOKIE CRISP'R/", "COOKIE CRISP'R"), br(),
        "COOKIE CRISP'R is a project developed by Nicolas Servant, Pierre Gestraud, Aurelien Bore and Clement Benoit from the bioinformatic and crispr'it platforms of the Institut Curie : ",
        br(),
        tags$a(href="https://science.curie.fr/plateformes/criblage-genetique-crispr-crisprit/","CRISPRIT platform"),
        br(),
        tags$a(href="https://science.curie.fr/recherche/developpement-cancer-genetique-epigenetique/dynamique-de-linformation-genetique/equipe-vallot/","Bioinformatic Platform"),br())
      )
    )
}