IGHV and Trisomy12

IGHV and trisomy12 are the most severe genomic modification on transcriptome level. How do they affect each other?

load packages

library(Biobase)
library(ggplot2)
library(genefilter)
library(DESeq2)
library(gridExtra)
library(reshape2)
library(dplyr)
library(geneplotter)
library(RColorBrewer)
library(ComplexHeatmap)
library(circlize)
library(piano)
library(ggpubr)
library(here)
library(clusterProfiler)
library(msigdbr)
library(org.Hs.eg.db)
library(enrichplot)
library(purrr)
library(data.table)
set.seed(1000)

load datasets

data_dir <- here("data")
output_dir <- here("output")
figure_dir <- here("output/figures")

#Countdata
load(paste0(data_dir, "/ddsrnaCLL_150218.RData"))

Heatmap of interacting genes

  #filter by variant
  sig_Genes <- rownames(resSig)
  
  #gene expression data
  geneExpression = assay(RNAnorm)[sig_Genes,]
  rownames(geneExpression) <- rowData(RNAnorm)$symbol[which(rownames(RNAnorm) %in% sig_Genes)]

  #scale and censor
  geneExpression_new <- log2(geneExpression)
  geneExpression_new<- t(scale(t(geneExpression_new)))
  geneExpression_new[geneExpression_new > 4] <- 4
  geneExpression_new[geneExpression_new < -4] <- -4

  mutnames <- c("IGHV-UM", "IGHV-M", "trisomy12", "both")
  
  mutStatus <- data.frame(colData(RNAnorm)) %>% mutate(IGHVnew = ifelse(IGHV %in% "M", 1, 0)) %>% 
    dplyr::select(-IGHV) %>% mutate(IGHV = IGHVnew) %>% 
    dplyr::select_("PatID", "IGHV", "trisomy12") %>% 
    mutate_("namA" = "IGHV", "namB" = "trisomy12") %>% 
    mutate(naA = as.numeric(as.character(namA))) %>% 
    mutate(naB = as.numeric(as.character(namB))) %>% 
    mutate(mut = factor(mutnames[1 + naA + 2 * naB], levels = mutnames)) %>% arrange(mut)

Warning: select_() is deprecated. 
Please use select() instead

The 'programming' vignette or the tidyeval book can help you
to program with select() : https://tidyeval.tidyverse.org
This warning is displayed once per session.

Warning: mutate_() is deprecated. 
Please use mutate() instead

The 'programming' vignette or the tidyeval book can help you
to program with mutate() : https://tidyeval.tidyverse.org
This warning is displayed once per session.

  geneExpression_new <- geneExpression_new[,mutStatus$PatID]


  #colors
  #colors <- colorRampPalette( rev(brewer.pal(11,"RdBu")) )(255)
  colors = colorRamp2(c(-4, -1, 0, 1, 4), c("#2166ac", "#4393c3", "#f7f7f7", "#d6604d", "#b2182b"))
  annocol <- get_palette("uchicago", 9)
  chromcol <- list(chromosome = c("12" = annocol[6], "other" = annocol[5]))
  
  annocolor <- list(Variant = c("IGHV-UM" = annocol[3], "IGHV-M" = annocol[5], "trisomy12" = annocol[7], "both" = annocol[9]))
  names(annocolor$Variant) <- c("IGHV-UM", "IGHV-M", "trisomy12", "both")
  mutationStatus <- data.frame(mutStatus$mut)
  rownames(mutationStatus) <- mutStatus$PatID
  colnames(mutationStatus) <- "Variant"

  #Column annotation
  ha_col = HeatmapAnnotation(df = mutationStatus, col = annocolor, annotation_height = unit(1.8, "cm"),
                             annotation_legend_param = list(title_gp = gpar(fontsize = 23), 
                                                            labels_gp = gpar(fontsize = 18),  
                                                            grid_height = unit(0.7, "cm"), 
                                                            grid_width = unit(0.3, "cm"), 
                                                            gap = unit(15, "mm")))

  #rowcluster
  geneExpression_dist <- dist(geneExpression_new)
  rowcluster = hclust(geneExpression_dist, method = "ward.D2")

  #heatmap
  h1 <- Heatmap(geneExpression_new, 
                col = colors,
                column_title = paste0("Gene interactions:", "IGHV", "-", "trisomy12"), 
                column_title_gp = gpar(fontsize = 23, fontface = "bold"), 
                heatmap_legend_param = list(title = "Expr", 
                                            title_gp = gpar(fontsize = 23), 
                                            grid_height = unit(0.7, "cm"), 
                                            grid_width = unit(0.3, "cm"), 
                                            gap = unit(15, "mm"), 
                                            labels_gp = gpar(fontsize = 18), 
                                            labels = c(-4, -1, 0, 1, 4)), 
                row_dend_width = unit(0.7, "cm"), 
                show_row_dend = T, 
                show_column_names =FALSE ,
                top_annotation = ha_col,
                show_row_names = FALSE, 
                show_column_dend = FALSE, 
                row_title_gp = gpar(fontsize = 0),
                cluster_columns = FALSE, 
                cluster_rows = rowcluster, 
                split = 4, gap = unit(0.2,"cm"), 
                column_order = mutStatus$PatID)

  #chromosome annotation
  #chromosome distribution
