Del17p13 signature

Differentially expressed genes

library(DESeq2)
library(tidyverse)
library(ggsci)
library(matrixStats)
library(piano)
library(reshape2)
library(genefilter)
library(Biobase)
library(ComplexHeatmap)
library(ggplot2)
library(gtable)
library(grid)
library(circlize)
library(gridExtra)
library(ggpubr)
library(RColorBrewer)
library(here)

data_dir <- here("data")
output_dir <- here("output")
figure_dir <- here("output/figures")

#dds data set. gene expression data + patmetadata
load(paste0(data_dir, "/ddsrnaCLL_150218.RData"))

variant <- "del17p13"
#filter for patients without NA in variant
ddsCLL <- ddsCLL[, !is.na(colData(ddsCLL)[,variant])]

#differentially expressed genes between del17p13 groups (see differential expression.html)
diff_all <- read.csv(file=paste0(output_dir, "/diff_genes/", variant, "_diffGenes.csv"))


rownames(diff_all) <- diff_all$X
diff_all <- diff_all[which(diff_all$padj < 0.01 ),-1]
diff <- diff_all[which( abs(diff_all$stat) > 4),]


mutStatus <- data.frame(colData(ddsCLL)) %>% arrange(del17p13)
colnames(ddsCLL) <-colData(ddsCLL)$PatID
ddsCLL <- ddsCLL[, mutStatus$PatID]

#expression data
ddsCLL <- estimateSizeFactors(ddsCLL)
RNAnorm <- varianceStabilizingTransformation(ddsCLL, blind = T)

#filter for sign. genes in variant
exprMat <- assay(RNAnorm)
exprVariant <- exprMat[rownames(diff),]
colnames(exprVariant) <- colData(ddsCLL)$PatID
exprVariant.new <- log2(exprVariant)
exprVariant.new <- t(scale(t(exprVariant.new)))
exprVariant.new[exprVariant.new > 4] <- 4
exprVariant.new[exprVariant.new < -4] <- -4
rownames(exprVariant.new) <- rowData(RNAnorm[rownames(diff),])$symbol

#colors
colors = colorRamp2(c(-4,-2,0,2,4), c("#2166ac","#4393c3", "#f7f7f7", "#d6604d","#b2182b"))
annocol <- get_palette("jco", 10)
annocolor <- list(del17p13 = c("1" = annocol[8], "0" = annocol[9]))
rowcolors <-colorRampPalette(brewer.pal(5, "Set1"))(5)
rowcolors[6] <- "white"



feature <- as.data.frame(colData(ddsCLL)[,c(variant)])
colnames(feature) <- c(variant)  

ha_col <- HeatmapAnnotation(df = feature, col = annocolor, annotation_height = unit(c(rep(1.9, 1)), "cm"), 
                            simple_anno_size = unit(1, "cm"),
                            annotation_name_gp = gpar(fontsize = 22, fontface = "bold"),
                            annotation_legend_param = list(title_gp = gpar(fontsize = 23), 
                                                           labels_gp = gpar(fontsize = 18),  
                                                           grid_height = unit(1.2, "cm"), 
                                                           grid_width = unit(1.2, "cm")))


#Annotate top 50 genes
diff <- diff_all[which(abs(diff_all$stat) > 4),]
sub_names <- unique(diff$Symbol)
geneIDs <- which(rownames(exprVariant.new) %in% sub_names)
rownames(exprVariant.new)[-geneIDs] <- ""

h1 <- Heatmap(exprVariant.new ,                                                     
              km = 2,
              gap = unit(0.5, "cm"),
              cluster_columns = F,
              clustering_distance_rows = "pearson",
              clustering_method_rows = "ward.D2",
              column_title = paste0("Gene signature: ", variant),                      
              col = colors,
              column_title_gp = gpar(fontsize = 25, fontface = "bold"), 
              heatmap_legend_param = list(title = "expr", 
                                          title_gp = gpar(fontsize = 23), 
                                          grid_height = unit(1.5, "cm"), 
                                          grid_width = unit(1.2, "cm"), 
                                          gap = unit(2, "cm"), 
                                          labels_gp = gpar(fontsize = 18)), 
              column_dend_height = unit(1, "cm"),
              show_row_dend = FALSE, 
              show_column_names = FALSE , 
              show_row_names = TRUE, 
              row_names_gp = gpar(fontsize = 12),
              top_annotation = ha_col)


