RNA_seq_coverage.R

plotRNAcoverage <- function(rs, tx, what=c("transcripts", "exons", "introns"), 
                            nsamples=100, main="RNA-seq bias plot", verbose="TRUE") {
    what <- match.arg(what)
    if (what=="transcripts") {
        genes <- exonsBy(txdb, "tx")
        values(genes) <- NULL
    } else if (what=="exons") {
        genes <- exons(txdb)
        values(genes) <- NULL
        genes <- split(genes, 1:length(genes))
    } else if (what=="introns") {
        genes <- unlist(intronsByTranscript(txdb))
        names(genes) <- NULL
        genes <- genes[!duplicated(as.data.frame(genes)),]
        genes <- split(genes, 1:length(genes))
    }
    if (verbose) message("Calculating sampling positions")    
    genesSamp <- rep(genes@unlistData[genes@partitioning@end], each=nsamples)
    values(genesSamp) <- NULL

    fixEnd <- function(x) {
        c(x[1], diff(x))
    }
    if (verbose) message("Determing gene lengths")    
    geneLengths <- tapply(ranges(genes)@unlistData@width, rep(1:length(genes), 
                        fixEnd(ranges(genes)@partitioning@end)), sum)
    sampBP <- unlist(sapply(geneLengths, function(u) {
        as.integer(seq(1, u, length.out=nsamples))
    }, simplify=FALSE, USE.NAMES=FALSE)) + rep(cumsum(c(0, geneLengths[-length(geneLengths)])), each=nsamples)
    ranges(genesSamp) <- IRanges(as.integer(unlist(ranges(genes)))[sampBP], width=1)
    plusSamp <- as(genesSamp[strand(genesSamp)=="+"], "RangesList")
    minusSamp <- as(genesSamp[strand(genesSamp)=="-"], "RangesList")

    if (verbose) message("Creating coverage objects")    
    plusCov <- coverage(grglist(rs[strand(rs)=="+"]), width=as.list(seqlengths(genes)))
    minusCov <- coverage(grglist(rs[strand(rs)=="-"]), width=as.list(seqlengths(genes)))

    if (verbose) message("Sampling sense coverage")    
    genesCov <- rbind(do.call(rbind, lapply(names(plusCov), function(chr) 
        matrix(plusCov[[chr]][plusSamp[[chr]], drop=TRUE], ncol=nsamples, byrow=TRUE))),
                      do.call(rbind, lapply(names(minusCov), function(chr) 
        matrix(minusCov[[chr]][minusSamp[[chr]], drop=TRUE], ncol=nsamples, byrow=TRUE))))/length(rs)*1e6

    if (verbose) message("Sampling antisense coverage")    
    genesACov <- rbind(do.call(rbind, lapply(names(plusCov), function(chr) 
        matrix(plusCov[[chr]][minusSamp[[chr]], drop=TRUE], ncol=nsamples, byrow=TRUE))),
                      do.call(rbind, lapply(names(minusCov), function(chr) 
        matrix(minusCov[[chr]][plusSamp[[chr]], drop=TRUE], ncol=nsamples, byrow=TRUE))))/length(rs)*1e6

    if (verbose) message("Plotting coverage")    
    matplot(x=1:nsamples/nsamples*100, y=cbind("Sense"=colMeans(genesCov), "Antisense"=colMeans(genesACov)), 
    type="b", col=c("black", "red"), pch=19, xlab=paste("% along", what, "length"), ylab="Normalised Coverage", main=main)
    legend("topright", fill=c("black", "red"), legend=c("Sense", "Antisense"))
    invisible(list(genesCov=genesCov, genesACov=genesACov))
}