Create Multiple Y Axes for an R Plot

R base graphics is a powerful tool for plotting data.  Sometimes it is convenient to visualize combinations of variables on the same plot.  Often different variables require different scales.  This can be facilitated by using different Y axes, such as a plot with a Y 2 axis on the right hand side.  In R adding a Y axis is very easy to do.  Here are two simple example scripts you can use to build off of.  The first two examples are two different ways of showing different scales for the same variable (i.e. temperature in Farenheit and Celsus).  The third example is two different variables overlaid on the same R plot with two Y axes used to show the scales of each variable.

Sort BAM File in Samtools

Samtools does a host of useful operations for .bam files.  One such operation is sorting.  Below is a simple example script to show how to use samtools to sort an unsorted BAM file.


This script will sort the unsorted.bam using 8 threads that allocates 12G of memory per thread.  The resulting sorted .bam file will be called sorted.bam.

To confirm a .bam file is sorted, check the header (samtools view -H sorted.bam) for the line:
SO:coordinate.

Additionally, here is the usage information for Samtools sort:

Usage:   samtools sort [options]

Options: -n        sort by read name
         -f        use as full file name instead of prefix
         -o        final output to stdout
         -l INT    compression level, from 0 to 9 [-1]
         -@ INT    number of sorting and compression threads [1]
         -m INT    max memory per thread; suffix K/M/G recognized [768M]

Zero Fill Coverage Gaps in Samtools Depth Output

Merging depth output from multiple .bam files can be difficult since Samtools only outputs depth counts for coordinates with non-zero coverage.  If you want to merge depth output from these .bam files you first need to fill in the base pair positions of no coverage with zero values so the depth output for all .bam files is the same length.  Then using a simple UNIX cat you can merge multiple .bam file depth output into one file for comparison and analysis.

Here is a simple Python script to zero fill Samtools depth output:



And below is an example of how to use it:

Samtools Coverage Depth for Multiple Regions of Interest

Here's some code to use Samtools to extract the sequencing coverage depth from a BAM file for multiple regions of interest as specified in a BED file.  The resulting filename_depth.txt gives coverage data for each base in the region that is in the BAM file.