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)

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("none", "IGHV", "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("jco", 10)
  chromcol <- list(chromosome = c("12" = annocol[6], "other" = annocol[5]))
  
  annocolor <- list(Variant = c("none" = annocol[1], annocol[2], annocol[3], "both" = annocol[4]))
  names(annocolor$Variant) <- c("none", "IGHV", "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 = 5, 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)

Genewise count distribution

#function to create stripchart plots for specific genes
gene_count <- function(gene_nam){
  geneEnsID <- rownames(ddsCLL)[which(rowData(ddsCLL)$symbol %in% gene_nam)]
  geneNum <- assay(RNAnorm)[geneEnsID,]
  mutPat <- as.data.frame(mutationStatus[colData(ddsCLL)$PatID,])
  colnames(mutPat) <- c("genotype")
  geneDat <- cbind(mutPat, geneNum)
  colnames(geneDat) <- c("genotype", "counts")
  
  p <- ggstripchart(geneDat, x = "genotype", y = "counts",
          color = "genotype",
          size = 3,
          palette = "jco",
          add = "mean_sd",
          title = paste(gene_nam),
          ylab = "normalized counts",
          font.x = 20, font.y = 20, font.legend = 20) + 
    font("xy.text", size = 20) + font("title", size = 25, 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")


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”

Clusterwise count distribution

cluster <- c("buffering_dn","sup_tri", "buffering_up" ,"buffering_up_strong", "inversion_up")

cluster_size <- lapply(c(1:5), 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”

saveRDS(distr, file = paste0(output_dir, "/figures/r_objects/epistasis/epistasis_gene_distr.rds"))