Summarizing, processing and quality control(QC) of SAM/BAM files

Many of these involve ​| 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

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

QC to get a (visual) summary of mapping statistics. For eg. coverage/distribution of mapped reads across the genome or transcriptome.

RSeQC: RNA-Seq quality control package for getting mapping statistics (eg. unique/multi-mapped reads)

bam_stat.py -i myFile.bam

Picard: CollectRnaSeqMetrics.jar to find coverage across gene body for 5' or 3' bias

java -jar /usr/local/share/picard-tools/CollectRnaSeqMetrics.jar INPUT=accepted_hits.bam REF_FLAT=refFlat.txt STRAND_SPECIFICITY=NONE OUTPUT=Out_RnaSeqMetrics.txt REFERENCE_SEQUENCE=hg19.fa CHART_OUTPUT=Out_RnaSeqMetrics.pdf

QualiMap: can be used on DNA or RNA-Seq to get summary of mapping and 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"
#rnaseq qc
bsub "qualimap rnaseq -bam myFile.bam -gtf Homo_sapiens.GRCh37.72.canonical.gtf -outdir output_qualimap_rnaseq -protocol non-strand-specific"