In vallotlab/ChromSCape: Analysis of single-cell epigenomics datasets with a Shiny App

options(knitr.duplicate.label = "allow")

MSIG.classes <- c("c1_positional","c2_curated","c3_motif","c4_computational",
          "c5_GO","c6_oncogenic","c7_immunologic","hallmark")[c(2,5,6,7,8)]

QC & Reduction

Column {data-width=250}

Cells

flexdashboard::valueBox(value = ncol(datamatrix), caption = "Number of cells", icon ="fa-circle")

Cells after filterings

if("filter" %in% run) flexdashboard::valueBox(value = ncol(scExp), caption = "Number of cells after filtering", icon ="fa-circle-notch")

Percentage of features passing filtering

if("filter" %in% run) flexdashboard::gauge(100*round(nrow(scExp)/ nrow(datamatrix),3), min = 0, max = 100, 
          symbol = '%', sectors = flexdashboard::gaugeSectors(
    success = c(20, 100), warning = c(5, 19.99), danger = c(0, 4.999)
  ))

Raw Cell count distribution

df = data.frame(coverage = sort(unname(Matrix::colSums(datamatrix)))) 
ggplot2::ggplot(df, ggplot2::aes(x = coverage)) + 
      ggplot2::geom_histogram(color="black", fill="steelblue", bins = 75) +
      ggplot2::labs(x="Log10(Reads per cell)", y = "nCells")  + 
      ggplot2::theme(panel.grid.major = ggplot2::element_blank(), panel.grid.minor= ggplot2::element_blank(), 
            panel.background = ggplot2::element_blank(), axis.line = ggplot2::element_line(colour="black"),
            panel.border = ggplot2::element_rect(colour="black", fill=NA)) +
      ggplot2::scale_x_log10()

Column {data-width=400}

UMAP - sample

if(any(c("cluster","DA") %in% run)) plot_reduced_dim_scExp(scExp_cf, color_by = "sample_id", reduced_dim = "UMAP")

UMAP - total counts

if(any(c("cluster","DA") %in% run)) plot_reduced_dim_scExp(scExp_cf, color_by = "total_counts", reduced_dim = "UMAP")

Column {data-width=400}

Most contributing features

if(any(c("filter") %in% run)) plot_most_contributing_features(scExp, component = colnames(SingleCellExperiment::reducedDim(scExp,"PCA"))[1])

Most contributing chromosomes

if(any(c("filter") %in% run)) plot_pie_most_contributing_chr(scExp, component = colnames(SingleCellExperiment::reducedDim(scExp,"PCA"))[1])

Correlation Clustering

Column {data-width=250}

Number of cluster

if(any(c("cluster", "DA") %in% run)) flexdashboard::valueBox(value = length(unique(scExp_cf$cell_cluster)), caption = "Number of clusters", icon ="fa-sitemap")

Percentage of well correlated cells

if(any(c("cluster", "DA") %in% run)) 
  flexdashboard::gauge(100*round(ncol(scExp_cf) / ncol(scExp),1), min = 0, max = 100, 
        symbol = '%', sectors = flexdashboard::gaugeSectors(
          success = c(80, 100), warning = c(50, 79.99), danger = c(0, 49.999)
        ))

Cell affectation

if(any(c("cluster", "DA") %in% run)) num_cell_in_cluster_scExp(scExp_cf)

Column {data-width=450}

Correlation

if(any(c("cluster", "DA") %in% run)){
  if("hc_cor" %in% names(scExp_cf@metadata)){
    plot_heatmap_scExp(scExp_cf)
  }
}

Consensus scores

if(any(c("consensus") %in% run)){
  if("icl" %in% names(scExp_cf@metadata)){
    plot_cluster_consensus_scExp(scExp_cf)
  }
}

Column {data-width=350}

UMAP - cell cluster

if(any(c("cluster", "DA") %in% run)){
  if("cell_cluster" %in% colnames(SummarizedExperiment::colData(scExp_cf))){
    plot_reduced_dim_scExp(scExp_cf, color_by = "cell_cluster", reduced_dim = "UMAP")
  }
}

Sample intra-correlation

if(any(c("cluster", "DA") %in% run)){
  if("Cor" %in% SingleCellExperiment::reducedDimNames(scExp_cf)){
    plot_intra_correlation_scExp(scExp_cf, by = "sample_id")
  }
}

DA & GSA

Column

Number of differential regions

if(any(c("DA") %in% run)){
  if( any(grepl("qval",colnames(SingleCellExperiment::rowData(scExp_cf)))) ){
    flexdashboard::valueBox(value = sum(summary_DA(scExp_cf)["differential",]), caption = "Number of differential regons", icon ="fa-cube")
  }
}

