1. Differential expression analysis

load packages

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(clusterProfiler)
library(msigdbr)
library(org.Hs.eg.db)
library(enrichplot)
library(here)
library(VennDiagram)

load data

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

#differentially expressed genes between TP53 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


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

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

Expression matrix

#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

Expression signature

#colors
colors = colorRamp2(c(-4,-2,0,2,4), c("#2166ac","#4393c3", "#f7f7f7", "#d6604d","#b2182b"))
annocol <- get_palette("jco", 10)
annocolor <- list(TP53 = 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) > 6),]
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 = 17),
              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 ) 

Past versions of “complex heatmap-1.png”

#dev.off()

saveRDS(h1, file = paste0(output_dir, "/figures/r_objects/TP53/TP53_heatmap.rds"))

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

#clusterProfiler
diff_res <- diff_res[-which(diff_res$Symbol %in% c("", NA)),]
gene_list <- diff_res$stat %>% set_names(diff_res$Symbol)
gene_list <- sort(gene_list, decreasing = TRUE)
gene_lfc <- diff_res$log2FoldChange %>% set_names(diff_res$Symbol)
gene_lfc <- sort(gene_lfc, decreasing = TRUE)
de_gene <- diff_res %>% filter(padj < 0.01) 
de_gene <- de_gene$Symbol

de_ens <- diff_res %>% filter(padj < 0.01)
de_ens <- de_ens$ID
#Get Gene IDs
gene_id <- bitr(de_ens, fromType = "ENSEMBL",
        toType = c("ENTREZID", "SYMBOL"),
        OrgDb = org.Hs.eg.db)

'select()' returned 1:1 mapping between keys and columns