chrom <- as.data.frame(rowData(RNAnorm[sig_Genes,])$chromosome)
rownames(chrom) <- sig_Genes
colnames(chrom ) <- "chromosome"
#select for chromosome12
chrom$chromosome <- ifelse(chrom$chromosome == 12, 12, "other")
  
  ha_chrom = rowAnnotation(df = chrom, 
                           col = chromcol, 
                           annotation_width = unit(0.8, "cm"), 
                           annotation_legend_param = list(ncol = 2, 
                                                          title_gp = gpar(fontsize = 23), 
                                                          labels_gp = gpar(fontsize = 18),  
                                                          grid_height = unit(0.7, "cm"), 
                                                          grid_width = unit(0.3, "cm")))
  
   
  #Annotate most significant genes
  top50 <-  rownames(resTab[which(abs(resTab$stat) > 5 ),])
  int_genes <- rowData(RNAnorm[top50,])$symbol

subset <- as.data.frame(rowData(RNAnorm[sig_Genes,]))
subset <- subset[-which(subset$symbol %in% ""),]
subset <- subset[-which(subset$symbol %in% NA),]

subset <- subset[which(subset$symbol %in% int_genes),]
rownames(subset) <- subset$symbol
geneIDs <- which(rownames(geneExpression_new) %in% rownames(subset))
labels <- rownames(geneExpression_new)[geneIDs]
ha_genes <- rowAnnotation(link = row_anno_link(at = geneIDs, labels = labels, 
                                               labels_gp = gpar(fontsize = 15)), 
                          width = unit(2.5, "cm"))

Warning: anno_link() is deprecated, please use anno_mark() instead.

  #svg(filename="~/git/figures_thesis/gene_expr/epistatsisTri12IGHV.svg", width=30, height=50)
  #pdf(file="/home/almut/Dokumente/git/Transcriptome_CLL/paper/figures/epistasis_Deseq.pdf", width=35, height=45)
p1 <- draw(h1 + ha_genes ) 

Past versions of “heatmap of interaction ighv-1.png”

  #dev.off()
  #draw(h1 + ha_chrom + ha_genes)
#saveRDS(p1, file = paste0(output_dir, "/figures/r_objects/epistasis/epistasis_heatmap.rds"))
IGHVTri12 <- list("sig_Genes" = sig_Genes, "geneExp" = geneExpression_new, "h1"= h1)
p1

