== SAM/BAM summarizing and processing ==

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

=== Get the official SAM/BAM file format description  ===

All fields in a SAM/BAM file are explained in the [http://samtools.github.io/hts-specs/SAMv1.pdf Sequence Alignment/Map Format Specification].

=== 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 -T foo > foo.sorted.bam
}}}
{{{
# 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
# Or on a folder of SAM files
for samFile in `/bin/ls *.sam`; do bsub /nfs/BaRC_Public/BaRC_code/Perl/SAM_to_BAM_sort_index/SAM_to_BAM_sort_index.pl $samFile ; done
}}}

=== Differences between SAM and BAM files ===

  * A BAM file is a binary version of a SAM file.
  * Both contain identical information about reads and their mapping.
  * A BAM file requires a header but a SAM file may not have one.  (Use 'samtools view -h reads.bam' to print the header with the mapped reads.)
  * Many operations (such as sorting and indexing) work only on BAM files.
  * For almost any application that requires SAM input, this can be created on the fly from a BAM file (using 'samtools view reads.bam |').
  * BAM files take up much less space than SAM files.
  * For archiving purposes, keep only the BAM file.  The SAM file can easily be regenerated (if ever needed).

=== Modify 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 > mapped_unmapped.flagstat.txt
}}}

=== 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
}}}

=== Keep only properly paired reads ===


{{{
samtools view -hS -f 2 mapped_unmapped.sam > mapped.properly_paired.sam
}}}


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

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

=== Get only uniquely mapped reads from sam/bam files generated by bowtie2 ===

{{{
# For sam files:
grep -v "XS:i" | grep "AS:i" foo.sam >| foo_uniquely_mapped.sam
# For bam files:
samtools view -h  foo.bam |grep -E -v "XS:i"|grep -E "@|AS:i" |samtools view -b - >| foo_uniquely_mapped.bam
}}}


=== 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
}}}


=== Get a list of multi-mapped reads, including the number of times each one was mapped ===

Tophat/bowtie mappers create the tag NH:i:XXX where XXX is the number of times the read has mapped.
{{{
bsub "samtools view accepted_hits.bam | grep -v NH:i:1 | perl -pe 's/AS.+(NH:i:\d+)/\$1/' | cut -f1,10,12 | perl -pe 's/NH:i://' | sort -u -k3,3nr > Multi-mapped.sorted.txt"
# Output format:
# read_ID<tab>read<tab>number times mapped
}}}

=== SAM flag explanation ===

   * An explanation of all flags is found in the [http://samtools.github.io/hts-specs/SAMv1.pdf Sequence Alignment/Map Format Specification].
   * Converting a flag into its components may be easiest with the Picard [https://broadinstitute.github.io/picard/explain-flags.html Decoding SAM flags] tool.



=== Split by strand by matched strand ===

{{{
# input: 	accepted_hits.bam
# output:	accepted_hits_negStrand.bam: mapped to negative strand
#		accepted_hits_posStrand.bam: mapped to positive strand

bsub "samtools view -f 16 -b accepted_hits.bam >| accepted_hits_negStrand.bam"
bsub "samtools view -F 16 -b accepted_hits.bam >| accepted_hits_posStrand.bam"

}}}

=== Split reads by pair ===
{{{
# input:	accepted_hits.bam
# output:	1st pair: accepted_hits_1stPair.bam
#		2nd pair: accepted_hits_2ndPair.bam
bsub "samtools view -b -f 0x0040 accepted_hits.bam > accepted_hits_1stPair.bam"
bsub "samtools view -b -F 0x0040 accepted_hits.bam > accepted_hits_2ndPair.bam"
}}}