Number unique gene sets

if(any(c("DA") %in% run)){
  if(any(grepl("qval",colnames(SingleCellExperiment::rowData(scExp_cf)))) & "enr" %in% names(scExp_cf@metadata)){
      groups = gsub(".*\\.","", grep("qval",colnames(SingleCellExperiment::rowData(scExp_cf)), value = TRUE))
    all_paths <- unique(as.character(unlist(sapply(1:length(groups), 
                                                   function(i) {
                                                     all = c()
                                                     if(length(scExp_cf@metadata$enr$Overexpressed)>=i) 
                                                       all = c(all, scExp_cf@metadata$enr$Overexpressed[[i]]$Gene.Set)
                                                     if(length(scExp_cf@metadata$enr$Underexpressed)>=i) 
                                                       all = c(all, scExp_cf@metadata$enr$Underexpressed[[i]]$Gene.Set)
                                                     if(length(scExp_cf@metadata$enr$Both)>=i) 
                                                       all = c(all, scExp_cf@metadata$enr$Both[[i]]$Gene.Set)
                                                     return(all)
                                                   }))))

    flexdashboard::valueBox(value = length(all_paths), caption = "Number of unique pathways", icon ="fa-cubes")
  }
}

Cell affectation

if(any(c("cluster", "DA") %in% run)){
  if("cell_cluster" %in% colnames(SummarizedExperiment::colData(scExp_cf))){
    num_cell_in_cluster_scExp(scExp_cf)
  }
}

Column

Differential Volcano plots

out = NULL
if(any(c("DA") %in% run)){
  if(any(grepl("qval",colnames(SingleCellExperiment::rowData(scExp_cf)))) & "enr" %in% names(scExp_cf@metadata)){
    diff_list = list()
    if(!is.null(scExp_cf)){
        groups = gsub(".*\\.","", grep("qval",colnames(SingleCellExperiment::rowData(scExp_cf)), value = TRUE))


        for(i in groups){

          diff = as.data.frame(
              SummarizedExperiment::rowData(scExp_cf)[,-grep("Rank|pval|ID|Count",colnames(SummarizedExperiment::rowData(scExp_cf)))]
          )
          if(nrow(diff)>0){

            diff = diff[,c("Gene","distanceToTSS", colnames(diff)[grep(i,colnames(diff))] )]
            diff = diff[which(!is.infinite(diff[,grep("logFC.",colnames(diff))])),]
            diff <- diff[order(diff[,paste0("qval.",i)]),]
            diff[,paste0("qval.",i)] = round(diff[,paste0("qval.",i)],8)
            diff[,paste0("logFC.",i)] = round(diff[,paste0("logFC.",i)],4)
            diff$Gene = substr(diff$Gene,1,20)
            diff_list[[i]] = diff
          }
          # collapse together all lines with newline separator
        }

    }

    out_volcano <- lapply(seq_along(groups), function(i) {

      a1 <- knitr::knit_expand(text = sprintf("### %s\n", groups[i])) 
      a2 <- knitr::knit_expand(text = "\n```r") # start r chunk
      a3 <- knitr::knit_expand(text = sprintf("\nplot_differential_volcano_scExp(scExp_cf, group = groups[%d])", i)) 
      a4 <- knitr::knit_expand(text = "\n```\n") 

      paste(a1, a2, a3, a4, collapse = '\n') # collapse together all lines with newline separator

    })

    out_table <- lapply(seq_along(groups), function(i) {

      a1 <- knitr::knit_expand(text = sprintf("### %s\n", groups[i])) 
      a2 <- knitr::knit_expand(text = "\n```r") # start r chunk
      a3 <- knitr::knit_expand(text = sprintf("\nDT::datatable(diff_list[[%d]][1:min(20,nrow(diff_list[[%d]])),], options = list(pageLength = 5, bPaginate = FALSE,  dom = 'tip'), class = 'display', rownames = FALSE)", i,i)) 
      a4 <- knitr::knit_expand(text = "\n```\n") 

      paste(a1, a2, a3, a4, collapse = '\n')
    })

    out = list()
    for(i in seq_along(groups)){
      out = c(out, out_volcano[[i]], out_table[[i]])  
    }
  }
}

r if(any(c("DA") %in% run)) paste(knitr::knit(text = paste(out, collapse = '\n')))

Column

options(DT.options = list(scrollY="300px",scrollX="300px", pageLength = 10, autoWidth = TRUE))

