TP53.Rmd

---
title: "TP53"
author: "almut"
date: "16 Juni 2019"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## TP53 signature

# Differentially expressed genes

### 1. Differential expression analysis
load packages
```{r packages, warning=FALSE, message=FALSE}
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)

```


load data
```{r datasets}
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
```{r var gene expression}

#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
```{r complex heatmap, fig.width= 15, fig.height=18}

#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 ) 
#dev.off()

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


# Sample and gene specific expression - top genes

```{r single gene counts, fig.width=6, fig.height=5}

gene_count <- function(gene_nam){
   geneEnsID <- rownames(ddsCLL)[which(rowData(ddsCLL)$symbol %in% gene_nam)]
  gc <- plotCounts(ddsCLL, gene = geneEnsID, intgroup = variant, returnData=TRUE)
  p <- ggboxplot(gc, x = variant, y = "count",
          color = variant,
          size = 1.2,
          palette = "jco",
          outlier.shape = NA, 
          add = "jitter",
          add.params = list(size = 2.5),
          yscale = "log10",
          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/TP53/de_genes/", gene_nam, ".rds"))
  p
}

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

lapply(geneList, gene_count)


```

## Gene set enrichment analysis

Gene sets
```{r 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)

gene_list_id <- bitr(diff_res$ID, fromType = "ENSEMBL",
        toType = c("ENTREZID", "SYMBOL"),
        OrgDb = org.Hs.eg.db)
names(gene_list_id) <- c("ID", "ENTREZID", "Symbol")
diff_id <- left_join(gene_list_id, diff_res)
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)
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
```{r clusterProfiler}
barplot(kk, showCategory=5)
barplot(em, showCategory=5)

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

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

ridgeplot(em2)
ridgeplot(kk2)

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

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
```{r network, fig.width=10, fig.height=10}
# 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")

p_net  

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

pnet_kegg

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
```{r heatplot, fig.width=20, fig.height=5}
heatplot(em2, foldChange=gene_lfc, showCategory = 3)
heatplot(kk2x, foldChange=gene_lfc, showCategory = 3 ) 
```