#svg(filename=paste0(figure_dir, "/", variant, "_gene_expr.svg"), width=30, height=45)
#pdf(file=paste0(figure_dir, "/", variant, "_gene_expr.pdf"), width=22, height=25)

draw(h1 ) 

#dev.off()

saveRDS(h1, file = paste0(output_dir, "/figures/r_objects/del17p13/del17p13_heatmap.rds"))

Sample and gene specific expression - top genes

#function to create stripchart plots for specific genes
gene_count <- function(gene_nam){
  geneEnsID <- rownames(ddsCLL)[which(rowData(ddsCLL)$symbol %in% gene_nam)]
  geneNum <- exprMat[geneEnsID,]
  mutPat <- as.data.frame(colData(ddsCLL)[, c(variant)])
  colnames(mutPat) <- c("genotype")
  geneDat <- cbind(mutPat, geneNum)
  colnames(geneDat) <- c("genotype", "counts")
  
  p <- ggstripchart(geneDat, x = "genotype", y = "counts",
          color = "genotype",
          size = 3,
          palette = "jco",
          add = "mean_sd",
          add.params = list(size = 1.2),
          title = paste(gene_nam),
          font.x = 20, font.y = 20, font.legend = 20, 
          ylab = "normalized counts") + font("xy.text", size = 20) + font("title", size = 20, face = "bold")
  saveRDS(p, file = paste0(output_dir, "/figures/r_objects/del17p13/de_genes/", gene_nam, ".rds"))
  p
}

diff <- diff_all[which(abs(diff_all$stat) > 5),]
geneList <- as.character(diff$Symbol)
geneList <- geneList[-which(geneList %in% "")]

lapply(geneList, gene_count)

[[1]]

Past versions of “single gene counts-1.png”

[[2]]

Past versions of “single gene counts-2.png”

[[3]]

Past versions of “single gene counts-3.png”

[[4]]

Past versions of “single gene counts-4.png”

[[5]]

Past versions of “single gene counts-5.png”

[[6]]

Past versions of “single gene counts-6.png”

[[7]]

Past versions of “single gene counts-7.png”

[[8]]

Past versions of “single gene counts-8.png”

[[9]]

Past versions of “single gene counts-9.png”

[[10]]

Past versions of “single gene counts-10.png”

[[11]]

Past versions of “single gene counts-11.png”

[[12]]

Past versions of “single gene counts-12.png”

[[13]]

Past versions of “single gene counts-13.png”

[[14]]

Past versions of “single gene counts-14.png”

[[15]]

Past versions of “single gene counts-15.png”

[[16]]

Past versions of “single gene counts-16.png”

[[17]]

Past versions of “single gene counts-17.png”

[[18]]

Past versions of “single gene counts-18.png”

[[19]]

Past versions of “single gene counts-19.png”

[[20]]

Past versions of “single gene counts-20.png”

[[21]]

Past versions of “single gene counts-21.png”

[[22]]

Past versions of “single gene counts-22.png”

[[23]]

Past versions of “single gene counts-23.png”

[[24]]

Past versions of “single gene counts-24.png”

[[25]]

Past versions of “single gene counts-25.png”

[[26]]

Past versions of “single gene counts-26.png”

[[27]]

Past versions of “single gene counts-27.png”

[[28]]

Past versions of “single gene counts-28.png”

[[29]]

Past versions of “single gene counts-29.png”

[[30]]

Past versions of “single gene counts-30.png”

[[31]]

Past versions of “single gene counts-31.png”

[[32]]

Past versions of “single gene counts-32.png”

[[33]]

Past versions of “single gene counts-33.png”

[[34]]

Past versions of “single gene counts-34.png”

[[35]]

Past versions of “single gene counts-35.png”

[[36]]

Past versions of “single gene counts-36.png”

[[37]]

Past versions of “single gene counts-37.png”

[[38]]

Past versions of “single gene counts-38.png”

[[39]]

Past versions of “single gene counts-39.png”

[[40]]

Past versions of “single gene counts-40.png”

[[41]]

Past versions of “single gene counts-41.png”

[[42]]

Past versions of “single gene counts-42.png”

[[43]]

Past versions of “single gene counts-43.png”

[[44]]