Heatmap filtered of interacting genes

 resSig <- subset(resOrdered, padj < 0.01 & abs(stat) > 4)
  #filter by variant
  sig_Genes <- rownames(resSig)
  
  #gene expression data
  geneExpression = assay(RNAnorm)[sig_Genes,]
  rownames(geneExpression) <- rowData(RNAnorm)$symbol[which(rownames(RNAnorm) %in% sig_Genes)]

  #scale and censor
  geneExpression_new <- log2(geneExpression)
  geneExpression_new<- t(scale(t(geneExpression_new)))
  geneExpression_new[geneExpression_new > 4] <- 4
  geneExpression_new[geneExpression_new < -4] <- -4

  mutnames <- c("IGHV-UM", "IGHV-M", "trisomy12", "both")
  
  mutStatus <- data.frame(colData(RNAnorm)) %>% mutate(IGHVnew = ifelse(IGHV %in% "M", 1, 0)) %>% 
    dplyr::select(-IGHV) %>% mutate(IGHV = IGHVnew) %>% 
    dplyr::select_("PatID", "IGHV", "trisomy12") %>% 
    mutate_("namA" = "IGHV", "namB" = "trisomy12") %>% 
    mutate(naA = as.numeric(as.character(namA))) %>% 
    mutate(naB = as.numeric(as.character(namB))) %>% 
    mutate(mut = factor(mutnames[1 + naA + 2 * naB], levels = mutnames)) %>% arrange(mut)

  geneExpression_new <- geneExpression_new[,mutStatus$PatID]


  #colors
  #colors <- colorRampPalette( rev(brewer.pal(11,"RdBu")) )(255)
  colors = colorRamp2(c(-4, -1, 0, 1, 4), c("#2166ac", "#4393c3", "#f7f7f7", "#d6604d", "#b2182b"))
  annocol <- get_palette("uchicago", 9)
  chromcol <- list(chromosome = c("12" = annocol[6], "other" = annocol[5]))
  
  annocolor <- list(Variant = c("IGHV-UM" = annocol[3], "IGHV-M" = annocol[5], "trisomy12" = annocol[7], "both" = annocol[9]))
  names(annocolor$Variant) <- c("IGHV-UM", "IGHV-M", "trisomy12", "both")
  mutationStatus <- data.frame(mutStatus$mut)
  rownames(mutationStatus) <- mutStatus$PatID
  colnames(mutationStatus) <- "Variant"

  #Column annotation
  ha_col = HeatmapAnnotation(df = mutationStatus, col = annocolor, simple_anno_size = unit(0.9, "cm"),
                             annotation_legend_param = list(title_gp = gpar(fontsize = 23), 
                                                            labels_gp = gpar(fontsize = 18),  
                                                            grid_height = unit(1, "cm"), 
                                                            grid_width = unit(0.3, "cm"), 
                                                            gap = unit(15, "mm")))

  #rowcluster
  geneExpression_dist <- dist(geneExpression_new)
  rowcluster = hclust(geneExpression_dist, method = "ward.D2")

  #heatmap
  h1 <- Heatmap(geneExpression_new, 
                col = colors,
                column_title = paste0("Gene interactions:", "IGHV", "-", "trisomy12"), 
                column_title_gp = gpar(fontsize = 23, fontface = "bold"), 
                heatmap_legend_param = list(title = "Expr", 
                                            title_gp = gpar(fontsize = 23), 
                                            grid_height = unit(0.7, "cm"), 
                                            grid_width = unit(0.3, "cm"), 
                                            gap = unit(15, "mm"), 
                                            labels_gp = gpar(fontsize = 18), 
                                            labels = c(-4, -1, 0, 1, 4)), 
                row_dend_width = unit(0.7, "cm"), 
                show_row_dend = T, 
                show_column_names =FALSE ,
                top_annotation = ha_col,
                show_row_names = FALSE, 
                show_column_dend = FALSE, 
                row_title_gp = gpar(fontsize = 0),
                cluster_columns = FALSE, 
                cluster_rows = rowcluster, 
                split = 4, gap = unit(0.2,"cm"), 
                column_order = mutStatus$PatID)


  
  #svg(filename="~/git/figures_thesis/gene_expr/epistatsisTri12IGHV.svg", width=30, height=50)
  #pdf(file="/home/almut/Dokumente/git/Transcriptome_CLL/paper/figures/epistasis_Deseq.pdf", width=35, height=45)
p1 <- draw(h1) 

  #dev.off()
  #draw(h1 + ha_chrom + ha_genes)
saveRDS(p1, file = paste0(output_dir, "/figures/r_objects/epistasis/epistasis_heatmap.rds"))

