R/annoByCpG.R

In bioinfo16/IRFinder: Vector integrated region analysis package (IRFinder)

#' @title Annotate the vector integrated site by CpG islands.
#' 
#' @description After load CpG island informaton from a ucsc data file, this function shows CpG island information
#'              related with vector integrated sites. User can get query sequence inserted in CpG site and distribution graph by this function.
#' 
#' @usage annoByCpG(hits, randomSet = NULL, mapTool = 'blast', organism = 'hg19', interval = 5000, range = c(-20000, 20000),
#`        outpath = '~', cpglen = 300, dbPath = paste0(.libPaths()[1], '/IRFinder/extdata'))
#' 
#' @param hits a GR object. This object made from *makeInputs* function.
#' @param randomSet a string vector. Type path to load a random set. 
#'                  If this value is null, random distribution analysis is not executed.
#' @param mapTool a character. Function serve two types of file such as outputs from BLAST and BLAT.
#'                Default is 'blast'. If you want to use BLAT output, use 'blat'.
#' @param organism a single character. This function serves 3 versions of organisms such as hg19, hg38 (Human)
#'                 and galGal6 (Chicken). Default is 'hg19'.
#' @param interval an integer vector. This number means interval number for distribution analysis. Default is 5000.
#' @param range an integer array. It means the range for highlight region of this analysis. Default range is c(-20000, 20000).
#' @param outpath an string vector. Plots are saved in this path. Default value is R home directory.
#' @param cpglen an integer vector. Validate length of CpG sites in islands.
#' @param dbPath a string vector. Directory path of database files.
#' 
#' @return Return a result list constituted by insertion table, distribution table and a GenomicRange object of CpG data.
#' 
#' 
#' @export