Past versions of “single gene counts-44.png”

[[45]]

Past versions of “single gene counts-45.png”

[[46]]

Past versions of “single gene counts-46.png”

[[47]]

Past versions of “single gene counts-47.png”

[[48]]

Past versions of “single gene counts-48.png”

[[49]]

Past versions of “single gene counts-49.png”

[[50]]

Past versions of “single gene counts-50.png”

[[51]]

Past versions of “single gene counts-51.png”

[[52]]

Past versions of “single gene counts-52.png”

[[53]]

Past versions of “single gene counts-53.png”

[[54]]

Past versions of “single gene counts-54.png”

[[55]]

Past versions of “single gene counts-55.png”

[[56]]

Past versions of “single gene counts-56.png”

[[57]]

Past versions of “single gene counts-57.png”

[[58]]

Past versions of “single gene counts-58.png”

[[59]]

Past versions of “single gene counts-59.png”

[[60]]

Past versions of “single gene counts-60.png”

[[61]]

Past versions of “single gene counts-61.png”

[[62]]

Past versions of “single gene counts-62.png”

Interesting top genes

A list of genes, with manual checked expression differences: Downregulated in Del17p13: + VPS53 + MPDU1 + DLV2 + SENP3 + TNK1 + ZZEF1 + C1QBP + ZBTB4 + GLOD4 + SMG6 + PFAS + TNFSF12 + MPRIP + RABEBPI + UBE2G1 + CTNS + FXR2 + more..

Upregulated in Del17p13: ..

Gene set enrichment analysis

Gene sets

#load gene set collection
#Hallmark
gsc <- loadGSC("/home/almut/Dokumente/masterarbeit/data/h.all.v6.0.symbols.gmt", type="gmt")
#Kegg
gsc_Kegg <- loadGSC("/home/almut/Dokumente/masterarbeit/data/c2.cp.kegg.v6.0.symbols.gmt", type="gmt")


#get all de outputs
load(paste0(output_dir,"/desRes_15112019.RData"))
difftab <- function(condition){
  dataTab <- data.frame(res_list[[condition]])
  dataTab$ID <- rownames(dataTab)
  #filter using pvalues
  dataTab <- dataTab %>%
    arrange(padj) %>%
    mutate(Symbol = rowData(ddsCLL[ID,])$symbol)# %>%
    #filter(abs(log2FoldChange) > 2)
  dataTab <- dataTab[!duplicated(dataTab$Symbol),]
  dataTab <- dataTab[!is.na(dataTab$Symbol),]
  rownames(dataTab) <- dataTab$ID
   dataTab
}

diff_res <- difftab(variant)

Run piano

gmtFile <- gsc_Kegg
diff_res <- diff_res[-which(diff_res$Symbol %in% c("", NA)),]

geneNam <- diff_res$Symbol 
pVal <- diff_res$padj
logFold <- diff_res$log2FoldChange
stat <- diff_res$stat
gsTab <- data.frame(gene = geneNam, stat = stat, logFold = logFold)
gsaTab <- data.frame(row.names = gsTab$gene, stat = gsTab$stat)
res <- runGSA(geneLevelStats = gsaTab,
                    geneSetStat = "gsea",
                    adjMethod = "fdr", gsc=gmtFile,
                    signifMethod = "geneSampling",
                    nPerm = 50000,
                    gsSizeLim=c(1, Inf))

Running gene set analysis:

Checking arguments...done!

*** Please note that running the GSEA-method may take a substantial amount of time! ***

Final gene/gene-set association: 4271 genes and 186 gene sets

  Details:

  Calculating gene set statistics from 4271 out of 21668 gene-level statistics

  Using all 21668 gene-level statistics for significance estimation

  Removed 995 genes from GSC due to lack of matching gene statistics

  Removed 0 gene sets containing no genes after gene removal

  Removed additionally 0 gene sets not matching the size limits

  Loaded additional information for 186 gene sets

Calculating gene set statistics...done!
Calculating gene set significance...done!
Adjusting for multiple testing...done!

saveRDS(res, file=paste0(output_dir,"/enrichment/", variant, "_gsea.rds"))
res <- readRDS(paste0(output_dir,"/enrichment/", variant, "_gsea.rds"))
  
Res_up <- arrange(GSAsummaryTable(res), `p adj (dist.dir.up)`)            
Res_dn <- arrange(GSAsummaryTable(res), `p adj (dist.dir.dn)`)
  

