R/writeBedFile.R

In drashley/InversionAnalysis_HGSVC: Find breakpoints in Strand-seq data

#' Generates a bedfile from inFile
#'
#' Write a bedfile from Breakpoint.R files for upload on to UCSC Genome browser
#'
#' @param fileName is used to index the bedfile(s)
#' @param dataDirectory 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 deltaWs A \code{\link{GRanges}} object with metadata column "deltaW" generated by \code{\link{deltaWCalculator}}
#' @param breakTrack A \code{\link{GRanges}} object with metadata "genoT" (e.g. newBreaks) will write a bedtrack with refined breakpoints
#' @author Ashley Sanders, David Porubsky, Aaron Taudt
#' @export

writeBedFile <- function(fileName, dataDirectory, fragments=NULL, deltaWs=NULL, breakTrack=NULL, bin=NA, 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)
	}

  index <- strsplit(fileName, '\\.')[[1]][1] # pulls out first element in fileName
	## Write read fragments to file
	if (!is.null(fragments)) {
		fragments <- insertchr(fragments)
		savefile.reads <- file.path(dataDirectory, paste0('AllChr_', basename(fileName), '_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')]
			write.table(bedfile,file=savefile.reads.gz, row.names=FALSE, col.names=F, quote=FALSE, append=T, sep='\t')
		}
		close(savefile.reads.gz)
	}

	## Write breakPoints to file
	if (!is.null(breakTrack)) {
		breakTrack <- insertchr(breakTrack)
		savefile.breakPoints <- file.path(dataDirectory, paste0('AllChr_', basename(fileName), '_breakPoints.bed.gz'))
		savefile.breakPoints.gz <- gzfile(savefile.breakPoints, 'w')
		head_br <- paste('track name=', index, '_BPs description=BedGraph_of_breakpoints_',basename(fileName),'_allChr visibility=dense color=75,125,180', sep="")
		write.table(head_br,  file=savefile.breakPoints.gz, row.names=FALSE, col.names=F, quote=FALSE, append=F, sep='\t')
		if (length(breakTrack)>0) {
			bpG<-  as.data.frame(breakTrack)[c('chromosome','start','end','genoT')]
			write.table(bpG, file=savefile.breakPoints.gz, row.names=FALSE, col.names=F, quote=FALSE, append=T, sep='\t')
		}
		close(savefile.breakPoints.gz)
	}

	## Write deltaWs to file
	if (!is.null(deltaWs)) {
		deltaWs <- insertchr(deltaWs)
		savefile.deltaWs <- file.path(dataDirectory, paste0('AllChr_', basename(fileName), '_deltaWs.bed.gz'))
		savefile.deltaWs.gz <- gzfile(savefile.deltaWs, 'w')
		head_dW<- paste('track type=bedGraph name=', index,'_dWs description=BedGraph_of_deltaWs_',basename(fileName), '_allChr visibility=dense color=200,100,10', sep="")
		write.table(head_dW, file=savefile.deltaWs.gz, row.names=FALSE, col.names=F, quote=FALSE, append=F, sep='\t')   
		if (length(deltaWs)>0) {
			bedG <- as.data.frame(deltaWs)[c('chromosome','start','end','deltaW')]
			write.table(bedG, file=savefile.deltaWs.gz, row.names=FALSE, col.names=F, quote=FALSE, append=T, sep='\t')
		}
		close(savefile.deltaWs.gz)
	}

}