== Mining, summarizing, and processing SAM/BAM files ==

Many of these involve [http://samtools.sourceforge.net/samtools.shtml | samtools]

=== Convert, sort, and/or index ===


{{{
Convert SAM to BAM:
 samtools view -bS -o foo.bam foo.sam
}}}

{{{
Convert BAM to SAM:
 samtools view -h -o foo.sam foo.bam
}}}
{{{
Sort BAM file (where ".bam" is added to "foo.sorted")
 samtools sort foo.bam foo.sorted
}}}
{{{
Index a sorted BAM file (which creates foo.sorted.bam.bai):
 samtools index foo.sorted.bam
    
# Both foo.sorted.bam and foo.sorted.bam.bai are needed for visualization.
 }}}

All three steps (SAM=>BAM, sorting, and indexing) can be merged into one command.  See
{{{
/nfs/BaRC_Public/BaRC_code/Perl/SAM_to_BAM_sort_index/SAM_to_BAM_sort_index.pl
}}}

=== Process a BAM file into another BAM file ===

In many cases, there's no need to create an intermediate SAM file.  For example, to extract selected (mapped to chrM) reads:
{{{
samtools index accepted_hits.bam   # Required if you want to select a genome region (like chrM)
samtools view -h accepted_hits.bam chrM | samtools view -bS - > accepted_hits.chrM_only.bam
}}}
We need to keep the header to convert back to BAM (hence the '-h' with 'samtools view' and the '$1 ~ ...' with awk).

=== Count the number of mapped reads ===

{{{
samtools flagstat mapped_unmapped.bam
}}}

=== Count the number of mapped reads by chromosome ===

{{{
Method 1 (all chromosomes)
1 - Index the BAM file: 
 samtools index mapped_reads.bam
2 - Get index statistics (including the number of mapped reads in the third column: 
 samtools idxstats mapped_reads.bam
}}}

{{{
Method 2 (one chromosome at a time, for example, chr2)
From SAM
 awk -F"\t" '$3 == "chr2" {print $1}' mapped_reads.sam | sort -u | wc -l
From BAM
 samtools view mapped_reads.bam chr2 | cut -f 1 | sort -u | wc -l
}}}

=== Remove unmapped reads ===


{{{
samtools view -hS -F 4 mapped_unmapped.sam > mapped_only.sam
}}}


=== How many multiple/uniquely mapped reads are in a bam/sam file? ===

{{{
bam_stat.py -i mapped_reads.bam >& bam_stat.out.txt
}}}

=== View alignment with samtools ===

{{{
-e: change identical bases to '='
samtools view -b accepted_hits.bam |samtools fillmd -e - /nfs/genomes/mouse_mm10_dec_11_no_random/fasta_whole_genome/mm10.fa|more
}}}


=== Use QualiMap to get (visual) summary of mapping statistics eg. coverage/distribution ===
{{{
Graphical interface: enter 'qualimap' on the command line

Command line:
unset DISPLAY  #needed for submitting to cluster
bsub "qualimap bamqc -bam myFile.bam -outdir output_qualimap"
}}}