R/exportUCSC.R

In daewoooo/BreakPointR: Find breakpoints in Strand-seq data

#' Export UCSC browser formated files
#'
#' Write a bedfile or bedgraph from a breakpointR object for upload on to the UCSC Genome browser.
#'
#' @param index A character used to name the bedfile(s).
#' @param outputDirectory Location to write bedfile(s).
#' @param fragments A \code{\link{GRanges-class}} object with strand and mapq metadata, such as that generated by \code{\link{readBamFileAsGRanges}}
#' @param deltaWs A \code{\link{GRanges-class}} object with metadata column "deltaW" generated by \code{\link{deltaWCalculator}}.
#' @param breakTrack A \code{\link{GRanges-class}} object with metadata "genoT" (e.g. newBreaks) will write a bedtrack with refined breakpoints.
#' @param confidenceIntervals A \code{\link{GRanges-class}} object with metadata "genoT" the same length as \code{breakTrack} (e.g. confint) will write a bedtrack with breakpoints confidence intervals.
#' @param breaksGraph A \code{\link{GRanges-class}} object.
#' @return \code{NULL}
#' @author Ashley Sanders, David Porubsky, Aaron Taudt
#' @importFrom utils write.table
#' @export
#' @examples
#'## Get an example file
#'exampleFolder <- system.file("extdata", "example_results", package="breakpointRdata")
#'exampleFile <- list.files(exampleFolder, full.names=TRUE)[1]
#'## Load the file 
#'brkpts <- get(load(exampleFile))
#'## Write results to BED files
#'breakpointr2UCSC(index='testfile', outputDirectory=tempdir(), breakTrack=brkpts$breaks)
         
breakpointr2UCSC <- function(index, outputDirectory, fragments=NULL, deltaWs=NULL, breakTrack=NULL, confidenceIntervals=NULL, breaksGraph=NULL) {
  
    ## Write read fragments to file
    if (!is.null(fragments) & length(fragments) > 0) {
        fragments <- insertchr(fragments)
        savefile.reads <- file.path(outputDirectory, paste0(index, '_reads.bed.gz'))
        ptm <- startTimedMessage("Writing to file ", savefile.reads, " ...")
        savefile.reads.gz <- gzfile(savefile.reads, 'w')
        col="103,139,139 243,165,97" # Set specific Strand-seq colors
        head_reads <- paste0('track name=', index, '_reads visibility=1 colorByStrand="', col, '"') 
        utils::write.table(head_reads, file=savefile.reads.gz, row.names=FALSE, col.names=FALSE, quote=FALSE, append=FALSE, 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='chr1', start=1, end=1, name=NA, mapq=0, strand='.')
        }
        utils::write.table(bedfile,file=savefile.reads.gz, row.names=FALSE, col.names=FALSE, quote=FALSE, append=TRUE, sep='\t')
        close(savefile.reads.gz)
        stopTimedMessage(ptm)
    }

    ## Write breakpoints and confidence intervals to file
    if (!is.null(breakTrack) & length(breakTrack) > 0) {
        breakTrack <- insertchr(breakTrack)
        savefile.breakPoints <- file.path(outputDirectory, paste0(index, '_breakPoints.bed.gz'))
        ptm <- startTimedMessage("Writing to file ", savefile.breakPoints, " ...")
        savefile.breakPoints.gz <- gzfile(savefile.breakPoints, 'w')
        head_br <- paste('track name=', index, '_BPs description=BedGraph_of_breakpoints_',index,'_allChr visibility=dense color=75,125,180', sep="")
        utils::write.table(head_br, file=savefile.breakPoints.gz, row.names=FALSE, col.names=FALSE, quote=FALSE, append=FALSE, sep='\t')
        if (length(breakTrack)>0) {
            bpG <-  as.data.frame(breakTrack)[c('chromosome','start','end','genoT')]
        } else {
            bpG <- data.frame(chromosome='chr1', start=1, end=1, genoT=NA)
        }
        if (!is.null(confidenceIntervals)) {
            bpG$score <- 0
            bpG$strand <- '.'
            bpG$thickStart <- bpG$start
            bpG$thickEnd <- bpG$end
            bpG$start <- start(confidenceIntervals)
            bpG$end <- end(confidenceIntervals)
        }
        utils::write.table(bpG, file=savefile.breakPoints.gz, row.names=FALSE, col.names=FALSE, quote=FALSE, append=TRUE, sep='\t')
        close(savefile.breakPoints.gz)
        stopTimedMessage(ptm)
    }

    ## Write deltaWs to file
    if (!is.null(deltaWs) & length(deltaWs) > 0) {
        deltaWs <- insertchr(deltaWs)
        savefile.deltaWs <- file.path(outputDirectory, paste0(index, '_deltaWs.bed.gz'))
        ptm <- startTimedMessage("Writing to file ", savefile.deltaWs, " ...")
        savefile.deltaWs.gz <- gzfile(savefile.deltaWs, 'w')
        head_dW <- paste('track type=bedGraph name=', index,'_dWs description=BedGraph_of_deltaWs_',index, '_allChr visibility=full color=200,100,10', sep="")
        utils::write.table(head_dW, file=savefile.deltaWs.gz, row.names=FALSE, col.names=FALSE, quote=FALSE, append=FALSE, sep='\t')   
        if (length(deltaWs)>0) {
            bedG <- as.data.frame(deltaWs)[c('chromosome','start','end','deltaW')]
        } else {
            bedG <- data.frame(chromosome='chr1', start=1, end=1, deltaW=NA)
        }
        utils::write.table(bedG, file=savefile.deltaWs.gz, row.names=FALSE, col.names=FALSE, quote=FALSE, append=TRUE, sep='\t')
        close(savefile.deltaWs.gz)
        stopTimedMessage(ptm)
    }

    ## Write breakpoint summary to file
    if (!is.null(breaksGraph) & length(breaksGraph) > 0) {
        breaksGraph <- insertchr(breaksGraph)
        savefile.breaksGraph <- file.path(outputDirectory, paste0(index, '_BreaksGraph.bed.gz'))
        ptm <- startTimedMessage("Writing to file ", savefile.breaksGraph, " ...")
        savefile.breaksGraph.gz <- gzfile(savefile.breaksGraph, 'w')
        head_breaks <- paste('track type=bedGraph name=', index,' description=BedGraph_of_breakpoints_',index, '_allChr visibility=dense color=125,35,190', sep="")
        utils::write.table(head_breaks, file=savefile.breaksGraph.gz, row.names=FALSE, col.names=FALSE, quote=FALSE, append=FALSE, sep='\t')
        if (length(breaksGraph)>0) {
            bedG <- as.data.frame(breaksGraph)[c('chromosome','start','end','hits')]
        } else {
            bedG <- data.frame(chromosome='chr1', start=1, end=1)
        }
        utils::write.table(bedG, file=savefile.breaksGraph.gz, row.names=FALSE, col.names=FALSE, quote=FALSE, append=TRUE, sep='\t')
        close(savefile.breaksGraph.gz)
        stopTimedMessage(ptm)
    }
}