Genewise count distribution

cluster <- c("buffering", "supression", "synergy", "inversion")
gene_by_cat <- lapply(1:4, function(clusternr){
  genes <- IGHVTri12$sig_Genes[row_order(IGHVTri12$h1)[[clusternr]]]
  gene_symbol <- rowData(ddsCLL)[genes, "symbol"]
}) %>% set_names(cluster)

#function to create stripchart plots for specific genes
gene_count <- function(gene_nam){
  gene_cat <- names(gene_by_cat) %>% map(function(cat){
    cat_gene <- "none"
    cat_gene <- ifelse(any(gene_nam %in% gene_by_cat[[cat]]), cat, cat_gene) 
    }) %>% unlist()
  gene_cat <- gene_cat[!gene_cat %in% "none"]
  gene_cat <- ifelse(gene_nam %in% "EBF1", "synergy", gene_cat)
  gene_cat <- ifelse(gene_nam %in% "FGF2", "suppression", gene_cat)
  ddsCLL$IGHV_tri12 <- mutationStatus[colData(ddsCLL)$PatID,]
  geneEnsID <- rownames(ddsCLL)[which(rowData(ddsCLL)$symbol %in% gene_nam)]
  gc <- plotCounts(ddsCLL, gene = geneEnsID, intgroup = "IGHV_tri12", returnData=TRUE)
  p <- ggboxplot(gc, x = "IGHV_tri12", y = "count",
          color = "IGHV_tri12",
          size = 1.2,
          palette = c(annocol[3], annocol[5], annocol[7], annocol[9]),
          add = "jitter",
          outlier.shape = NA,
          add.params = list(size = 2.5),
          yscale = "log10",
          title = paste0(gene_nam, ": ", gene_cat),
          font.x = 20, font.y = 20, font.legend = 20, 
          ylab = "normalized counts") + font("xy.text", size = 20) + font("title", size = 20, face = "bold")
   #ggsave(file=paste0("/home/almut/Dokumente/git/Transcriptome_CLL/paper/figures/epi_genes/genetic_interaction_", gene_nam, ".svg"),      plot=p, width=7, height=5)
  saveRDS(p, file = paste0(output_dir, "/figures/r_objects/epistasis/de_genes/", gene_nam, ".rds"))
  p
}


#interesting genes
diff <- resTab[which(abs(resTab$stat) > 6 ),]
geneList <- diff$symbol
geneList <- geneList[-which(geneList %in% "")]
geneList <- c(geneList, "LEF1", "TIMELESS", "CHAD", "BCL2A1", "EML6", "PPP1R14A", "EPHB6", "GEN1", "EBF1", 
              "EBF4", "SLC4A8", "CAMK2N1", "FGF2")


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

[[12]]

[[13]]

[[14]]

[[15]]

[[16]]

[[17]]

[[18]]

[[19]]

[[20]]

[[21]]

[[22]]

[[23]]

[[24]]

[[25]]

[[26]]

[[27]]

[[28]]

Clusterwise count distribution

cluster_size <- lapply(c(1:4), function(clusternr){
  size <- length(row_order(IGHVTri12$h1)[[clusternr]])
  name <- cluster[clusternr]
  clust <- c(name, as.numeric(size))
}) %>% do.call(rbind,.) 

cluster_size <- as.data.frame(cluster_size)
colnames(cluster_size) <- c("name", "size")
cluster_size$size <- as.numeric(as.character(cluster_size$size))
distr <- ggbarplot(cluster_size, "name", "size",
   fill = "name", color = "name",
   palette = "jco",
   ylab = "Number of genes",
   xlab = "cluster")
 
 #ggsave(file="/home/almut/Dokumente/masterarbeit/workinprogress/distrib_ofEpiTypes.svg", plot=distr, width=8, height=5)
distr

Past versions of “number of genes within cluster-1.png”

Gene set enrichment analysis

Gene sets

