Methylation groups

Analyse genes associated with methylation groups

libraries

library(tidyverse)
library(ggplot2)
library(DESeq2)
library(ggpubr)
library(ComplexHeatmap)
library(RColorBrewer)
library(circlize)
library(here)
library(piano)

Data preprocessing

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

#arrange columns
mutStatus <- data.frame(colData(ddsCLL)) %>% arrange(Methylation) %>% filter(!is.na(Methylation)) %>% mutate("IP" = ifelse(Methylation %in% c("LP", "HP"), 0, 1))
colnames(ddsCLL) <-colData(ddsCLL)$PatID
ddsCLL <- ddsCLL[, mutStatus$PatID]
colData(ddsCLL)$IP <- mutStatus$IP
table(colData(ddsCLL)$IP)

  0   1 
141  32 

Normalize

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

Exploratory data analysis

Clustering on the most variable genes already split methylation groups. How do methylation groups affect highly variable genes? Can we distinguish all 3 groups?

PCA

#Plot PCA
exprMat <- assay(RNAnorm)

#top 5000 most variant genes
sds <- rowSds(exprMat)
exprMat <- exprMat[order(sds, decreasing = T)[1:150],]


#Calculate PCA
pcaRes <- prcomp(t(exprMat), scale =T)
varExp <- (pcaRes$sdev^2 / sum(pcaRes$sdev^2)) * 100
pcaTab <- data.frame(pcaRes$x[,c(1:10)])
names(varExp) <- colnames(pcaRes$x)

#add background information
pcaTab <- cbind(pcaTab, data.frame(colData(RNAnorm)))


#plot PCA and color samples based on annotations
annocol <- get_palette("jco", 10)
  
p <- ggscatter(pcaTab, x = "PC1", y = "PC2", color = "Methylation", palette = c( annocol[7], annocol[5], annocol[6]),
  ylab = sprintf("PC2 (%2.1f%%)",varExp[2]), xlab = sprintf("PC1 (%2.1f%%)",varExp[1]), legend = "right", main = "PCA Methylation groups") + coord_fixed()

p

Past versions of pca-1.png

#ggsave(file=paste0(figure_dir, "/pca_Meth_IP_all_top150.svg"), plot=p, width=5, height=5)

Gene expression

###Deseq
ddsCLL <- estimateSizeFactors(ddsCLL)

# deseq2 function
diff <- function(cond){
  ddsCLL_new <- ddsCLL[,!is.na(colData(ddsCLL)[,cond])]
  design(ddsCLL_new) <- as.formula(paste("~ ", paste(cond)))
  rnaRaw <- DESeq(ddsCLL_new, betaPrior = FALSE)
  res <- results(rnaRaw)
  resOrdered <- res[order(res$pvalue),]
}

#diff_meth <- diff("IP")

#saveRDS(diff_meth, file= paste0(output_dir,"/diff_meth_IP_vs_all.rds"))

diff_meth <- readRDS(paste0(output_dir,"/diff_meth_IP_vs_all.rds"))

dataTab <- data.frame(diff_meth)
dataTab$ID <- rownames(dataTab)

#filter using pvalues
dataTab <- dataTab %>%
    arrange(padj) %>%
    mutate(Symbol = rowData(ddsCLL[ID,])$symbol)
dataTab <- dataTab[!duplicated(dataTab$Symbol),]
dataTab <- dataTab[!is.na(dataTab$Symbol),]
rownames(dataTab) <- dataTab$ID
write.csv(dataTab, file=paste0(output_dir, "/diff_genes/meth_IP_vs_all_diffGenes.csv"))

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


#Differentially expressed genes
genes <- dataTab %>% filter(padj <= 0.01, abs(stat) > 6)  
exprMat <- assay(RNAnorm)
exprMat<- exprMat[genes$ID,]

#scale gene expression
colnames(exprMat) <- colData(ddsCLL)$PatID
exprMat.new <- log2(exprMat)
exprMat.new <- t(scale(t(exprMat.new)))
exprMat.new[exprMat.new > 4] <- 4
exprMat.new[exprMat.new < -4] <- -4