#Plot
resPlot <- Res_up[, c(1:3,5,8,9)]
resPlot_dn <- Res_dn[, c(1:3,7,8,9)]
colnames(resPlot) <- c("pathway", "gene_number", "stat", "p.adj","genes_up" , "genes_dn")
colnames(resPlot_dn) <- c("pathway", "gene_number", "stat", "p.adj","genes_up" , "genes_dn")


enrichPlot <- resPlot_dn[c(1:6),] %>% mutate(log10Padj = -log10(p.adj)) %>% arrange(-row_number())
enrichPlot$pathway <- gsub("KEGG_", "",  enrichPlot$pathway)
enrichPlot$pathway <- gsub("_", " ",  enrichPlot$pathway)
enrichPlot$pathway <- tolower(enrichPlot$pathway)

#plot
p <- ggbarplot(enrichPlot, x = "pathway", y = "gene_number",
          fill = "p.adj",   
          orientation = "vertical",
          label = enrichPlot$pathway, lab.col = "white", lab.size = 4, lab.vjust = 0, lab.hjust = 1.01,
          sort.val = "none",         
          sort.by.groups = FALSE,     
          xlab = "Pathway",
          ylab = "gene number",
          legend.title = "p.adj",
          rotate = TRUE, font.x = 20, font.y = 20, font.legend = 20, legend = "right", 
          title = "del17p13 - Kegg down",
          ggtheme = theme_pubr()) + 
  font("x.text", size = 16) + 
  font("y.text", size = 0) +
  font("title", size = 20, face = "bold") + 
  gradient_color(c("blue", "white", "red"))
  
#ggsave(file=paste0(figure_dir,"/GSEA_", variant, "_Kegg.svg"), plot=p, width=14, height=7)
  p

Past versions of piano-1.png

saveRDS(p, file = paste0(output_dir, "/figures/r_objects/del17p13/del17p13_enrichment_dn.rds"))


enrichPlot <- resPlot[c(1:6),] %>% mutate(log10Padj = -log10(p.adj)) %>% arrange(-row_number())
enrichPlot$pathway <- gsub("KEGG_", "",  enrichPlot$pathway)
enrichPlot$pathway <- gsub("_", " ",  enrichPlot$pathway)
enrichPlot$pathway <- tolower(enrichPlot$pathway)

#plot
p <- ggbarplot(enrichPlot, x = "pathway", y = "gene_number",
          fill = "p.adj",   
          orientation = "vertical",
          label = enrichPlot$pathway, lab.col = "white", lab.size = 4, lab.vjust = 0, lab.hjust = 1.01,
          sort.val = "none",         
          sort.by.groups = FALSE,     
          xlab = "Pathway",
          ylab = "gene number",
          legend.title = "p.adj",
          rotate = TRUE, font.x = 20, font.y = 20, font.legend = 20, legend = "right", 
          title = "del17p13 - Kegg up",
          ggtheme = theme_pubr()) + 
  font("x.text", size = 16) + 
  font("y.text", size = 0) +
  font("title", size = 20, face = "bold") + 
  gradient_color(c("blue", "white", "red"))
  
#ggsave(file=paste0(figure_dir,"/GSEA_", variant, "_Kegg.svg"), plot=p, width=14, height=7)
  p

saveRDS(p, file = paste0(output_dir, "/figures/r_objects/del17p13/del17p13_enrichment_up.rds"))

Session information

sessionInfo()

R version 3.6.3 (2020-02-29)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.7 LTS

Matrix products: default
BLAS:   /usr/lib/libblas/libblas.so.3.6.0
LAPACK: /usr/lib/lapack/liblapack.so.3.6.0

locale:
 [1] LC_CTYPE=de_DE.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=de_DE.UTF-8        LC_COLLATE=de_DE.UTF-8    
 [5] LC_MONETARY=de_DE.UTF-8    LC_MESSAGES=de_DE.UTF-8   
 [7] LC_PAPER=de_DE.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
 [1] grid      parallel  stats4    stats     graphics  grDevices utils    
 [8] datasets  methods   base     

