Subsetting/ extracting reads from bam files

1. Regions (make an index first, then select your regions)
samtools index bam_file 
samtools view -b file.bam chr1 > file_chr1.bam 
samtools view -b file.bam chr1:1-10000 > file_chr1_1_10000.bam 

2. Regions given in bed file
samtools view -bL regions.bed file.bam > file_regions.bam 

3. First n reads
samtools view -h file.bam | head -n 10000 | samtools view -bS - > file_n10000.bam

Adding the optional IH tag to SAM files

One of the major complaints about the 2 most often used aligners BWA and Bowtie is its failure to report the NH or IH tag.

The IH tag is an indicator of the number of stored alignments in the SAM file that contains the current query (i.e. the read). This is meaningful for multi-mapped reads if you want to know to how many locations the same read has been mapped (eg. assuming your Bowtie parameter "k" has been set to more than 1).

I've written an awk oneliner that will add this tag to your SAM file. What it does is to iterate the file twice, first to tabulate counts, and second to write the extra tag.

Merging RNA-seq replicates from Tophat

Sometimes, you might want to present your replicates as one consensus sequencing track on viewers like IGV or IGB.

To merge bam files : Use samtools merge.
To merge junctions.bed files : Use the tbedjuncs and combineJuncs scripts at seqscripts

cat junctions_file1.bed junction_file_2.bed junction_file_3.bed | tbed2juncs | combineJuncs > combined.junctions.bed

Ref: http://seqanswers.com/forums/showthread.php?t=3381

SAM file directionality

Understanding direction flags and how to extract them

• http://biobeat.wordpress.com/2013/03/14/separating-bam-files-to-forward-and-reverse-strand-files/
• http://biobeat.wordpress.com/2013/04/29/directional-rna-seq-part-1-extract-strand-information-from-sam-file/

 Paired-end specific discussion:

• http://ppotato.wordpress.com/2010/08/25/samtool-bitwise-flag-paired-reads/
• http://www.novocraft.com/wiki/tiki-index.php?page=Extracting+unmapped+reads+from+a+BAM+file+produced+by+novoalign

Useful samflag utility
http://picard.sourceforge.net/explain-flags.html

Trick : seqMINER

seqMINER is a Chip-seq density clustering and plotting software. A trick to keep clustering from previous computations and plot the same ordering on a new set of bam/bed files: Load your bam file, right click and choose view heatmap (bypassing the "cluster" button). Voila!