out = NULL
gsa_over = gsa_under = list()
if(any(c("GSA") %in% run)){
  if(!is.null(scExp_cf@metadata$enr)){

    for(i in seq_along(scExp_cf@metadata$diff$groups)){
      if(length(scExp_cf@metadata$enr$Overexpressed)>=i){
        if(!is.null(scExp_cf@metadata$enr$Overexpressed[[i]])){
          over_tab = table_enriched_genes_scExp(scExp_cf, set = 'Overexpressed',  group = paste0('C',i), enr_class_sel = MSIG.classes)
          colnames(over_tab)[3] = "N_genes"
          colnames(over_tab)[6] = "Genes"
          over_tab$p.value =   round(over_tab$p.value,7)
          over_tab$adj.p.value =   round(over_tab$adj.p.value,7)
          over_tab$Gene_set = ifelse(nchar(over_tab$Gene_set)>25, paste0(substr(over_tab$Gene_set,1,25),"..."), over_tab$Gene_set)
          over_tab$Genes = ifelse(nchar(over_tab$Genes)>20, paste0(substr(over_tab$Genes,1,20),"..."), over_tab$Genes)
          gsa_over[[i]] = over_tab
        }
      }
      if(length(scExp_cf@metadata$enr$Underexpressed)>=i){
        if(!is.null(scExp_cf@metadata$enr$Underexpressed[[i]])){
          under_tab = table_enriched_genes_scExp(scExp_cf, set = 'Underexpressed',  group = paste0('C',i), enr_class_sel = MSIG.classes)
          colnames(under_tab)[3] = "N_genes"
          colnames(under_tab)[6] = "Genes"
          under_tab$p.value =   round(under_tab$p.value,7)
          under_tab$adj.p.value =   round(under_tab$adj.p.value,7)
          under_tab$Gene_set = ifelse(nchar(under_tab$Gene_set)>25, paste0(substr(under_tab$Gene_set,1,25),"..."), under_tab$Gene_set)
          under_tab$Genes = ifelse(nchar(under_tab$Genes)>20, paste0(substr(under_tab$Genes,1,20),"..."), under_tab$Genes)
          gsa_under[[i]] = under_tab
        }
      }
    }
  }

  out_gsa_table_up <- lapply(seq_along(scExp_cf@metadata$diff$groups), function(i) {
    if(!is.null(scExp_cf@metadata$enr)){
      if(length(scExp_cf@metadata$enr$Overexpressed)>=i){
        if(!is.null(scExp_cf@metadata$enr$Overexpressed[[i]])){
          if(nrow(gsa_over[[i]])>0){
            a1 <- knitr::knit_expand(text = sprintf("### Gene sets associated with enriched mark - %s \n", scExp_cf@metadata$diff$groups[i])) 
            a2 <- knitr::knit_expand(text = "\n```r") # start r chunk
            a3 <- knitr::knit_expand(text = sprintf("\nDT::datatable(gsa_over[[%d]], options = list(pageLength = 10,  dom = 'tip', bPaginate = FALSE), class = 'display', rownames = FALSE)", i))
            a4 <- knitr::knit_expand(text = "\n```\n") 

            paste(a1, a2, a3, a4, collapse = '\n')
          }

        }

      }
    }
  })

  out_gsa_table_dn <- lapply(seq_along(scExp_cf@metadata$diff$groups), function(i) {
    if(!is.null(scExp_cf@metadata$enr)){
      if(length(scExp_cf@metadata$enr$Underexpressed)>=i){
        if(!is.null(scExp_cf@metadata$enr$Underexpressed[[i]])){
          if(nrow(gsa_under[[i]])>0){
            a1 <- knitr::knit_expand(text = sprintf("### Gene sets associated with depleted mark - %s \n", scExp_cf@metadata$diff$groups[i])) 
            a2 <- knitr::knit_expand(text = "\n```r") # start r chunk
            a3 <- knitr::knit_expand(text = sprintf("\nDT::datatable(gsa_over[[%d]], options = list(pageLength = 10,  dom = 'tip', bPaginate = FALSE), class = 'display', rownames = FALSE)", i))
            a4 <- knitr::knit_expand(text = "\n```\n") 

            paste(a1, a2, a3, a4, collapse = '\n')
          }
        }
      }
    }
  })

  out = list()
  for(i in seq_along(scExp_cf@metadata$diff$groups)){

    out = c(out, out_gsa_table_up[[i]], out_gsa_table_dn[[i]])  
  }

}

r if(any(c("GSA") %in% run)) paste(knitr::knit(text = paste(out, collapse = '\n')))

Coverage

Column

---

