-
Notifications
You must be signed in to change notification settings - Fork 0
Rsamtools workflow
In R we can also the documentation using tkWidgets, in order to do this we must install the R package;
Now to interact with vignette documentation and code chunks
library(tkWidgets) vExplorer()
Select the package Rsamtools => Rsamtools-Overview.Rnw, and click on view PDF. Code chunks can be executed in sequential order while reading the documentation.
Go to the Rsamtools workflow directory;
cd $HOME/VISG-course-2012/supplementary_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.
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)
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)
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)
The function scanBamHeader() retrieves the names and widths of all references in the bam file header information;
ft <- scanBamHeader(bamName)[[1]][["targets"]] print(ft)
Specifying references of interest and their genomics ranges with GRanges();
which <- GRanges(names(ft), IRanges(1, ft)) print(which)
Merging the which and what filtering criterion into a set of parameters using ScanBamParam;
param <- ScanBamParam(which=which, what=what)
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)
This code snippet illustrates that memory issues associated with extremely large bam files can be overcome by iterating through the entire file processing a set number of records at a time;
## See help(BamFile), 'yieldSize' argument alters the number of records ## to yield each time the file is read from using 'scanBam' allowing ## memory efficient loading
bfl <- open(BamFile(bamName, yieldSize=1000))
while (nrec <- length((bam <- scanBam(bfl))[[1]][[1]])) {
cat("records:", nrec, "\n")
## Further Processing of object bam...
}
close(bfl)
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)
names(bam)
[1] "CO_Pool1_contig00004:1-1928" "EDH2_Pool1_contig00005:1-853" [3] "FT1_Pool1_contig00001:1-3360"
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)
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")
}
We can select a reference of interest by name, or numeric index;
identical(bam[["CO_Pool1_contig00004:1-1928"]], bam[[1]])
Returning the first few vector elements of the cigar string;
head(bam[[1]][["cigar"]])
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)
Printing all the cigar strings which do not contain INDELS;
print(bam[[1]][["cigar"]][noINDELS])
Calculating the nucleotide frquencies using the function alphabetByCycle, this returns an n by p matrix matrix of alphabets versus cycle;
abc <- alphabetByCycle(bam[[1]][["seq"]])
Investigating the non ambiguous IUPAC nucleotide codes, A, C, G, T;
print(abc[,1:4])
Creating a plot of nucleotide frequency versus cycle for each nucleotide A, C, G, T. The variable PNG=TRUE/FALSE in the conditional statements will create a portable network graph (png) file if PNG=TRUE in ./plots/alphabetPlot.png;
if(PNG) png(file.path(plotDir, "alphabetPlot.png")) matplot(t(abc[1:4, ]), type="l", lty=1, lwd=1, ylab="Nuclotide frequency", main=names(bam)[1]) if(PNG) dev.off()
Same alphabet plot, but this time focusing on the first 50 cycles;
if(PNG) png(file.path(plotDir, "alphabetPlot_50cycles.png")) matplot(t(abc[1:4,1:50]), type="l", lty=1, lwd=1, ylab="Nuclotide frequency", main=names(bam)[1]) if(PNG) dev.off()
Finding the indices of the + (plus=forward) an - (minus = reverse complement) strands;
indPos <- which(bam[[1]][["strand"]] == "+") indNeg <- which(bam[[1]][["strand"]] == "-")
Alphabet by cycle for each strand;
abc.P <- alphabetByCycle(bam[[1]][["seq"]][indPos]) abc.N <- alphabetByCycle(bam[[1]][["seq"]][indNeg])
Repeating the alphabet plot for the subset of reads for + and - strands
This is an illustration extracting bam elements for the first reference to create a ShortReadQ object (See Quality control workflow);
fq <- ShortReadQ(bam[[1]][["seq"]], bam[[1]][["qual"]], BStringSet(bam[[1]][["qname"]]))
We can generate a quality report using qa() directly from the bam file or from the ShortReadQ object generated in (22);
## QC directly from bam File qaSummary <- qa(bamName, type="BAM") ## , lane="Roche 454")
## QC from ShortReadQ object qaSummaryR <- qa(fq, lane="Roche 454") fq <- ShortReadQ(bam[[1]][["seq"]], bam[[1]][["qual"]], BStringSet(bam[[1]][["qname"]]))
Using the perCycle component to generate a Quality by cycle plot from ShortRead;
perCycle <- qaSummary[["perCycle"]] head(perCycle[["baseCall"]]) head(perCycle[["quality"]])
if(PNG) png(file.path(plotDir, "cycleQuality.png")) ShortRead:::.plotCycleQuality(perCycle$quality, main=names(bam)[1]) if(PNG) dev.off()
Creating a coverage plot (function from chipseq package) using IRanges start and end coordinates;
IRanges <- IRanges(start = bam[[1]][["pos"]], width=bam[[1]][["qwidth"]]) Cov <- coverage(IRanges) Peaks <- slice(Cov, 0)
if(PNG) png(file.path(plotDir, "coverage.png")) coverageplot(Peaks, main=names(bam)[1]) if(PNG) dev.off()
Using the indices evaluated in 17, to calculate coverage information for +/- strands;
IRangesF <- IRanges(start = bam[[1]][["pos"]][indPos], width=bam[[1]][["qwidth"]][indPos]) CovF <- coverage(IRangesF) PeaksF <- slice(CovF, 0)
IRangesR <- IRanges(start = bam[[1]][["pos"]][indNeg], width=bam[[1]][["qwidth"]][indNeg]) CovR <- coverage(IRangesR) PeaksR <- slice(CovR, 0)
Repeating coverage plot visualizing + strand coverage depth is positive on the y-axis, - strand coverage is negative on he y-axis;
if(PNG) png(file.path(plotDir, "coverageStrands.png")) coverageplot(PeaksF, PeaksR, main=names(bam)[1]) if(PNG) dev.off()