annoByCpG = function(hits, randomSet = NULL, mapTool = 'blast', organism = 'hg19', interval = 5000, range = c(-20000, 20000),
                     outpath = '~', cpglen = 300, dbPath = paste0(.libPaths()[1], '/IRFinder/extdata')){
  library(GenomicRanges)
  library(stringr)
  library(grDevices)
  library(regioneR)

  cat('---------- Annotation integrated sites : CpG islands ----------\n')
  cat(paste0('Start time : ', date(), '\n'))

  if(length(which(c('hg19', 'hg38', 'galGal6') %in% organism)) == 0){
    return(cat("You can use hg19/hg38/galGal6 data only. ( Input : ", paste(organism, collapse = ','), ")\n",
               '---------- Annotation process is halted. ----------\nFinish time : ', date(), '\n'))
  } else {}

  cat('---------- Loading a gene table ----------\n')
  #### 01. Load a gene table
  tab_loc = paste0(dbPath, '/cpgIslandExtUnmasked.txt.gz')
  #tab_loc = system.file('extdata', 'cpgIslandExtUnmasked.txt.gz', package = 'IRFinder')
  tab_gz = gzfile(description = tab_loc, open = "r")
  suppressWarnings(dataTable <- read.delim(file = tab_gz, header = FALSE, stringsAsFactors = FALSE))
  colnames(dataTable) = c('bin', 'chrom', 'start', 'end', 'name',
                          'length', 'cpgNum', 'gcNum', 'perCpg',
                          'perGc', 'obsExp')
  close(tab_gz)
  dataTable = subset(dataTable, dataTable$length >= cpglen)
  dataTable = subset(dataTable, str_detect(dataTable$chrom, '_') == FALSE)
  dataTable[,3] = dataTable[,3]+1

  cat('Done.\n')
  cat('---------- Creating a GRanges object  ----------\n')

  #### 02. Make GR object by a gene table
  gr_cpgs = regioneR::toGRanges(dataTable[,c(2,3,4,7,9)])

  #### 03. Make interval GR objects of a gene table
  ranges = seq(from = range[1], to = range[2], interval)
  gr_cpg_dist = vector("list", length = (length(ranges)-1))
  for(x in 1:(length(ranges)-1)){
    tmp_start = gr_cpgs@ranges@start + ranges[x]
    tmp_end = gr_cpgs@ranges@start + ranges[x+1]-1
    gr_cpg_dist[[x]] = GRanges(seqnames = as.character(gr_cpgs@seqnames),
                                 ranges = IRanges(start = tmp_start,
                                                  end = tmp_end),
                                 strand = '*')
  }

  if(!is.null(randomSet)){
    #### 04. Make random GR object
    tmp = read.delim(file = randomSet, header = TRUE, stringsAsFactors = FALSE)
    gr_random = regioneR::toGRanges(tmp[,c(2,3,3)])
  } else {}
  
  cat('Done.\n')
  cat('---------- Annotating integrated regions  ----------\n')

  #### 05. Make gene information
  inside_cpg = as.data.frame(findOverlaps(hits, gr_cpgs, type = 'any',
                                          ignore.strand = TRUE),
                             stringsAsFactors = FALSE)
  a = as.data.frame(hits[inside_cpg$queryHits,], stringsAsFactors = FALSE)
  b = dataTable[inside_cpg$subjectHits,]
  inside_tab = cbind(a,b)
  
  if(!is.null(randomSet)){
    inside_cpg_ran = as.data.frame(findOverlaps(gr_random, gr_cpgs, type = 'any',
                                                ignore.strand = TRUE),
                                   stringsAsFactors = FALSE)
    a = as.data.frame(gr_random[inside_cpg_ran$queryHits,], stringsAsFactors = FALSE)
    b = dataTable[inside_cpg_ran$subjectHits,]
    inside_ran_tab = cbind(a,b)
    
    dist_cpg_ran = vector("list", length = length(ranges)-1)
    dist_cpg_ran_tab = data.frame(stringsAsFactors = FALSE)
  } else {}

  dist_cpg = vector("list", length = length(ranges)-1)
  dist_cpg_tab = data.frame(stringsAsFactors = FALSE)
  
  for(x in 1:(length(ranges)-1)){
    dist_cpg[[x]] = as.data.frame(findOverlaps(query = hits,
                                                subject = gr_cpg_dist[[x]],
                                                type = "any", ignore.strand = TRUE),
                                   stringsAsFactors = FALSE)
    a = as.data.frame(hits[dist_cpg[[x]]$queryHits,], stringsAsFactors = FALSE)
    b = dataTable[dist_cpg[[x]]$subjectHits,]
    tmp1 = cbind(a,b)

    dist_cpg_tab = rbind(dist_cpg_tab, tmp1)
    
    if(!is.null(randomSet)){
      ## random
      dist_cpg_ran[[x]] = as.data.frame(findOverlaps(query = gr_random,
                                                     subject = gr_cpg_dist[[x]],
                                                     type = "any", ignore.strand = TRUE),
                                        stringsAsFactors = FALSE)
      a = as.data.frame(gr_random[dist_cpg_ran[[x]]$queryHits,], stringsAsFactors = FALSE)
      b = dataTable[dist_cpg_ran[[x]]$subjectHits,]
      tmp2 = cbind(a,b)
      
      dist_cpg_ran_tab = rbind(dist_cpg_ran_tab, tmp2)
    } else {}
  }
  
  #### 06. Count the number of hit in each range
  count_hit_cpgs = vector("list", length = length(ranges)-1)
  
  if(!is.null(randomSet)){
    count_hit_cpgs_ran = vector("list", length = length(ranges)-1)
  } else {}

  for(x in 1:(length(ranges)-1)){
    count_hit_cpgs[[x]] = countOverlaps(query = hits, subject = gr_cpg_dist[[x]], type = "any", ignore.strand = TRUE)
    if(!is.null(randomSet)){
      count_hit_cpgs_ran[[x]] = countOverlaps(query = gr_random, subject = gr_cpg_dist[[x]], type = "any", ignore.strand = TRUE)
    } else {}
  }

  hits_cpgs = as.numeric(apply(as.data.frame(count_hit_cpgs, stringsAsFactors = FALSE), MARGIN = 2, sum))
  if(!is.null(randomSet)){
    hits_cpgs_ran = as.numeric(apply(as.data.frame(count_hit_cpgs_ran, stringsAsFactors = FALSE), MARGIN = 2, sum))
  } else {}

  #### 07. Drawing histograms
  ymax_c = max(hits_cpgs)
  mid = vector("numeric", length = (length(ranges)-1))
  for(x in 1:(length(ranges)-1)){
    mid[x] = (ranges[x] + ranges[x+1])/2
  }
  
  freq_sites_cpgs = rep(x = mid, times = hits_cpgs)
  if(!is.null(randomSet)){
    freq_sites_cpgs_ran = rep(x = mid, times = hits_cpgs_ran)
  } else {}

  cat('Done.\n')
  cat('---------- Drawing histograms ----------\n')

  ## CpG
  png(filename = paste0(outpath, '/Distribution_CpGs_', organism, '.png'),
      width = 1200, height = 700)
  hist_frequency = hist(freq_sites_cpgs/1000, ylim = c(0, ceiling(ymax_c * 1.1)), breaks = ranges/1000,
                        ylab = '#Integration events', xlab = "kbs",
                        probability = FALSE, main = NULL, col = 'cornflowerblue',
                        axes = FALSE, cex.lab = 1.5)
  axis(side = 1, at = ranges/1000, cex.axis = 1.5)
  axis(side = 2, at = seq(0, ceiling(ymax_c * 1.1), 1), cex.axis = 1.5)
  text(0, ceiling(ymax_c * 1.1), labels = 'CpG island', cex = 2, font = 2)
  arrows(0,0,0,ceiling(ymax_c * 1.1)*0.95, length = 0.15, code = 1, lwd = 3)
  arrows(min(ranges/1000), ceiling(ymax_c * 1.1)*0.95, range[1]/1000*0.025, ceiling(ymax_c * 1.1)*0.95, length = 0.1, code = 1)
  arrows(max(ranges/1000), ceiling(ymax_c * 1.1)*0.95, range[2]/1000*0.025, ceiling(ymax_c * 1.1)*0.95, length = 0.1, code = 1)
  text(range[2]/2000, ceiling(ymax_c * 1.1)*0.9, 'Upstream', cex = 1.5, col = 'black')
  text(-(range[2]/2000), ceiling(ymax_c * 1.1)*0.9, 'Downstream', cex = 1.5, col = 'black')
  dev.off()
  
  if(!is.null(randomSet)){
    count_sites_cpgs = plyr::count(freq_sites_cpgs)[,2]
    count_sites_cpgs_ran = plyr::count(freq_sites_cpgs_ran)[,2]
    
    count_data = rbind(count_sites_cpgs, count_sites_cpgs_ran)
    png(paste0(outpath, '/Random_distribution_CpGs_', organism, '.png'), width = 1200, height = 750)
    barplot(count_data, beside = TRUE, ylim = c(0, (max(count_data)+2)),
            main = "Random distribution (CpGs)", xlab = 'Intervals (Kbs)', ylab = "#Integration events",
            col = c('Dodgerblue', 'skyblue'), names.arg = ranges[c(1:((length(ranges)-1)/2), ((length(ranges)+3)/2):length(ranges))])
    dev.off()
  } else {}

  if(!is.null(randomSet)){
    result_list = vector("list", length = 4)
    result_list[[1]] = inside_tab
    result_list[[2]] = dist_cpg_tab
    result_list[[3]] = dist_cpg_ran_tab
    result_list[[4]] = gr_cpgs
    names(result_list) = c('CpG_inside_hits', 'CpG_exp_distribution', 'CpG_random_distribution', 'CpG_data')
  } else {
    result_list = vector("list", length = 3)
    result_list[[1]] = inside_tab
    result_list[[2]] = dist_cpg_tab
    result_list[[4]] = gr_cpgs
    names(result_list) = c('CpG_inside_hits', 'CpG_exp_distribution', 'CpG_data')
  }

  cat('---------- Annotation process is finished. ----------\n')
  cat(paste0('Finish time : ', date(), '\n'))

  return(result_list)
}