other attached packages:
 [1] here_0.1                    RColorBrewer_1.1-2         
 [3] ggpubr_0.2                  magrittr_1.5               
 [5] gridExtra_2.3               circlize_0.4.6             
 [7] gtable_0.3.0                ComplexHeatmap_2.0.0       
 [9] genefilter_1.66.0           reshape2_1.4.3             
[11] piano_2.0.2                 ggsci_2.9                  
[13] forcats_0.4.0               stringr_1.4.0              
[15] dplyr_0.8.1                 purrr_0.3.2                
[17] readr_1.3.1                 tidyr_0.8.3                
[19] tibble_2.1.3                ggplot2_3.1.1              
[21] tidyverse_1.2.1             DESeq2_1.24.0              
[23] SummarizedExperiment_1.14.0 DelayedArray_0.10.0        
[25] BiocParallel_1.18.0         matrixStats_0.54.0         
[27] Biobase_2.44.0              GenomicRanges_1.36.0       
[29] GenomeInfoDb_1.20.0         IRanges_2.18.1             
[31] S4Vectors_0.22.0            BiocGenerics_0.30.0        

loaded via a namespace (and not attached):
  [1] readxl_1.3.1           backports_1.1.4        Hmisc_4.2-0           
  [4] fastmatch_1.1-0        workflowr_1.4.0        plyr_1.8.4            
  [7] igraph_1.2.4.1         lazyeval_0.2.2         shinydashboard_0.7.1  
 [10] splines_3.6.3          digest_0.6.19          htmltools_0.3.6       
 [13] gdata_2.18.0           checkmate_1.9.3        memoise_1.1.0         
 [16] cluster_2.1.1          limma_3.40.2           annotate_1.62.0       
 [19] modelr_0.1.4           colorspace_1.4-1       blob_1.1.1            
 [22] rvest_0.3.4            haven_2.1.0            xfun_0.7              
 [25] crayon_1.3.4           RCurl_1.95-4.12        jsonlite_1.6          
 [28] survival_2.44-1.1      glue_1.3.1             zlibbioc_1.30.0       
 [31] XVector_0.24.0         GetoptLong_0.1.7       shape_1.4.4           
 [34] scales_1.0.0           DBI_1.0.0              relations_0.6-8       
 [37] Rcpp_1.0.1             xtable_1.8-4           htmlTable_1.13.1      
 [40] clue_0.3-57            foreign_0.8-76         bit_1.1-14            
 [43] Formula_1.2-3          DT_0.17                htmlwidgets_1.3       
 [46] httr_1.4.0             fgsea_1.10.0           gplots_3.0.1.1        
 [49] acepack_1.4.1          pkgconfig_2.0.2        XML_3.98-1.20         
 [52] nnet_7.3-15            locfit_1.5-9.1         labeling_0.3          
 [55] tidyselect_0.2.5       rlang_0.3.4            later_0.8.0           
 [58] AnnotationDbi_1.46.0   munsell_0.5.0          cellranger_1.1.0      
 [61] tools_3.6.3            visNetwork_2.0.7       cli_1.1.0             
 [64] generics_0.0.2         RSQLite_2.1.1          broom_0.5.2           
 [67] evaluate_0.14          yaml_2.2.0             knitr_1.23            
 [70] bit64_0.9-7            fs_1.3.1               caTools_1.17.1.2      
 [73] nlme_3.1-152           whisker_0.3-2          mime_0.7              
 [76] slam_0.1-45            xml2_1.2.0             compiler_3.6.3        
 [79] rstudioapi_0.10        png_0.1-7              marray_1.62.0         
 [82] geneplotter_1.62.0     stringi_1.4.3          lattice_0.20-38       
 [85] Matrix_1.3-2           shinyjs_1.0            pillar_1.4.1          
 [88] GlobalOptions_0.1.0    data.table_1.12.2      bitops_1.0-6          
 [91] httpuv_1.5.1           R6_2.4.0               latticeExtra_0.6-28   
 [94] promises_1.0.1         KernSmooth_2.23-15     gtools_3.8.1          
 [97] assertthat_0.2.1       rprojroot_1.3-2        rjson_0.2.20          
[100] withr_2.1.2            GenomeInfoDbData_1.2.1 hms_0.4.2             
[103] rpart_4.1-15           rmarkdown_1.13         git2r_0.25.2          
[106] sets_1.0-18            shiny_1.3.2            lubridate_1.7.4       
[109] base64enc_0.1-3