R/exportSignals.R

Defines functions exportSignals

Documented in exportSignals

#' output signals to file
#'
#' Output signals to bedgraph, bed, wig, etc, for track viewer
#'
#' @param dat An object of \link[GenomicRanges]{GRanges}, 
#' or \link[SummarizedExperiment]{RangedSummarizedExperiment}
#' with assays of raw counts, ratios, background correct ratios,
#' smoothed ratios and z-scores. It should be an element of output of 
#' \link{smoothRatiosByChromosome}
#' @param assayName character(1). Assay name for 
#' \link[SummarizedExperiment]{RangedSummarizedExperiment}
#' @param colName character(1). Column name of metadata of dat or assay of dat 
#' for coverage weight, see
#' \link[GenomicRanges]{coverage}, 
#' \link[SummarizedExperiment]{RangedSummarizedExperiment}.
#' @param con The connection to which data is saved. If this is a character
#' vector, it is assumed to be a filename and a corresponding file connection
#' is created and then closed after exporting the object.
#' If missing, a \link[IRanges]{SimpleRleList} will be returned.
#' @param format The format of the output. see \link[rtracklayer]{export}.
#' @param ... Parameters to pass to \link[rtracklayer]{export}
#' @import GenomicAlignments
#' @import SummarizedExperiment
#' @import GenomicRanges
#' @import rtracklayer
#' @export
#' @return If con is missing, a \link[IRanges]{SimpleRleList} will be returned.
#' Otherwise, nothing is returned.
#' @examples
#' gr <- GRanges("chr1", IRanges(seq_len(100), 201:300), reads=rep(1, 100))
#' myTrackLine <- new("TrackLine", name="my track",
#'                     description="description of my track",
#'                     color=col2rgb("red")[, 1],
#'                     visibility="full")
#' exportSignals(gr, colName="reads", 
#'               con="test.bedGraph", trackLine=myTrackLine)
#' data(triplicates.counts)
#' exportSignals(triplicates.counts, "counts", 
#'               "G26.bam", "test.bw", format="bigWig")
#'
exportSignals <- function(dat, assayName, colName, 
                          con, format="bedGraph", ...){
    stopifnot(inherits(dat, c("GRanges", "RangedSummarizedExperiment")))
    if(class(dat)=="GRanges"){
        gr <- dat
    }else{## RangedSummarizedExperiment
        stopifnot(!missing(assayName))
        stopifnot(assayName %in% names(assays(dat)))
        gr <- rowRanges(dat)
        mcols(gr) <- assays(dat)[[assayName]]
    }
    if(!colName %in% colnames(mcols(gr))){
        stop("colName in not valid.")
    }
    seqlevels(gr) <- seqlevels(gr)[seqlevels(gr) %in% unique(seqnames(gr))]
    dat1 <- coverage(gr, weight = mcols(gr)[, colName])
    dat2 <- coverage(gr)
    dat <- dat1/dat2
    dat[is.na(dat)] <- 0
    if(missing(con)){
        return(dat)
    }else{
        return(export(object=dat, con=con, format=format, ...))
    }
}
Bioconductor-mirror/NADfinder documentation built on July 28, 2017, 4:27 a.m.