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

#differentially expressed genes between SF3B1 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$log2FoldChange) > 2 & abs(diff_all$stat) > 4),]


mutStatus <- data.frame(colData(ddsCLL)) %>% arrange(SF3B1)
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(SF3B1 = 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 = 25),
              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/SF3B1/SF3B1_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",
          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/SF3B1/de_genes/", gene_nam, ".rds"))
  p
}

geneList <- as.character(diff$Symbol[1:20])


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”

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 <- 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: 114 genes and 122 gene sets

  Details:

  Calculating gene set statistics from 114 out of 509 gene-level statistics

  Using all 509 gene-level statistics for significance estimation

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

  Removed 64 gene sets containing no genes after gene removal

  Removed additionally 0 gene sets not matching the size limits

  Loaded additional information for 122 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 = "SF3B1 - 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/SF3B1/SF3B1_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