out_coverage = NULL
if(any(c("coverage") %in% run)){
    if(length(coverages)>0){
          plot_list = coverage_chr = coverage_start = coverage_end = list()
          coverage_color = setNames(unique(scExp_cf$cell_cluster_color), unique(scExp_cf$cell_cluster))

          for(i in seq_along(genes_to_plot)){
            Gene = genes_to_plot[i]
            eval(parse(text = paste0("data(", ref_genome, ".GeneTSS)")))
            gene_annot = eval(parse(text = paste0(ref_genome, ".GeneTSS")))
            if(grepl(":", Gene) ){
              coverage_chr[[i]] = gsub(":.*","",Gene)
              coverage_start[[i]] = as.numeric(gsub(".*:","",gsub("-.*","",Gene)))
              coverage_end[[i]] = as.numeric(gsub(".*-","",Gene))
            } else{
              if(Gene %in% gene_annot$Gene){
                coverage_chr[[i]] = gene_annot$chr[which(gene_annot$Gene == Gene)]
                coverage_start[[i]] = gene_annot$start[which(gene_annot$Gene == Gene)] - 25000
                coverage_end[[i]] = gene_annot$start[which(gene_annot$Gene == Gene)] + 25000
              }
            }
            if(length(coverage_chr)>=i && !is.null(coverage_chr[[i]])){
              plot_coverage_BigWig(coverages, coverage_color, chrom = coverage_chr[[i]], start = coverage_start[[i]], end = coverage_end[[i]], ref = ref_genome)
              plot_list[[i]] <- grDevices::recordPlot()
              plot.new() ## clean up device
            }
          }

          out_coverage <- lapply(seq_along(genes_to_plot), function(i) {

            if(length(coverage_chr)>=i && !is.null(coverage_chr[[i]])){
              a1 <- knitr::knit_expand(text = sprintf("### Coverage - %s \n", genes_to_plot[i]))
              a2 <- knitr::knit_expand(text = "\n```r") # start r chunk
              a3 <- knitr::knit_expand(text = sprintf("\nprint(plot_list[[%d]])", i))
              a4 <- knitr::knit_expand(text = "\n```\n")
              paste(a1, a2, a3, a4, collapse = '\n')
            }
          })
    }
}

r if(any(c("coverage") %in% run) & length(coverages)>0) paste(knitr::knit(text = paste(out_coverage, collapse = '\n')))

Copy Number

Column {data-width=250}

Percentage of copy number gains in non-controls

if(any(c("CNA") %in% run)){
  if(!is.null(control_samples_CNA)){
    gain_or_loss_cytoBand = reducedDim(scExp, "gainOrLoss_cytoBand")
    gain_or_loss_cytoBand = gain_or_loss_cytoBand[which(!scExp$sample_id %in% control_samples_CNA),]
    percent_gol = round(100 * sum(abs(gain_or_loss_cytoBand[gain_or_loss_cytoBand==1])) / (ncol(gain_or_loss_cytoBand) * nrow(gain_or_loss_cytoBand)),3)
    flexdashboard::valueBox(value = percent_gol, caption = "% of copy number gains in non-controls", icon ="fa-chevron-circle-up",  color = "danger")
  }
}

Percentage of copy number loss in non-controls

if(any(c("CNA") %in% run)){
  if(!is.null(control_samples_CNA)){
    gain_or_loss_cytoBand = reducedDim(scExp, "gainOrLoss_cytoBand")
    gain_or_loss_cytoBand = gain_or_loss_cytoBand[which(!scExp$sample_id %in% control_samples_CNA),]
    percent_gol = round(100 * sum(abs(gain_or_loss_cytoBand[gain_or_loss_cytoBand==-1])) / (ncol(gain_or_loss_cytoBand) * nrow(gain_or_loss_cytoBand)),3)
    flexdashboard::valueBox(value = percent_gol, caption = "% of copy number loss in non-controls", icon ="fa-chevron-circle-down", color = "success")
  }
}

Column {data-width=400}

Gain or Loss barplots of top 20 most variable cytoBands

if(any(c("CNA") %in% run)){
  if(!is.null(control_samples_CNA)){
    plot_gain_or_loss_barplots(scExp, cells = scExp$cell_id[which(!scExp$sample_id %in% control_samples_CNA)])
  }
}

Column {data-width=400}

UMAPs of top 1 most variable cytoBands

if(any(c("CNA") %in% run)){
    top_variable_cyto = get_most_variable_cyto(scExp, top = 4)
    plot_reduced_dim_scExp_CNA(scExp, top_variable_cyto$cytoBand[1])
}

UMAPs of top 2 most variable cytoBands

if(any(c("CNA") %in% run)){
  if(!is.null(control_samples_CNA)){
    plot_reduced_dim_scExp_CNA(scExp, top_variable_cyto$cytoBand[2])
  }
}

vallotlab/ChromSCape documentation built on Oct. 15, 2023, 1:47 p.m.