wiki:SOPs/miningSAMBAM

Version 5 (modified by gbell, 12 years ago) ( diff )

--

Mining, summarizing, and processing 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 view -h accepted_hits.bam | awk -F "\t" '$3 ~ /chrM/ || $1 ~ /^@/ {print $0}' | 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, without a BAM index)
From SAM
 awk -F"\t" '$3 == "chr2" {print $1}' mapped_reads.sam | sort -u | wc -l
From BAM
 samtools view mapped_reads.bam | awk -F"\t" '$3 == "chr2" {print $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
Note: See TracWiki for help on using the wiki.