#colors
colors <- colorRampPalette(rev( brewer.pal(10,"RdBu")) )(20)
annocol <- get_palette("jco", 10)
annocolor <- list(Methylation = c("IP" = annocol[5], "LP" = annocol[6], "HP" =  annocol[7]), IGHV = c("M" = annocol[1], "U" = annocol[2])) 

#Annotation
feature <- as.data.frame(colData(ddsCLL)[,c("Methylation", "IGHV")]) 
colnames(feature) <- c("Methylation", "IGHV") 

#gene symbol as rownames
rownames(exprMat.new) <- rowData(RNAnorm[rownames(exprMat),])$symbol

ha_col <- HeatmapAnnotation(df = feature, col = annocolor, annotation_height = unit(c(rep(1.3, 2)), "cm"), annotation_legend_param = list(title_gp = gpar(fontsize = 40), labels_gp = gpar(fontsize = 35),  grid_height = unit(1.9, "cm"), grid_width = unit(1.9, "cm")))


h1 <- Heatmap(exprMat.new ,                                                    
              km = 2,
              cluster_columns = F,
              clustering_distance_rows = "pearson",
              clustering_method_rows = "ward.D2",
              column_title ="Gene signature methylation groups ",                      
              col = colors,
              column_title_gp = gpar(fontsize = 50, fontface = "bold"), 
              heatmap_legend_param = list(title = "expr", 
                                          title_gp = gpar(fontsize = 40), 
                                          grid_height = unit(1.9, "cm"), 
                                          grid_width = unit(1.9, "cm"), 
                                          gap = unit(2, "cm"), 
                                          labels_gp = gpar(fontsize = 35)), 
              show_row_dend = FALSE, 
              show_column_names = FALSE , 
              show_row_names = FALSE, 
              row_names_gp = gpar(fontsize = 21),
              top_annotation = ha_col)



#Annotate top 50 genes
sub_names <- genes[1:50,"Symbol"]
sub_names <- sub_names[-which(sub_names %in% "")]
geneIDs <- which(rownames(exprMat.new) %in% sub_names)
labels <- rownames(exprMat.new)[geneIDs]
ha_genes <- rowAnnotation(link = row_anno_link(at = geneIDs, labels = labels, labels_gp = gpar(fontsize = 35)), width = unit(9, "cm"))

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

#svg(filename=paste0(figure_dir, "/gene_expr_Methylation_IP_all.svg"), width=30, height=35)
#pdf(file=paste0(figure_dir,"/gene_expr_Methylation_IP_all.pdf"), width=30, height=45)
p1 <- draw( h1 + ha_genes)

#dev.off()

saveRDS(p1, file = paste0(output_dir, "/figures/r_objects/Methylation_IP/meth_ip_heatmap.rds"))

Gene and Sample specific expression - top genes

#function to create stripchart plots for specific genes
gene_count <- function(gene_nam){
  geneEnsID <- rownames(ddsCLL)[which(rowData(ddsCLL)$symbol %in% gene_nam)]
  geneNum <- exprMat[geneEnsID,]
  mutPat <- as.data.frame(colData(ddsCLL)[, c("Methylation")])
  colnames(mutPat) <- c("genotype")
  geneDat <- cbind(mutPat, geneNum)
  colnames(geneDat) <- c("genotype", "counts")
  
  p <- ggstripchart(geneDat, x = "genotype", y = "counts",
          color = "genotype",
          palette = "jco",
          add = "mean_sd",
          title = paste(gene_nam),
          ylab = "normalized counts")
 # ggsave(file=paste0(figure_dir, "/methylation_IP_all/genetic_interaction_", gene_nam, ".svg"), plot=p, width=6, height=5)
  saveRDS(p, file = paste0(output_dir, "/figures/r_objects/Methylation_IP/de_genes/", gene_nam, ".rds"))
  p
}

geneList <- sub_names[1:30]
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]]

Past versions of “single gene counts-11.png”

[[12]]

Past versions of “single gene counts-12.png”

[[13]]

Past versions of “single gene counts-13.png”

[[14]]

Past versions of “single gene counts-14.png”

[[15]]

Past versions of “single gene counts-15.png”

[[16]]

Past versions of “single gene counts-16.png”

[[17]]

Past versions of “single gene counts-17.png”

[[18]]

Past versions of “single gene counts-18.png”

