BinCounts: Get mean read coverage from 5’ to 3’ across the most highly...

In njmadrid/RNAseqQuality: Provides a set of functions for summarizing the quality of RNA-seq data.

Description

Note that the genes/CDS in TopCDSDF are only considered to be suggestions of genes with high expression. For a variety of reasons (e.g. expression of certain CDS may be very rare in some genes), there may not be any expression in any of the bam files for a CDS within a gene. Those CDS will be ignored.

Usage

BinCounts(TopCDSDF, PileUpPerFile, NumBins = 100)

Arguments

TopCDSDF	Data frame of CDS returned from select(Homo.sapiens).
PileUpPerFile	List of Rsamtools pileups per file returned from GetPileUpPerGene().
NumBins	Number of bins to use when parsing the counts for each gene. Default value is 100.

Value

A list of the median, mean, and 3rd quartile bin counts.

Author(s)

Nathaniel J. Madrid, Jason Byars

Examples

BamFilePaths <- dir(path="/home/ubuntu/BamFiles", pattern=".*.bam$", full.names=T)

TopGenesDF <- data.frame(gene = c("NM_016824", "NM_019903", "NP_001112"))
cdsColumns <- c("CDSCHROM", "CDSSTART", "CDSEND", "CDSSTRAND", "SYMBOL")
TopCDSDF <- select(Homo.sapiens, keys=as.character(TopGenesDF$gene), columns=cdsColumns, keytype="REFSEQ")
TopCDSDF$CDSCHROM <- sub("chr", "", TopCDSDF$CDSCHROM)
TopCDSDF <- TopCDSDF[which(TopCDSDF$CDSCHROM %in% c(1:22, "X", "Y")), ]

# BIN THE COUNTS

PileUpPerFile <- lapply(BamFilePaths[1:3], function(f) GetPileUpPerGene(f, TopCDSDF)) 

BinnedCounts <- BinCounts(TopCDSDF, PileUpPerFile, NumBins=100)
BinnedCountsMedian <- BinnedCounts[[1]]
BinnedCountsMean <- BinnedCounts[[2]]
BinnedCounts3rdQ <- BinnedCounts[[3]]

# CREATE PLOTS

g <- ggplot(BinnedCountsMedian, aes(x=BinID, y=RescaledCount, color=FileID)) + geom_line()
g <- g + ylab(label="Normalized Read Coverage") + xlab("Normalized Distance Along Transcript 5' -> 3' (%)")
g <- g + ggtitle(paste(format(length(unique(TopGenesDF$gene)), big.mark=","), " Genes", sep=""))
g

njmadrid/RNAseqQuality documentation built on May 20, 2019, 3:32 p.m.