Rsamtools workflow

Table of Contents Resources Rsamtools workflow 0. Setup 1. sam to bam (also sorts and makes an .bai index file) 2. Examining the bam file: w hat columns are we interested in? 3. Retrieving header information in a BAM file 4. Which features to extract? -requires an indexed BAM

Resources

• Workshop code
• Rsamtools overview
• Rsamtools Overview Rcode

Rsamtools workflow

The help manual pages can be invoked for any core/contributed function run within this workflow e.g.

 library(Rsamtools)
 ## Display the help page for the scanBam() function
 help(scanBam)

The Rsamtools package provides an interface in R to BAM files produced by samtools and other software, and represent a flexible format for storing ‘short’ reads aligned to reference genomes.

This workflow covers importing BAM files into R. The two paradigms are to specify which references (and potentially their sequence ranges) are of specific interest to import, and what columns of the bam file are of interest. Once imported, specific attributes can be filtered further by subseting specific records of interest.

0. Setup

The first command clears out any existing R objects from your current workspace

 # Clear working directory objects
 rm(list=ls())

Any contributed libraries need to be installed once, but loaded each time they are used with either require(), or library();

 require(ShortRead)
 require(Rsamtools)
 require(chipseq)

Loading help documentation vignette;

 # Getting vignette help
 if(interactive()) vignette("Rsamtools-Overview")

The next commands source some global variables in the file ~/VISG-course-2012/supplementary_Rsamtools_usage/globals.R

# Globals
cat("[ Code to source() ]\n")
cat(readLines("globals.R"), sep="\n")

## source globals code
source("globals.R")

cat("[ path to bam:",fl, "]\n")
print(fl)

1. sam to bam (also sorts and makes an .bai index file)

The Rsamtools function asBam() is an interface to execute command line utilities samtools sort and samtools index and generates a sorted, indexed bam file;

bamDest <- "aln" bamName <- asBam(fl, bamDest, overwrite=TRUE)

2. Examining the bam file: w hat columns are we interested in?

We can specify any valid columns of interest from within the samtools specifications, the function scanBamWhat(), lists all available column names;

 what <- scanBamWhat()
 print(what)

3. Retrieving header information in a BAM file

The function scanBamHeader() retrieves the names and widths of all references in the bam file header information;

 ft    <- scanBamHeader(bamName)[[1]][["targets"]]
 print(ft)

4. Which features to extract? -requires an indexed BAM

Specifying references of interest and their genomics ranges with GRanges();

 which <- GRanges(names(ft), IRanges(1, ft))
 print(which)