[[19]]

Past versions of “single gene counts-19.png”

[[20]]

Past versions of “single gene counts-20.png”

[[21]]

Past versions of “single gene counts-21.png”

[[22]]

Past versions of “single gene counts-22.png”

[[23]]

Past versions of “single gene counts-23.png”

[[24]]

Past versions of “single gene counts-24.png”

[[25]]

Past versions of “single gene counts-25.png”

[[26]]

Past versions of “single gene counts-26.png”

[[27]]

Past versions of “single gene counts-27.png”

[[28]]

Past versions of “single gene counts-28.png”

[[29]]

Past versions of “single gene counts-29.png”

[[30]]

Past versions of “single gene counts-30.png”

gene_count("SOX11")

Past versions of “single gene counts-31.png”

Gene set enrichment

variant <- "IP"
gmtFile <- loadGSC(paste0(data_dir,"/c2.cp.kegg.v6.0.symbols.gmt"), type="gmt")

diff_res <- dataTab
diff_res  <- diff_res[-which(diff_res$Symbol %in% c("", NA)),]
#get genes and pvalues 
geneNam <- diff_res$Symbol 
pVal <- diff_res$padj
logFold <- diff_res$log2FoldChange
stat <- diff_res$stat
gsTab <- data.frame(gene = geneNam, stat = stat)
gsaTab <- data.frame(row.names = gsTab$gene, stat = gsTab$stat)
res <- runGSA(geneLevelStats = gsaTab,
                      geneSetStat = "gsea",
                     adjMethod = "fdr", gsc=gmtFile,
                     signifMethod = "geneSampling",
                      nPerm = 50000,
                      gsSizeLim=c(1, Inf))

Running gene set analysis:

Checking arguments...done!

*** Please note that running the GSEA-method may take a substantial amount of time! ***

Final gene/gene-set association: 4271 genes and 186 gene sets

  Details:

  Calculating gene set statistics from 4271 out of 21668 gene-level statistics

  Using all 21668 gene-level statistics for significance estimation

  Removed 995 genes from GSC due to lack of matching gene statistics

  Removed 0 gene sets containing no genes after gene removal

  Removed additionally 0 gene sets not matching the size limits

  Loaded additional information for 186 gene sets

Calculating gene set statistics...done!
Calculating gene set significance...done!
Adjusting for multiple testing...done!

Res_up <- arrange(GSAsummaryTable(res), `p adj (dist.dir.up)`)                 
Res_dn <- arrange(GSAsummaryTable(res), `p adj (dist.dir.dn)`)
  

#Plot
resPlot <- Res_dn[, c(1:3,7,8,9)]
colnames(resPlot) <- c("pathway", "gene_number", "stat", "p.adj","genes_up" , "genes_dn")


enrichPlot <- resPlot %>% filter(p.adj < 0.1) %>% mutate(log10Padj = -log10(p.adj)) #%>% mutate(genes = ifelse(gene_number > 5, ">5", "<=5"))
enrichPlot$log10Padj[which(enrichPlot$log10Padj == Inf)] <- 5

p <- ggbarplot(enrichPlot, x = "pathway", y = "log10Padj",
          fill = "gene_number",         
          color = "white",           
          palette =  "gsea",            
          sort.val = "asc",          
          sort.by.groups = FALSE,     
          ylab = "-log10(padj)",
          legend.title = "#diff.genes",
          rotate = TRUE,
          font.x = 20, font.y = 20, font.legend = 20, legend = "right", 
          title = "Methylation groups - Kegg",
          ggtheme = theme_pubr()) +
  font("xy.text", size = 16) + 
  font("title", size = 20, face = "bold")
  
#ggsave(file=paste0(figure_dir, "/GSEA_Meth_IP_vs_all_Kegg.svg"), plot=p, width=14, height=7)

p

Past versions of gsea-1.png

saveRDS(p, file = paste0(output_dir, "/figures/r_objects/Methylation_IP/meth_ip_enrichment.rds"))

Session information

sessionInfo()

R version 3.6.0 (2019-04-26)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.6 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] piano_2.0.2                 here_0.1                   
 [3] circlize_0.4.6              RColorBrewer_1.1-2         
 [5] ComplexHeatmap_2.0.0        ggpubr_0.2                 
 [7] magrittr_1.5                DESeq2_1.24.0              
 [9] SummarizedExperiment_1.14.0 DelayedArray_0.10.0        
