MethylAnalysis_bulk_v3.R

MethylAnalysis_bulk_v3 <- function(methcall.folder = NULL, # input folder with .cov files (bismark)
                           outfolder = NULL, # output folder (if not exist it will be created)
                           proj.id = "myproject", # project id to prefix in putput files
                           samplesheet = NULL,
                           sheetnum = 1,
                           design.var = "Disease", # variable to subset samplesheet
                           case.var = "Case", # 1 in treatment vector
                           control.var = "Control", # 0 in treatment vector
                           idcol = "SampleID",
                           mincoverage = 10,
                           lowperc = 5,
                           lowcount = 10,
                           assem="hg38",
                           do.subset = T,
                           chr.subset = "chr9",
                           start.subset = 136668000,
                           end.subset = 136671000,
                           pipeline.meth = "bismarkCoverage",
                           plot.covariates = c("Condition","Batch","MRD",
                                               "Torelapse","Chemorefractory",
                                               "Sex","CytoA","Tissue"),
                           idstoremove = NULL,
                           delta.meth = 10,
                           plot.categorical.vars = c("Condition","Batch","MRD","Torelapse","Chemorefractory","Sex","CytoA"),
                           plot.continuous.vars = c("miR-126","egfl7"),
                           plot.height = 9,
                           plot.width = 12,
                           cellw = 15, 
                           cellh = 15,
                           fontnumsize = 5,
                           fontsize = 8
                           ){
  
  # load libraries
  
  require(methylKit)
  require(ggplot2)
  require(reshape2)
  require(plyr)
  require(ggrepel)
  library(GenomicRanges)
  library(openxlsx)
  library(pheatmap)
  
  # check vars
  
  if(is.null(methcall.folder)){
    stop("Please set methcall.folder")
  }
  
  if(is.null(outfolder)){
    stop("Please set outfolder")
  }
  
  if(is.null(samplesheet)){
    stop("Please set samplesheet")
  }
  
  # Setup variables
  
  infolder = paste0(methcall.folder, "/")
  dir.create(infolder, showWarnings = F)
  
  outfolder = paste0(outfolder, "/")
  dir.create(outfolder, showWarnings = F)
  
  qcfolder <- paste0(outfolder, "QC/")
  dir.create(qcfolder, showWarnings = F)
  
  
  # reading sample sheet with metadata
  ssheet <- read.xlsx(samplesheet, sheet = sheetnum)
  
  if(!is.null(idstoremove)){
    ssheet <- subset.data.frame(x = ssheet, subset = !ssheet[[idcol]] %in% idstoremove)
  }

  # defining design vector according to variable
  
  ssheet <- subset.data.frame(x = ssheet, subset = ssheet[[design.var]] %in% c(case.var, control.var))
  
  d.vec <- ssheet[[design.var]]
  d.vector <- ifelse(d.vec == case.var, 1, 0)
  
  # initialzing coverage .cov files list
  covs <- list()
  sampleids <- as.list(ssheet[[idcol]])
  names(sampleids) <- ssheet[[idcol]]
  
  for (i in ssheet[[idcol]]) {
    covs[[i]] <- paste0(infolder, "/", i, ".bismark.cov")
  }
  
  saveRDS(object = covs, file = "covs_object.rds")
  
  # creating object 
  
  myobj <- methRead(location = covs,
                    sample.id = sampleids,
                    assembly = assem,
                    pipeline = pipeline.meth,
                    treatment = d.vector,
                    context="CpG", 
                    mincov = mincoverage)
  
  saveRDS(myobj, file = "Myobject.rds")
  
  names(myobj) <- ssheet[[idcol]]
  
  saveRDS(myobj, file = "Myobject_with_names.rds")

  # subsetting if declared
  
  if(isTRUE(do.subset)){
    my.win = GRanges(seqnames = chr.subset, ranges = IRanges(start = start.subset, end = end.subset))
    myobj <- selectByOverlap(myobj,my.win)
  }

  
  saveRDS(myobj, file = "Myobject_with_names_aftersubset.rds")
  
  # filtering on minimum coverage 
  myobj <- filterByCoverage(methylObj = myobj, lo.count=lowcount, lo.perc = lowperc)
  
  # Normalization
  myobj <- normalizeCoverage(obj = myobj)
  
 # saveRDS(object = myobj, file = "Initial_object.rds")
  
  # Calculate basic stats and PCs

  metricsfolder <- paste0(qcfolder, "Metrics/")
  dir.create(path = metricsfolder, showWarnings = F)
    
  for (id in ssheet[[idcol]]) {
    
    png(filename = paste0(metricsfolder,proj.id,"_CpG_pct_methylation_sample_", id, ".png"), 
        width = 9, height = 6, units = "in", res = 96)
    print(getMethylationStats(myobj[[id]],plot=TRUE,both.strands=FALSE))
    dev.off()
    png(filename = paste0(metricsfolder, proj.id,"_Coverage_stats_sample_", id, ".png"), 
        width = 9, height = 6, units = "in", res = 96)
    print(getCoverageStats(myobj[[id]],plot=TRUE,both.strands=FALSE))
    dev.off()
  }
  
  # create meth obj 
  
  #meth <- unite(object = myobj, destrand=FALSE)
  meth <- unite(object = myobj, destrand=FALSE) # for debugging
  saveRDS(meth, "savemeth.tmp.rds")
  # Perform correlation
  
  sink(paste0(qcfolder, proj.id, "_Correlations.txt"))
  getCorrelation(meth,plot=FALSE)
  sink()
  
  if(length(ssheet[[idcol]]) < 15){
    png(filename = paste0(qcfolder,proj.id, "_Correlations_pearson_pairwise.png"), 
        width = 9, height = 6, units = "in", res = 96)
    print(getCorrelation(meth,plot=TRUE))
    dev.off()
  }
  
  png(filename = paste0(qcfolder, proj.id, "_Clustering.png"), 
      width = 9, height = 6, units = "in", res = 96)
  clusterSamples(meth, dist="euclidean", plot=TRUE, method = "ward.D2")
  dev.off()
  
  # Re-plotting PCs (custom chart)
  
  # compute PCs and store in object
  
  pca_compt <- PCASamples(meth, obj.return = T, screeplot = F)
  
  # extract PCs components
  
  pcafolder <- paste0(qcfolder, "PCA/")
  dir.create(path = pcafolder, showWarnings = F)
  
  pca_pc1_2 <- as.data.frame(x = pca_compt$x[,1:2])
  
  for(myvars in plot.covariates){
    pca_pc1_2$condition <- as.factor(ssheet[[myvars]])
    png(filename = paste0(pcafolder, proj.id, "_PCA_",myvars,".png"), 
        width = 9, height = 6, units = "in",res=96)
    print(ggplot(data = pca_pc1_2, 
                 mapping = aes(x = PC1, y=PC2, col=condition, label=rownames(pca_pc1_2))) + 
            geom_point(size=3) + geom_text_repel(size=3) + ggtitle(label = "Principal component analysis", subtitle = myvars) + 
            theme(plot.title = element_text(size = 16, face = "bold", hjust = 0.5)) +
            theme(plot.subtitle=element_text(size=12, hjust=0.5, face="italic", color="black")) +
            theme(axis.title = element_text(size=12, hjust=0.5, face="bold", color="black")) +
            theme(legend.text = element_text(size=8, hjust=0.5)) +
            theme(legend.title = element_blank()) +
            theme(axis.text = element_text(size=12, hjust=0.5, color="black")))
    dev.off()
    
  }
  
  # retrieve and store % of methylation 
  
  perc.meth <- percMethylation(meth)
  
  saveRDS(perc.meth,"pctmethly.rds")
  
  base::rownames(perc.meth) <- paste0(meth$chr, "_", meth$start)
  
  # Perform diff methylation
  
  myDiff=calculateDiffMeth(meth)
  
  write.table(myDiff,paste0(outfolder,proj.id,"_DiffMeth_single_CpG.txt"), row.names = F)
  difftest <- read.table(paste0(outfolder,proj.id,"_DiffMeth_single_CpG.txt"), header = T)

  difftest$comparison <- proj.id
  difftest$qvalue_r <- as.character(cut(x = difftest$qvalue,
                                        breaks = c(-1, 1e-100, 1e-10, 1e-02, 1),
                                        labels = c("***","**","*","ns")))
  
  myindex <- abs(difftest$meth.diff) < delta.meth
  difftest$meth.diff <- abs(difftest$meth.diff)
  difftest$qvalue_r[myindex] <- "ns"
  write.table(difftest,paste0(outfolder,proj.id,"_DiffMeth_single_CpG.txt"), row.names = F)
  

  # Adding color list and annotations

  library(RColorBrewer)

  color_list <- list()
  annrows <- NULL
  anncols <- subset.data.frame(difftest, select = c("qvalue_r","meth.diff"))
  base::rownames(anncols) <- difftest$start
  
  rows.annot.vars.cat = plot.categorical.vars
  rows.annot.vars.con = plot.continuous.vars

  # if(!is.null(rows.annot.vars.cat) | !is.null(rows.annot.vars.con)){
  #   rownames(ssheet) <- ssheet$SampleID
  #   annrows <- subset.data.frame(x = ssheet, select = c(rows.annot.vars.cat, rows.annot.vars.con))
  #   concolors <- RColorBrewer::brewer.pal(n = 9, name = "Set1")
  #   catcolors <- NULL
  #   for (varcon in 1:length(rows.annot.vars.con)) {
  #     color_list[[rows.annot.vars.con[varcon]]] <- colorRampPalette(c("lightgrey", concolors[varcon]))(10)
  #   }
  #   for (varcat in 1:length(rows.annot.vars.cat)) {
  #     ncolor <- 1
  #     for (val in unique(ssheet[[rows.annot.vars.cat[varcat]]])[order(unique(ssheet[[rows.annot.vars.cat[varcat]]]))]){
  #       if(length(unique(ssheet[[rows.annot.vars.cat[varcat]]])) > 9){
  #         catcolors <- colorRampPalette(RColorBrewer::brewer.pal(n = 9, name = "Set3"))(length(unique(ssheet[[rows.annot.vars.cat[varcat]]])))
  #       }
  #       else{
  #         catcolors <- RColorBrewer::brewer.pal(n = 9, name = "Set3")
  #       }
  #       color_list[[rows.annot.vars.cat[varcat]]][[val]] <- catcolors[ncolor]
  #       ncolor <- ncolor + 1
  #     }
  #   }
  # }
  
  color_list[["qvalue_r"]] = c("*" = "#6497b1", "**" = "#03396c", "***" = "#011f4b", ns= "#c4cacf")
  color_list[["meth.diff"]] = colorRampPalette(brewer.pal(n = 11, "Reds"))(100)
  color_list[["miR-126"]] <- brewer.pal(n = 9, "PuRd")
  
  # Heatmap
  
  pctmeth_matrix <- t(perc.meth)
  base::colnames(pctmeth_matrix) <- gsub(x = base::colnames(pctmeth_matrix), pattern = "chr[0-9]+_", replacement = "")
  
  pheatmap(mat = pctmeth_matrix, main = gsub(x = proj.id, pattern = "_", replacement = " "), 
           filename = paste0(outfolder, proj.id,"_CpG_percent_methylation_matrix_pheatmap.pdf"), width = plot.width, height = plot.height,
           na_col = "black", 
           cluster_cols = FALSE, 
           cluster_rows = TRUE, 
           annotation_row = annrows, 
           cellwidth = cellw, 
           cellheight = cellh, 
           display_numbers = T, 
           fontsize = fontsize,
           fontsize_number = fontnumsize, 
           number_format = "%.0f",
           #border_color = "#CBBEB5", 
           annotation_col = anncols, 
           #labels_row = lrow,
           #labels_col = lcol,
           #gaps_row = gaps.row,
           gaps_col = c(8),
           annotation_colors = color_list,
           color = c("#F5F5F5","#EEEEEE","#CCCCCC","#999999", "#666666","#333333","#000000"), breaks = c(0,10,20,30,50,70,90,100)
  )
  
  saveRDS(object = pctmeth_matrix, paste0(outfolder, proj.id,"_CpG_percent_methylation_matrix_pheatmap.rds"))
  saveRDS(object = list(anncol = anncols, annrow = annrows, anncolors = color_list), paste0(outfolder, proj.id,"_annotations_matrix_pheatmap.rds"))
  saveRDS(object = difftest, paste0(outfolder, proj.id,"_CpG_differential_methylation.rds"))
  saveRDS(object = perc.meth, paste0(outfolder, proj.id,"_CpG_percent_methylation.rds"))
  
  write.table(x = perc.meth, file = paste0(outfolder, proj.id,"_CpG_percent_methylation.txt"))
  write.table(x = difftest, file = paste0(outfolder, proj.id,"_CpG_differential_methylation.txt"))

  saveRDS(myobj, file = paste0(outfolder,proj.id,"_methylkit.rds"))
  saveRDS(myobj, file = paste0(outfolder,proj.id,"_methylkit_meth.rds"))

}