Warning in bitr(de_ens, fromType = "ENSEMBL", toType = c("ENTREZID",
"SYMBOL"), : 11.54% of input gene IDs are fail to map...

gene_list_id <- bitr(diff_res$ID, fromType = "ENSEMBL",
        toType = c("ENTREZID", "SYMBOL"),
        OrgDb = org.Hs.eg.db)

'select()' returned 1:many mapping between keys and columns

Warning in bitr(diff_res$ID, fromType = "ENSEMBL", toType = c("ENTREZID", :
18.06% of input gene IDs are fail to map...

names(gene_list_id) <- c("ID", "ENTREZID", "Symbol")
diff_id <- left_join(gene_list_id, diff_res)

Joining, by = c("ID", "Symbol")

gene_list_id <- diff_id$stat %>% set_names(diff_id$ENTREZID)
gene_list_id <- sort(gene_list_id, decreasing = TRUE)
gene_lfc_id <- diff_id$log2FoldChange %>% set_names(diff_id$ENTREZID)
gene_lfc_id <- sort(gene_lfc_id, decreasing = TRUE)

#convert gsc
m_t2g <- msigdbr(species = "Homo sapiens", category = "H") %>% 
  dplyr::select(gs_name, human_gene_symbol)

#Hallmark
em2 <- GSEA(gene_list, TERM2GENE = m_t2g, pvalueCutoff = 0.1)

preparing geneSet collections...

GSEA analysis...

Warning in fgsea(pathways = geneSets, stats = geneList, nperm = nPerm, minSize = minGSSize, : There are ties in the preranked stats (0% of the list).
The order of those tied genes will be arbitrary, which may produce unexpected results.

leading edge analysis...

done...

em <- enricher(de_gene, TERM2GENE = m_t2g)

#Kegg
kk <- enrichKEGG(gene_id$ENTREZID,
                 organism     = 'hsa',
                 pvalueCutoff = 0.2)

kk2 <- gseKEGG(geneList     = gene_list_id,
               organism     = 'hsa',
               nPerm        = 1000,
               minGSSize    = 50,
               pvalueCutoff = 0.2,
               verbose      = FALSE)

kk2x <- setReadable(kk2, 'org.Hs.eg.db', 'ENTREZID')

Visualize ClusterProfiler results

barplot(kk, showCategory=5)

Past versions of clusterProfiler-1.png

barplot(em, showCategory=5)

Past versions of clusterProfiler-2.png

dot1 <- clusterProfiler::dotplot(em2, showCategory=10) + ggtitle("GSEA for TP53") +
  theme_pubr() +
  theme(legend.position="right") + 
  theme(plot.title = element_text(face = "bold")) 

wrong orderBy parameter; set to default `orderBy = "x"`

dot1

Past versions of clusterProfiler-3.png

dotplot(em, showCategory=10) + ggtitle("Enrichment for TP53")

wrong orderBy parameter; set to default `orderBy = "x"`

Past versions of clusterProfiler-4.png

dotplot(kk2, showCategory=10) + ggtitle("GSEA for TP53")

wrong orderBy parameter; set to default `orderBy = "x"`

Past versions of clusterProfiler-5.png

dot2 <- clusterProfiler::dotplot(kk, showCategory=10) + ggtitle("Enrichment for TP53") +
  theme_pubr() +
  theme(legend.position="right") +
  theme(plot.title = element_text(face = "bold"))

wrong orderBy parameter; set to default `orderBy = "x"`

dot2

Past versions of clusterProfiler-6.png

ridgeplot(em2)

Picking joint bandwidth of 0.302

Past versions of clusterProfiler-7.png

ridgeplot(kk2)

Picking joint bandwidth of 0.287

Past versions of clusterProfiler-8.png

gseaplot2(em2, geneSetID = 3, title = em2$Description[3])

Past versions of clusterProfiler-9.png

gseaplot2(kk2, geneSetID = 2, title = kk2$Description[2])

Past versions of clusterProfiler-10.png

saveRDS(dot1, file = paste0(output_dir, "/figures/r_objects/TP53/enrich_dot_hm.rds"))
saveRDS(dot2, file = paste0(output_dir, "/figures/r_objects/TP53/enrich_dot2.rds"))

network plot

# Networks Hallmark
 em2_sub <- em2
 em2_sub@result <- em2@result[which(em2@result$Description %in% c("HALLMARK_OXIDATIVE_PHOSPHORYLATION",
                                                                    "HALLMARK_DNA_REPAIR",
                                                                    "HALLMARK_G2M_CHECKPOINT")),]
p_net <- cnetplot(em2_sub, categorySize="pvalue", foldChange=gene_lfc) + 
  scale_colour_gradientn(colors = c("#581845", "#900C3F", "#C70039", "#FF5733", "#FFC300", "#DAF7A6")) + 
  guides(size = FALSE) + 
  labs(color = "logFC")

Scale for 'colour' is already present. Adding another scale for
'colour', which will replace the existing scale.

p_net  

Past versions of network-1.png

# Networks KEGG
 kk2_sub <- kk2x
 kk2_sub@result <- kk2x@result[which(kk2x@result$Description %in% c("p53 signaling pathway",
                                                                    "Oxidative phosphorylation",
                                                                    "Transcriptional misregulation in cancer"
                                                                    )),]

pnet_kegg <- cnetplot(kk2_sub, categorySize="pvalue", foldChange=gene_lfc) + 
  scale_color_gradient(high="blue", low="red") +
  guides(size = FALSE) + 
  labs(color = "logFC")

Scale for 'colour' is already present. Adding another scale for
'colour', which will replace the existing scale.

pnet_kegg

Past versions of network-2.png

saveRDS(pnet_kegg, file = paste0(output_dir, "/figures/r_objects/TP53/enrich_net_kegg.rds"))
saveRDS(p_net, file = paste0(output_dir, "/figures/r_objects/TP53/enrich_net_hm.rds"))

heatplot

heatplot(em2, foldChange=gene_lfc, showCategory = 3)

Past versions of heatplot-1.png

heatplot(kk2x, foldChange=gene_lfc, showCategory = 3 ) 

Past versions of heatplot-2.png

Overlap with del17

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

rownames(del17) <- del17$X
del17 <- del17[which(del17$padj < 0.01 ),-1]
del17<- del17[,-1]
TP53 <- diff_all
TP53  <- TP53[-which(TP53$Symbol %in% ""),]

venntab <- unique(c(as.vector(del17$Symbol), as.vector(TP53$Symbol)))
vennTab_new <- as.data.frame(venntab) %>% mutate(del17p13 = ifelse(as.character(venntab) %in% del17$Symbol, 1, 0)) %>% mutate(TP53 = ifelse(as.character(venntab) %in% TP53$Symbol, 1, 0))

TP53_only <- vennTab_new %>% filter(del17p13 == 0 & TP53 == 1) %>% dplyr::select(venntab)
del17_only <- vennTab_new %>% filter(del17p13 == 1 & TP53 == 0) %>% dplyr::select(venntab)


vennTab_new <- vennTab_new[,-1]
rownames(vennTab_new) <- venntab
del17_area <- vennTab_new %>% filter(del17p13 == 1)
TP53_area <- vennTab_new %>% filter(TP53 == 1)
Both_area <- vennTab_new %>% filter(del17p13 == 1 & TP53 == 1)

#vennList <- list("del17p13" = as.vector(del17$Symbol), "TP53" = as.vector(TP53$Symbol))
#venn <- venn.diagram(vennList, filename = NULL)

venn_pair <- draw.pairwise.venn(nrow(del17_area), nrow(TP53_area), nrow(Both_area), c("del17p13", "TP53"), fill = c("#00AFBB", "#E7B800"), cex = rep(2, 3), cat.cex = rep(2, 2),  ind = TRUE)


grid.draw(venn_pair)

Past versions of overlap-1.png

grid.newpage()
saveRDS(venn_pair, file = paste0(output_dir, "/figures/r_objects/TP53/venn_pair.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] VennDiagram_1.6.20          futile.logger_1.4.3        
 [3] here_0.1                    enrichplot_1.4.0           
 [5] org.Hs.eg.db_3.8.2          AnnotationDbi_1.46.0       
 [7] msigdbr_7.0.1               clusterProfiler_3.12.0     
 [9] RColorBrewer_1.1-2          ggpubr_0.2                 
[11] magrittr_1.5                gridExtra_2.3              
[13] circlize_0.4.6              gtable_0.3.0               
[15] ComplexHeatmap_2.0.0        genefilter_1.66.0          
[17] reshape2_1.4.3              piano_2.0.2                
[19] ggsci_2.9                   forcats_0.4.0              
[21] stringr_1.4.0               dplyr_0.8.1                
[23] purrr_0.3.2                 readr_1.3.1                
[25] tidyr_0.8.3                 tibble_2.1.3               
[27] ggplot2_3.1.1               tidyverse_1.2.1            
[29] DESeq2_1.24.0               SummarizedExperiment_1.14.0
[31] DelayedArray_0.10.0         BiocParallel_1.18.0        
[33] matrixStats_0.54.0          Biobase_2.44.0             
[35] GenomicRanges_1.36.0        GenomeInfoDb_1.20.0        
[37] IRanges_2.18.1              S4Vectors_0.22.0           
[39] BiocGenerics_0.30.0        

loaded via a namespace (and not attached):
  [1] shinydashboard_0.7.1   tidyselect_0.2.5       RSQLite_2.1.1         
  [4] htmlwidgets_1.3        munsell_0.5.0          DT_0.17               
  [7] withr_2.1.2            colorspace_1.4-1       GOSemSim_2.10.0       
 [10] knitr_1.23             rstudioapi_0.10        DOSE_3.10.2           
 [13] labeling_0.3           git2r_0.25.2           slam_0.1-45           
 [16] urltools_1.7.3         GenomeInfoDbData_1.2.1 polyclip_1.10-0       
 [19] bit64_0.9-7            farver_2.0.3           rprojroot_1.3-2       
 [22] generics_0.0.2         lambda.r_1.2.3         xfun_0.7              
 [25] sets_1.0-18            R6_2.4.0               clue_0.3-57           
 [28] graphlayouts_0.6.0     locfit_1.5-9.1         bitops_1.0-6          
 [31] fgsea_1.10.0           gridGraphics_0.5-0     assertthat_0.2.1      
 [34] promises_1.0.1         scales_1.0.0           ggraph_2.0.2          
 [37] nnet_7.3-16            tidygraph_1.1.2        workflowr_1.4.0       
 [40] rlang_0.3.4            GlobalOptions_0.1.0    splines_3.6.3         
 [43] lazyeval_0.2.2         acepack_1.4.1          broom_0.5.2           
 [46] europepmc_0.3          checkmate_1.9.3        BiocManager_1.30.4    
 [49] yaml_2.2.0             modelr_0.1.4           backports_1.1.4       
 [52] httpuv_1.5.1           qvalue_2.16.0          Hmisc_4.2-0           
 [55] tools_3.6.3            relations_0.6-8        ggplotify_0.0.5       
 [58] gplots_3.0.1.1         ggridges_0.5.2         Rcpp_1.0.1            
 [61] plyr_1.8.4             base64enc_0.1-3        visNetwork_2.0.7      
 [64] progress_1.2.2         zlibbioc_1.30.0        RCurl_1.95-4.12       
 [67] prettyunits_1.0.2      rpart_4.1-15           GetoptLong_0.1.7      
 [70] viridis_0.5.1          cowplot_0.9.4          haven_2.1.0           
 [73] ggrepel_0.8.1          cluster_2.1.2          fs_1.3.1              
 [76] futile.options_1.0.1   data.table_1.12.2      DO.db_2.9             
 [79] triebeard_0.3.0        whisker_0.3-2          hms_0.4.2             
 [82] shinyjs_1.0            mime_0.7               evaluate_0.14         
 [85] xtable_1.8-4           XML_3.98-1.20          readxl_1.3.1          
 [88] shape_1.4.4            compiler_3.6.3         KernSmooth_2.23-15    
 [91] crayon_1.3.4           htmltools_0.3.6        later_0.8.0           
 [94] Formula_1.2-3          geneplotter_1.62.0     lubridate_1.7.4       
 [97] DBI_1.0.0              formatR_1.7            tweenr_1.0.1          
[100] MASS_7.3-54            Matrix_1.3-4           cli_1.1.0             
[103] marray_1.62.0          gdata_2.18.0           igraph_1.2.4.1        
[106] pkgconfig_2.0.2        rvcheck_0.1.8          foreign_0.8-76        
[109] xml2_1.2.0             annotate_1.62.0        XVector_0.24.0        
[112] rvest_0.3.4            digest_0.6.19          rmarkdown_1.13        
[115] cellranger_1.1.0       fastmatch_1.1-0        htmlTable_1.13.1      
[118] shiny_1.3.2            gtools_3.8.1           rjson_0.2.20          
[121] nlme_3.1-152           jsonlite_1.6           viridisLite_0.3.0     
[124] limma_3.40.2           pillar_1.4.1           lattice_0.20-38       
[127] httr_1.4.0             survival_2.44-1.1      GO.db_3.8.2           
[130] glue_1.3.1             UpSetR_1.4.0           png_0.1-7             
[133] bit_1.1-14             ggforce_0.3.1          stringi_1.4.3         
[136] blob_1.1.1             latticeExtra_0.6-28    caTools_1.17.1.2      
[139] memoise_1.1.0