Del13q14 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 <- "del13q14"
#filter for patients without NA in variant
ddsCLL <- ddsCLL[, !is.na(colData(ddsCLL)[,variant])]

#differentially expressed genes between del13q14 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) > 5) ,]


mutStatus <- data.frame(colData(ddsCLL)) %>% arrange(del13q14)
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(del13q14 = 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"), annotation_legend_param = list(title_gp = gpar(fontsize = 50), labels_gp = gpar(fontsize = 45),  grid_height = unit(1.9, "cm"), grid_width = unit(1.9, "cm")))


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 = 60, fontface = "bold"), 
              heatmap_legend_param = list(title = "expr", 
                                          title_gp = gpar(fontsize = 50), 
                                          grid_height = unit(1.9, "cm"), 
                                          grid_width = unit(1.9, "cm"), 
                                          gap = unit(2, "cm"), 
                                          labels_gp = gpar(fontsize = 45)), 
              column_dend_height = unit(2.5, "cm"),
              show_row_dend = FALSE, 
              show_column_names = FALSE , 
              show_row_names = TRUE, 
              row_names_gp = gpar(fontsize = 40),
              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/del13q14/del13q14_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 <- counts(ddsCLL)[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",
          palette = "jco",
          add = "mean_sd",
          title = paste(gene_nam),
          font.x = 18, font.y = 18, font.legend = 16, 
          ylab = "normalized counts") + font("xy.text", size = 15) + font("title", size = 20, face = "bold")
  #ggsave(file=paste0(figure_dir, "/tri12/genetic_interaction_", gene_nam, ".svg"), plot=p, width=6, height=5)
  saveRDS(p, file = paste0(output_dir, "/figures/r_objects/del13q14/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]]

[[60]]

[[61]]

[[62]]

[[63]]

[[64]]

[[65]]

[[66]]

[[67]]

[[68]]

[[69]]

[[70]]

[[71]]

[[72]]

[[73]]

[[74]]

[[75]]

[[76]]

[[77]]

[[78]]

[[79]]

[[80]]

[[81]]

[[82]]

[[83]]

[[84]]

[[85]]

[[86]]

[[87]]

[[88]]

[[89]]

[[90]]

[[91]]

[[92]]

[[93]]

[[94]]

[[95]]

[[96]]

[[97]]

[[98]]

[[99]]

[[100]]

[[101]]

[[102]]

[[103]]

[[104]]

[[105]]

[[106]]

[[107]]

[[108]]

[[109]]

[[110]]

[[111]]

[[112]]

[[113]]

[[114]]

[[115]]

[[116]]

[[117]]

[[118]]

[[119]]

[[120]]

[[121]]

[[122]]

[[123]]

[[124]]

[[125]]

[[126]]

[[127]]

[[128]]

[[129]]

[[130]]

[[131]]

[[132]]

[[133]]

[[134]]

[[135]]

[[136]]

[[137]]

[[138]]

[[139]]

[[140]]

[[141]]

[[142]]

[[143]]

[[144]]

[[145]]

[[146]]

[[147]]

[[148]]

[[149]]

[[150]]

[[151]]

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")

Run piano

gmtFile <- gsc_Kegg

diff_res <- diff_all
diff_res$chromosome <- rowData(RNAnorm)[rownames(diff_res),]$chromosome
diff_res <- diff_res[-which(diff_res$Symbol %in% c("", NA)),]
#Filter chromosome 13 - we are not interested in the dosage effect 
diff_res <- diff_res[-which(diff_res$chromosome %in% c("13")),]

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: 188 genes and 150 gene sets

  Details:

  Calculating gene set statistics from 188 out of 786 gene-level statistics

  Using all 786 gene-level statistics for significance estimation

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

  Removed 36 gene sets containing no genes after gene removal

  Removed additionally 0 gene sets not matching the size limits

  Loaded additional information for 150 gene sets

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

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)) 

#plot
p <- ggbarplot(enrichPlot, x = "pathway", y = "log10Padj",
          fill = "gene_number",         
          color = "white",            
          palette =  "gsea",            
          sort.val = "asc",         
          sort.by.groups = FALSE,     
          ylab = "-log10(padj)",
          legend.title = "#diff.genes",
          rotate = TRUE,
          font.x = 20, font.y = 20, font.legend = 20, legend = "right", 
          title = "del13q14 - Kegg",
          ggtheme = theme_pubr()) + 
  font("xy.text", size = 16) + 
  font("title", size = 20, face = "bold")
  
#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/del13q14/del13q14_enrichment.rds"))

Session information

sessionInfo()

R version 3.6.0 (2019-04-26)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.6 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.0          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.0          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-71         bit_1.1-14            
 [43] Formula_1.2-3          DT_0.7                 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-12            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.0            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-140           whisker_0.3-2          mime_0.7              
 [76] slam_0.1-45            xml2_1.2.0             compiler_3.6.0        
 [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.2-17          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