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 5. Create a parameter object for scanning BAM files 6. Load sorted BAM file into R 7. bam file: class, length, element classes 8. Each element of the list corresponds to a range specified by 'which' 9. First bam[[1]] list component 10. Each component is a list containing the elements specified by 'what' 11. Referencing by name or index is the same 12. Cigar string see help(cigar-utils) for utility functions 13. grep on cigar string containing INDEL characters I or D 14. Reads not containing INDELS in alignment

Resources

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 select what fields that will be parsed. All valid column names of interest from within the samtools specifications, are listed by the function scanBamWhat();

 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)

5. Create a parameter object for scanning BAM files

Merging the which and what filtering criterion into a set of parameters using ScanBamParam;

 param <- ScanBamParam(which=which, what=what)

6. Load sorted BAM file into R

The 'scanBam' function parses binary BAM files;

 bam <- scanBam(bamName, param = param)
 ## Large files will hit memory limits e.g. 2^32 for 32-bit architecture
 object.size(bam)

7. bam file: class, length, element classes

A bam file loaded into R is essentially a list of lists, the first list level equals the number of references selected by which, the second level contains the what field columns of interest. The following functions query the first level;

 class(bam)
 length(bam)
 sapply(bam, class)

8. Each element of the list corresponds to a range specified by 'which'

names(bam)

 [1] "CO_Pool1_contig00004:1-1928"  "EDH2_Pool1_contig00005:1-853"
 [3] "FT1_Pool1_contig00001:1-3360"

9. First bam[[1]] list component

The following functions query the second level of what field columns for the first reference;

 class(bam[[1]])
 length(bam[[1]])
 sapply(bam[[1]], class)

10. Each component is a list containing the elements specified by 'what'

This code snippet prints a the first few elements of each second level list element;

 for(j in seq_len(length(bam[[1]]))) {
   cat("[", names(bam[[1]])[j], "list element ]\n")
   print(head(bam[[1]][[j]], n=10))
   cat("...\n\n")
 }

11. Referencing by name or index is the same

We can select a reference of interest by name, or numeric index;

 identical(bam[["CO_Pool1_contig00004:1-1928"]], bam[[1]])

12. Cigar string see help(cigar-utils) for utility functions

Returning the first few vector elements of the cigar string;

 head(bam[[1]][["cigar"]])

13. grep on cigar string containing INDEL characters I or D

Using grep, searching for matches containing I (insertion), or D (deletion) within each element of a character vector and returning their indices;

 noINDELS <- grep("[ID]", bam[[1]][["cigar"]], invert=TRUE)

14. Reads not containing INDELS in alignment

Printing all the cigar strings which do not contain INDELS;

 print(bam[[1]][["cigar"]][noINDELS])

Navigation Menu

Rsamtools workflow

Table of Contents

Resources