#convert names
convert_names <- function(nam_vec){
  nam_vec <- gsub("HALLMARK_", "", nam_vec)
  nam_vec<- gsub("_", " ", nam_vec)
  nam_vec <- tolower(nam_vec)
  nam_vec <- gsub("tnfa signaling via nfkb", "TNFA signaling via NFKB", nam_vec)
  nam_vec <- gsub("nfkb", "NFKB", nam_vec)
  nam_vec <- gsub("myc", "Myc", nam_vec)
  nam_vec <- gsub("il2 stat5", "IL2 STAT5", nam_vec)
  nam_vec <- gsub("uv response", "UV response", nam_vec)
  nam_vec <- gsub("e2f targets", "E2F targets", nam_vec)
  nam_vec
}

#load gene set collection
#Hallmark
gsc <- loadGSC("/home/almut/Dokumente/masterarbeit/data/h.all.v6.0.symbols.gmt", type="gmt")

#names(gsc$gsc) <- convert_names(names(gsc$gsc))
#gsc$addInfo[,1] <- convert_names(gsc$addInfo[,1])

#Kegg
gsc_Kegg <- loadGSC("/home/almut/Dokumente/masterarbeit/data/c2.cp.kegg.v6.0.symbols.gmt", type="gmt")


diff_res <- resOrderedTab
diff_res$ID <- rownames(diff_res)

#clusterProfiler
diff_res <- diff_res[-which(diff_res$symbol %in% c("", NA)),]
gene_list <- diff_res$stat %>% set_names(diff_res$symbol)
dup <- names(gene_list)[duplicated(names(gene_list))]
gene_list <- gene_list[-which(names(gene_list) %in% dup)]
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"), : 9.63% 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.07% 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...

leading edge analysis...

done...

em <- enricher(de_gene, TERM2GENE = m_t2g, pvalueCutoff = 0.2)

#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

#em2@result$ID <- convert_names(em2@result$ID)
#em2@result$Description <- convert_names(em2@result$Description)
barplot(kk, showCategory=5)

barplot(em, showCategory=5)

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

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

dot1

clusterProfiler::dotplot(em, showCategory=10) + ggtitle("Enrichment for epistatsis trisomy12 and IGHV")

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

clusterProfiler::dotplot(kk2, showCategory=10) + ggtitle("GSEA for epistatsis trisomy12 and IGHV")

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

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

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

dot2

ridgeplot(em2)

Picking joint bandwidth of 0.282

ridgeplot(kk2)

Picking joint bandwidth of 0.283

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/epistasis/enrich_dot_hm.rds"))
saveRDS(dot2, file = paste0(output_dir, "/figures/r_objects/epistasis/enrich_dot2.rds"))

network plot

# Networks Hallmark
#  em2_sub <- em2
#  em2_sub@result <- em2@result[which(em2@result$Description %in% c("TNFA signaling via NFKB",
#                                                                   "IL2 STAT5 signaling",
#                                                                     "Myc targets v2")),]
# 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("Cytokine-cytokine receptor interaction",
                                                                    "Hematopoietic cell lineage",
                                                                    "Antigen processing and presentation",
                                                                    "Apoptosis")),]