#' Generates a bedfile from an input GRanges file
#'
#' Write a bedfile from Breakpoint.R files for upload on to UCSC Genome browser
#'
#' @param gr A \code{\link{GRanges-class}} object with genomic ranges to be exported into UCSC format.
#' @inheritParams breakpointr2UCSC
#' @param colorRGB An RGB color to be used for submitted ranges.
#' @return \code{NULL}
#' @author David Porubsky
#' @importFrom utils write.table
#' @export
#' @examples
#'## Get an example file
#'exampleFolder <- system.file("extdata", "example_results", package="breakpointRdata")
#'exampleFile <- list.files(exampleFolder, full.names=TRUE)[1]
#'## Load the file 
#'counts <- get(load(exampleFile))[['counts']]
#'## Export 'wc' states into a UCSC formated file
#'ranges2UCSC(gr=counts[counts$states == 'wc'], index='testfile', outputDirectory=tempdir())

ranges2UCSC <- function(gr, outputDirectory=".", index="bedFile", colorRGB='0,0,0') {
  
    ## Insert 'chr' before chromosome number if missing
    gr <- insertchr(gr)
  
    ## Prepare file for bed export
    savefile<- file.path(outputDirectory, paste0(index, '.bed.gz'))
    savefile.gz <- gzfile(savefile, 'w')
    ptm <- startTimedMessage("Writing to file ", savefile, " ...")
  
    ## Write a header to the file
    header <- paste('track name=', index, ' description=Bed_regions_',index,' visibility=dense color=',colorRGB, sep="")
    utils::write.table(header, file=savefile.gz, row.names=FALSE, col.names=FALSE, quote=FALSE, append=FALSE, sep='\t')
  
    ## Write the rest of the file  
    if (length(gr)>0) {
        bedF <- as.data.frame(gr)[c('chromosome','start','end')]
    } else {
        bedF <- data.frame(chromosome='chr1', start=1, end=1)
    }
    utils::write.table(bedF, file=savefile.gz, row.names=FALSE, col.names=FALSE, quote=FALSE, append=TRUE, sep='\t')
    close(savefile.gz)
  
    stopTimedMessage(ptm)
}