R/exportUCSC.R

In daewoooo/HapSCElocatoR: Find SCEs in haploid Strand-seq data

#' Generates a UCSC fromated bed file
#'
#' Write a bedfile from yeastSCElocatoR files for upload into a UCSC Genome browser
#'
#' @param index is used to index the bedfile(s)
#' @param outputDirectory is location to write bedfile(s)
#' @param fragments A \code{\link{GRanges}} object with strand and mapq metadata, such as that generated by \code{\link{bam2GRanges}}
#' @param segments A data.frame object with localized regions of 0 and 1 or mixture 0,1 
#' @param recombs A data.frame object with coordinates of every state transition
#' @author Ashley Sanders, David Porubsky
#' @export

exportUCSC <- function(index, outputDirectory, fragments=NULL, segments=NULL, recombs=NULL, col="103,139,139 243,165,97") {
  
  ## Insert chromosome for in case it's missing
  insertchr <- function(gr) {
    mask <- which(!grepl('chr', seqnames(gr)))
    mcols(gr)$chromosome <- as.character(seqnames(gr))
    mcols(gr)$chromosome[mask] <- sub(pattern='^', replacement='chr', mcols(gr)$chromosome[mask])
    mcols(gr)$chromosome <- as.factor(mcols(gr)$chromosome)
    return(gr)
  }
  
  ## Write read fragments to file
  if (!is.null(fragments)) {
    fragments <- insertchr(fragments)
    savefile.reads <- file.path(outputDirectory, paste0(index, '_reads.bed.gz'))
    savefile.reads.gz <- gzfile(savefile.reads, 'w')
    head_reads <- paste0('track name=', index, '_reads visibility=1 colorByStrand="', col, '"') 
    write.table(head_reads, file=savefile.reads.gz, row.names=FALSE, col.names=F, quote=FALSE, append=F, sep='\t')   
    if (length(fragments)>0) {
      bedfile <- as.data.frame(fragments)[c('chromosome','start','end', 'mapq', 'strand')]
      bedfile$name <- 0 # adds a column of 0 as the'name' in the bedfile
      bedfile <- bedfile[c('chromosome','start','end','name','mapq', 'strand')]
    } else {
      bedfile <- data.frame(chromosome='chrI', start=1, end=1, name=NA, mapq=0, strand='.')
    }
    write.table(bedfile,file=savefile.reads.gz, row.names=FALSE, col.names=F, quote=FALSE, append=T, sep='\t')
    close(savefile.reads.gz)
  }
  
  if (!is.null(segments)) {
    mask <- which(!grepl('chr', segments$chromosome))
    segments$chromosome <- sub(pattern='^', replacement='chr', segments$chromosome[mask])
    
    savefile.segments <- file.path(outputDirectory, paste0(index, '_segments.bed.gz'))
    savefile.segments.gz <- gzfile(savefile.segments, 'w')
    head_reads <- paste0('track name=', index, '_segments visibility=1 itemRgb=On') 
    write.table(head_reads, file=savefile.segments.gz, row.names=FALSE, col.names=F, quote=FALSE, append=F, sep='\t')   
    if (nrow(segments)>0) {
      bedfile <- as.data.frame(segments)[c('chromosome','start','end')]
      bedfile$name <- 0 # adds a column of 0 as the'name' in the bedfile
      bedfile$score <- 0
      bedfile$strand <- '.'
      bedfile$thickStart <- 0
      bedfile$thickEnd <- 0
      segments$direction <- factor(segments$direction)
      levels(segments$direction) <- c("103,139,139", "243,165,97", "0,102,0")
      bedfile$itemRgb <- segments$direction
      bedfile <- bedfile[c('chromosome','start','end','name','score','strand','thickStart','thickEnd','itemRgb')]
    } else {
      bedfile <- data.frame(chromosome='chrI', start=1, end=1, name=NA, mapq=0, strand='.')
    }
    write.table(bedfile,file=savefile.segments.gz, row.names=FALSE, col.names=F, quote=FALSE, append=T, sep='\t')
    close(savefile.segments.gz)
  }
  
  if (!is.null(recombs)) {
    mask <- which(!grepl('chr', recombs$chromosome))
    recombs$chromosome <- sub(pattern='^', replacement='chr', recombs$chromosome[mask])
    
    savefile.recombs <- file.path(outputDirectory, paste0(index, '_recombs.bed.gz'))
    savefile.recombs.gz <- gzfile(savefile.recombs, 'w')
    head_br <- paste('track name=', index, '_recombs visibility=dense color=75,125,180', sep="")
    write.table(head_br,  file=savefile.recombs.gz, row.names=FALSE, col.names=F, quote=FALSE, append=F, sep='\t')
    if (nrow(recombs)>0) {
      bpG <-  as.data.frame(recombs)[c('chromosome','start','end')]
    } else {
      bpG <- data.frame(chromosome='chrI', start=1, end=1)
    }
    write.table(bpG, file=savefile.recombs.gz, row.names=FALSE, col.names=F, quote=FALSE, append=T, sep='\t')
    close(savefile.recombs.gz)
  }
}