pnet_kegg <- cnetplot(kk2_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.

pnet_kegg

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

cluster_gsea <- function(clusternr){
  gene_in <- IGHVTri12$sig_Genes[row_order(IGHVTri12$h1)[[clusternr]]] 
  gene_id <- bitr(gene_in, fromType = "ENSEMBL",
        toType = c("ENTREZID", "SYMBOL"),
        OrgDb = org.Hs.eg.db)

  names(gene_id) <- c("ID", "ENTREZID", "symbol")

  #Hallmark
  em <- enricher(gene_id$symbol, TERM2GENE = m_t2g, pvalueCutoff = 0.5, qvalueCutoff = 0.5)

  #Kegg
  kk <- enrichKEGG(gene_id$ENTREZID,
                 organism     = 'hsa',
                 pvalueCutoff = 0.5,
                 qvalueCutoff = 0.5)
  plot(barplot(em, showCategory=5))
  plot(barplot(kk, showCategory=5))
}

lapply(1:4, cluster_gsea)

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

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

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

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

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

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

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

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

[[1]]

[[2]]

[[3]]

[[4]]

epitype_gsea <- function(clusternr){
   gene_in <- IGHVTri12$sig_Genes[row_order(IGHVTri12$h1)[[clusternr]]] 
   gene_out <- IGHVTri12$sig_Genes[!IGHVTri12$sig_Genes %in% gene_in]
   diff_res <- resOrderedTab
   diff_res$ID <- rownames(diff_res)

   #clusterProfiler
   #filter genes
   diff_res <- diff_res[-which(diff_res$symbol %in% c("", NA)),]
   diff_res <- diff_res[!duplicated(diff_res$symbol),]
   diff_res <- diff_res[!rownames(diff_res) %in% gene_out,]
   gene_list <- diff_res$stat %>% set_names(diff_res$symbol)
   gene_list <- sort(gene_list, decreasing = TRUE)
   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)

    #Hallmark
    em2 <- GSEA(gene_list, TERM2GENE = m_t2g, pvalueCutoff = 0.1)
    #em2@result$ID <- convert_names(em2@result$ID)
    #em2@result$Description <- convert_names(em2@result$Description)
    #Kegg
    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')
    dot1 <- clusterProfiler::dotplot(em2, showCategory=12) + ggtitle(paste0("GSEA epistasis ", cluster[clusternr])) +
      theme_pubr() +
      theme(legend.position="right") + 
      theme(plot.title = element_text(face = "bold")) 
    plot(dot1)
    saveRDS(dot1, file = paste0(output_dir, "/figures/r_objects/epistasis/enrich_dot_", cluster[clusternr],".rds"))
    dot2 <- clusterProfiler::dotplot(kk2, showCategory=10) + ggtitle(paste0("GSEA epistasis ", cluster[clusternr])) +
    theme_pubr() +
    theme(legend.position="right") +
    theme(plot.title = element_text(face = "bold"))
    plot(dot2)

    plot(ridgeplot(em2))
    plot(ridgeplot(kk2))
    em2
}

em_list <- lapply(1:4, epitype_gsea) %>% set_names(cluster)

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

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

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

preparing geneSet collections...

GSEA analysis...

leading edge analysis...

done...

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

Picking joint bandwidth of 0.217

Picking joint bandwidth of 0.216

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

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

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

preparing geneSet collections...

GSEA analysis...

leading edge analysis...

done...

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

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

Picking joint bandwidth of 0.218

Picking joint bandwidth of 0.217

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

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

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

preparing geneSet collections...

GSEA analysis...

leading edge analysis...

done...

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

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

Picking joint bandwidth of 0.217

Picking joint bandwidth of 0.241

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

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

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

preparing geneSet collections...

GSEA analysis...

leading edge analysis...

done...

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

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

Picking joint bandwidth of 0.258

Picking joint bandwidth of 0.267

#get unique pathways per type
sig_pw <- lapply(1:4, function(cluster_nr){
  pw_unique <- em_list[[cluster[cluster_nr]]]@result %>% filter(p.adjust < 0.05) %>% dplyr::select(ID)
  pw_unique$ID
}) %>% set_names(cluster)

dup_pw <- unlist(sig_pw)[duplicated(unlist(sig_pw))]
all_pw <- Reduce(intersect, sig_pw)
unique_pw <- sig_pw %>% map(function(pw){pw <- pw[!pw %in% all_pw]})
unique_pw 

$buffering
[1] "HALLMARK_E2F_TARGETS" "HALLMARK_APOPTOSIS"  

$supression
[1] "HALLMARK_APOPTOSIS"

$synergy
[1] "HALLMARK_G2M_CHECKPOINT"          "HALLMARK_INFLAMMATORY_RESPONSE"  
[3] "HALLMARK_E2F_TARGETS"             "HALLMARK_IL6_JAK_STAT3_SIGNALING"

$inversion
[1] "HALLMARK_INFLAMMATORY_RESPONSE"    
[2] "HALLMARK_INTERFERON_GAMMA_RESPONSE"
[3] "HALLMARK_IL6_JAK_STAT3_SIGNALING"  
[4] "HALLMARK_APOPTOSIS"                
[5] "HALLMARK_G2M_CHECKPOINT"           
[6] "HALLMARK_E2F_TARGETS"