[11] BiocParallel_1.18.0         matrixStats_0.54.0         
[13] Biobase_2.44.0              GenomicRanges_1.36.0       
[15] GenomeInfoDb_1.20.0         IRanges_2.18.1             
[17] S4Vectors_0.22.0            BiocGenerics_0.30.0        
[19] forcats_0.4.0               stringr_1.4.0              
[21] dplyr_0.8.1                 purrr_0.3.2                
[23] readr_1.3.1                 tidyr_0.8.3                
[25] tibble_2.1.3                ggplot2_3.1.1              
[27] tidyverse_1.2.1            

loaded via a namespace (and not attached):
  [1] readxl_1.3.1           backports_1.1.4        Hmisc_4.2-0           
  [4] fastmatch_1.1-0        workflowr_1.4.0        plyr_1.8.4            
  [7] igraph_1.2.4.1         lazyeval_0.2.2         shinydashboard_0.7.1  
 [10] splines_3.6.0          digest_0.6.19          htmltools_0.3.6       
 [13] gdata_2.18.0           checkmate_1.9.3        memoise_1.1.0         
 [16] cluster_2.1.0          limma_3.40.2           annotate_1.62.0       
 [19] modelr_0.1.4           colorspace_1.4-1       blob_1.1.1            
 [22] rvest_0.3.4            haven_2.1.0            xfun_0.7              
 [25] crayon_1.3.4           RCurl_1.95-4.12        jsonlite_1.6          
 [28] genefilter_1.66.0      survival_2.44-1.1      glue_1.3.1            
 [31] gtable_0.3.0           zlibbioc_1.30.0        XVector_0.24.0        
 [34] GetoptLong_0.1.7       shape_1.4.4            scales_1.0.0          
 [37] DBI_1.0.0              relations_0.6-8        Rcpp_1.0.1            
 [40] xtable_1.8-4           htmlTable_1.13.1       clue_0.3-57           
 [43] foreign_0.8-71         bit_1.1-14             Formula_1.2-3         
 [46] DT_0.7                 htmlwidgets_1.3        httr_1.4.0            
 [49] fgsea_1.10.0           gplots_3.0.1.1         acepack_1.4.1         
 [52] pkgconfig_2.0.2        XML_3.98-1.20          nnet_7.3-12           
 [55] locfit_1.5-9.1         labeling_0.3           tidyselect_0.2.5      
 [58] rlang_0.3.4            later_0.8.0            AnnotationDbi_1.46.0  
 [61] visNetwork_2.0.7       munsell_0.5.0          cellranger_1.1.0      
 [64] tools_3.6.0            cli_1.1.0              generics_0.0.2        
 [67] RSQLite_2.1.1          broom_0.5.2            evaluate_0.14         
 [70] yaml_2.2.0             knitr_1.23             bit64_0.9-7           
 [73] fs_1.3.1               caTools_1.17.1.2       nlme_3.1-140          
 [76] whisker_0.3-2          mime_0.7               slam_0.1-45           
 [79] xml2_1.2.0             compiler_3.6.0         rstudioapi_0.10       
 [82] png_0.1-7              marray_1.62.0          geneplotter_1.62.0    
 [85] stringi_1.4.3          lattice_0.20-38        Matrix_1.2-17         
 [88] ggsci_2.9              shinyjs_1.0            pillar_1.4.1          
 [91] GlobalOptions_0.1.0    data.table_1.12.2      bitops_1.0-6          
 [94] httpuv_1.5.1           R6_2.4.0               latticeExtra_0.6-28   
 [97] promises_1.0.1         KernSmooth_2.23-15     gridExtra_2.3         
[100] gtools_3.8.1           assertthat_0.2.1       rprojroot_1.3-2       
[103] rjson_0.2.20           withr_2.1.2            GenomeInfoDbData_1.2.1
[106] hms_0.4.2              rpart_4.1-15           rmarkdown_1.13        
[109] git2r_0.25.2           sets_1.0-18            shiny_1.3.2           
[112] lubridate_1.7.4        base64enc_0.1-3