IGHV.Rmd

---
title: "IGHV"
author: "aluetge"
date: "2019-06-29"
output: html_document
---


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

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

```


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

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


mutStatus <- data.frame(colData(ddsCLL)) %>% arrange(IGHV)
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.height= 18, fig.width= 14.22}

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


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

ha_col <- HeatmapAnnotation(df = feature, col = annocolor, annotation_height = unit(c(rep(1.9, 1)), "cm"), annotation_legend_param = list(title_gp = gpar(fontsize = 17), labels_gp = gpar(fontsize = 15),  grid_height = unit(0.7, "cm"), grid_width = unit(0.3, "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 = 20, fontface = "bold"), 
              heatmap_legend_param = list(title = "expr", 
                                          title_gp = gpar(fontsize = 20), 
                                          grid_height = unit(0.7, "cm"), 
                                          grid_width = unit(0.3, "cm"), 
                                          gap = unit(2, "cm"), 
                                          labels_gp = gpar(fontsize = 15)), 
              column_dend_height = unit(1, "cm"),
              show_row_dend = FALSE, 
              show_column_names = FALSE , 
              row_title_gp = gpar(fontsize = 0),
              show_row_names = FALSE, 
              row_names_gp = gpar(fontsize = 0),
              top_annotation = ha_col)

#Annotate top 50 genes
sub_names <- diff[1:50,"Symbol"]
sub_names <- sub_names[-which(sub_names %in% "")]
sub_names <- c(as.character(sub_names), "CD38", "LPL", "ZAP70", "SEPT10", "ADAM29", "PEG10")
#
geneIDs <- which(rownames(exprVariant.new) %in% sub_names)
labels <- rownames(exprVariant.new)[geneIDs]
ha_genes <- rowAnnotation(link = row_anno_link(at = geneIDs, labels = labels, labels_gp = gpar(fontsize = 13)), width = unit(2.5, "cm"))


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

p1 <- draw(h1 + ha_genes ) 
#dev.off()
saveRDS(p1, file = paste0(output_dir, "/figures/r_objects/ighv/ighv_heatmap.rds"))
```


## Top genes

```{r single gene counts, fig.width=6, fig.height=5}
#function to create stripchart plots for specific genes
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",
          add = "jitter",
          add.params = list(size = 2.5),
          outlier.shape = NA, 
          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/ighv/de_genes/", gene_nam, ".rds"))
  p
}

diff <- diff_all[which(abs(diff_all$stat) > 10),]
geneList <- as.character(diff$Symbol)
geneList <- geneList[-which(geneList %in% "")]
geneList <- c(geneList, "CD38", "LPL", "ZAP70", "SEPT10", "ADAM29", "PEG10", "EGR3", "NFAT5", "BCAT1")

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 IGHV") +
  theme_pubr() +
  theme(legend.position="right") + 
  theme(plot.title = element_text(face = "bold")) 
dot1
dotplot(em, showCategory=10) + ggtitle("Enrichment for IGHV")

dotplot(kk2, showCategory=10) + ggtitle("GSEA for IGHV")
dot2 <- clusterProfiler::dotplot(kk, showCategory=10) + ggtitle("Enrichment for IGHV") +
  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/ighv/enrich_dot_hm.rds"))
saveRDS(dot2, file = paste0(output_dir, "/figures/r_objects/ighv/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_TNFA_SIGNALING_VIA_NFKB",
                                                                    "HALLMARK_KRAS_SIGNALING_DN",
                                                                    "HALLMARK_P53_PATHWAY",
                                                                    "HALLMARK_ANGIOGENESIS")),]
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("B cell receptor signaling pathway",
                                                                    "Chemokine signaling pathway",
                                                                    "T cell receptor signaling pathway"
                                                                    )),]

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/ighv/enrich_net_kegg.rds"))
saveRDS(p_net, file = paste0(output_dir, "/figures/r_objects/ighv/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)
```






## IG gene status

```{r IG variants}

#load object with variant status
load(paste0(data_dir, "/2018-03-05_IGHV.RData"))
#sort patIDs by ddsCll object
IG_gene <- as.tibble(IG_gene_analysis_Uppsala) %>% mutate(IGHV1_69 = ifelse(Vgene %in% "IGHV1-69", 1, 0), PatID = patientID) %>% dplyr::select(IGHV1_69, PatID)

colData_tib <- as.tibble(colData(ddsCLL)) 

colDat_integrated <- left_join(colData_tib,IG_gene)

colData(ddsCLL)$IGHV1_69 <- as.factor(colDat_integrated$IGHV1_69)


geneEnsID <- rownames(ddsCLL)[which(rowData(ddsCLL)$symbol %in% "IGHV1-69")]
geneNum <- counts(ddsCLL)[geneEnsID,]
mutPat <- as.data.frame(colData(ddsCLL)[, c("IGHV", "IGHV1_69")])
colnames(mutPat) <- c("genotype", "IGHV1_69_variant")
geneDat <- cbind(mutPat, geneNum)
colnames(geneDat) <- c("genotype", "IGHV1_69_variant","counts")
geneDat <- geneDat[!is.na(geneDat$IGHV1_69_variant),]

# geneTib <- as.tibble(geneDat) %>% mutate(PatID = rownames(geneDat)) %>%
#   dplyr::filter(IGHV1_69_variant == 1) %>%
#   arrange(counts)


p <- ggboxplot(geneDat, x = "IGHV1_69_variant", y = "counts",
          color = "IGHV1_69_variant",
          size = 1.2,
          palette = c("#00AFBB", "#FC4E07"),
          add = "jitter",
          add.params = list(size = 2.5),
          outlier.shape = NA,
          yscale = "log10",
          title = "IGHV1-69 gene expression",
          font.x = 20, font.y = 20, font.legend = 20, 
          ylab = "normalized counts", xlab = "IGHV1-69 variant", 
          legend = "right") + font("xy.text", size = 20) + font("title", size = 20, face = "bold")

p

saveRDS(p, file = paste0(output_dir, "/figures/r_objects/ighv/ighv_gene_status.rds"))

```