Get GC content from BAM files in R

The requested R packages are:

library( Rsamtools )
library( TxDb.Hsapiens.UCSC.hg19.knownGene )
library( ShortRead )

We start from a set of aligned BAM files. The BAM files will be loaded using:

bamfls <- BamFileList( list.files( "...", pattern = "bam$", full.names = TRUE ) )

The reference comes as:

ex   <- unlist( exonsBy( TxDb.Hsapiens.UCSC.hg19.knownGene, "gene" ) )

And the function to extract the GC content is defined as:

getGCcontent <- function( bam_path, ex ) {
    if( missing( ex ) ) {
        stop( "Needed reference as 'ex'. ")
    }
    aln   <- readGappedReads( bam_path )
    olap  <- findOverlaps( aln, ex, ignore.strand = TRUE )
    hits  <- tabulate( queryHits( olap ), length( aln ) )
    keep  <- queryHits( olap ) %in% which( hits == 1 )
    reads <- qseq( aln )[ queryHits( olap )[ keep ] ]
    readGc <- rowSums( alphabetFrequency( reads, as.prob = TRUE )[ , c( "C", "G" ) ] )

    return( mean( readGc ) )
}
Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: