learning_bam_file.Rmd

---
title: "Learning the BAM format"
output: github_document
---

```{r setup, include=FALSE}
Sys.setenv(PATH=paste0(Sys.getenv("PATH"), ":", getwd()))
knitr::opts_chunk$set(echo = TRUE)
```

## Introduction

![Build README](https://github.com/davetang/learning_bam_file/actions/workflows/create_readme.yml/badge.svg)

SAMtools provides various (sub)tools for manipulating alignments in the SAM/BAM format. The SAM (Sequence Alignment/Map) format (BAM is just the binary form of SAM) is currently the _de facto_ standard for storing large nucleotide sequence alignments. If you are working with high-throughput sequencing data, at some point you will probably have to deal with SAM/BAM files, so familiarise yourself with them! For the latest information on SAMtools, please refer to the [release notes](https://github.com/samtools/samtools/releases).

The examples in this README use the `ERR188273_chrX.bam` BAM file (stored in the `eg` folder) generated as per https://github.com/davetang/rnaseq using the HISAT2 + StringTie2 RNA-seq pipeline. This README is generated using the `create_readme.sh` script; if you want to generate this file yourself, please use [this Docker image](https://hub.docker.com/repository/docker/davetang/r_build) and the `Makefile` in this directory. For example:

```bash
# clone this repo
git clone https://github.com/davetang/learning_bam_file.git
cd learning_bam_file

docker pull davetang/r_build:4.1.2
docker run --rm -it -v $(pwd):/work davetang/r_build:4.1.2 /bin/bash

# inside the Docker container
make
```

## Installing SAMtools

For installing SAMtools, I recommend using `Conda` and the [Bioconda samtools package](https://anaconda.org/bioconda/samtools). I also recommend using [Miniconda](https://docs.conda.io/en/latest/miniconda.html) instead of Anaconda because Anaconda comes with a lot of tools/packages that you will probably not use. I wrote a [short introduction to Conda](https://davetang.github.io/reproducible_bioinformatics/conda.html) if you want to find learn more.

Once you have installed Miniconda, you can install SAMtools as follows:

```bash
conda install -c bioconda samtools
```

Otherwise you can download the source and compile it yourself; change `dir` to the location you want `samtools` to be installed. `samtools` will be installed in `${dir}/bin`, so make sure this is in your `$PATH`.

```bash
#!/usr/bin/env bash

set -euo pipefail

ver=1.15
tool=samtools
url=https://github.com/samtools/${tool}/releases/download/${ver}/${tool}-${ver}.tar.bz2
dir=${HOME}/local

wget ${url}
tar xjf ${tool}-${ver}.tar.bz2
cd ${tool}-${ver}
./configure --prefix=${dir}
make && make install
cd ..

rm -rf ${tool}-${ver} ${tool}-${ver}.tar.bz2

>&2 echo Done
exit 0
```

## Basic usage

If you run `samtools` on the terminal without any parameters or with `--help`, all the available utilities are listed:

```{bash engine.opts='-l'}
samtools --help
```

## Viewing

Use [bioSyntax](https://github.com/bioSyntax/bioSyntax) to prettify your output.

```bash
samtools view aln.bam | sam-less
```

![bioSyntax](img/sam_less.png)

## Converting a SAM file to a BAM file

A BAM file is just a SAM file but stored in binary format; you should always convert your SAM files into BAM format since they are smaller in size and are faster to manipulate.

I don't have a SAM file in the example folder, so let's create one and check out the first ten lines. Note: remember to use `-h` to ensure the SAM file contains the sequence header information. Generally, I recommend storing only sorted BAM files as they use even less disk space and are faster to process.

```{bash engine.opts='-l'}
samtools view -h eg/ERR188273_chrX.bam > eg/ERR188273_chrX.sam
```

Notice that the SAM file is much larger than the BAM file.

Size of SAM file.

```{bash engine.opts='-l'}
ls -lh eg/ERR188273_chrX.sam
```

Size of BAM file.

```{bash engine.opts='-l'}
ls -lh eg/ERR188273_chrX.bam
```

We can use `head` to view a SAM file.

```{bash engine.opts='-l'}
head eg/ERR188273_chrX.sam
```

The lines starting with an "@" symbol contains the header information. The @SQ tag is the reference sequence dictionary; SN refers to the reference sequence name and LN refers to the reference sequence length. If you don't see lines starting with the "@" symbol, the header information is probably missing. You can generate this information again by running the command below, where `ref.fa` is the reference FASTA file used to map the reads.

```bash
samtools view -bT sequence/ref.fa aln.sam > aln.bam
```

If the header information is available, we can convert a SAM file into BAM by using `samtools view -b`. In newer versions of SAMtools, the input format is auto-detected, so we no longer need the `-S` parameter.

```{bash engine.opts='-l'}
samtools view -b eg/ERR188273_chrX.sam > eg/my.bam
```

## Converting a BAM file to a CRAM file

The CRAM format is even more compact. Use `samtools view` with the `-T` and `-C` arguments to convert a BAM file into CRAM.

```{bash engine.opts='-l'}
samtools view -T genome/chrX.fa -C -o eg/ERR188273_chrX.cram eg/ERR188273_chrX.bam

ls -lh eg/ERR188273_chrX.[sbcr]*am
```

You can use `samtools view` to view a CRAM file just as you would for a BAM file.

```{bash engine.opts='-l'}
samtools view eg/ERR188273_chrX.cram | head
```

I have an [old blog post](https://davetang.org/muse/2014/09/26/bam-to-cram/) on the CRAM format.

## Sorting a SAM/BAM file

Many downstream tools require sorted BAM files and since they are slightly more compact than unsorted BAM files, you should always sorted BAM files. In SAMtools version 1.3 or newer, you can directly generate a sorted BAM file from a SAM file.

```{bash engine.opts='-l'}
samtools sort eg/ERR188273_chrX.sam -o eg/sorted.bam
ls -l eg/ERR188273_chrX.bam
ls -l eg/sorted.bam
```

You should use use additional threads (if they are available) to speed up sorting; to use four threads, use `-@ 4`.

Time taken using one thread (default).

```{bash engine.opts='-l'}
time samtools sort eg/ERR188273_chrX.sam -o eg/sorted.bam
```

Time taken using four threads.

```{bash engine.opts='-l'}
time samtools sort -@ 4 eg/ERR188273_chrX.sam -o eg/sorted.bam
```

Many of the SAMtools subtools can use additional threads, so make use of them if you have the resources!

## Creating a BAM index file

Various tools require BAM index files, such as IGV, which is a tool that can be used for visualising BAM files.

```{bash engine.opts='-l'}
samtools index eg/ERR188273_chrX.bam
```

## Adding read groups

Some tools like GATK and Picard require [read groups](https://gatk.broadinstitute.org/hc/en-us/articles/360035890671-Read-groups) (RG). You can add or replace read groups using `samtools addreplacerg`.

```{bash engine.opts='-l'}
samtools addreplacerg -r "@RG\tID:ERR188273\tSM:ERR188273\tPL:illumina" -o eg/ERR188273_chrX_rg.bam eg/ERR188273_chrX.bam
samtools head eg/ERR188273_chrX_rg.bam
```

If you want to replace existing read groups, just use the same command.

```{bash engine.opts='-l'}
samtools addreplacerg -r "@RG\tID:ERR188273_2\tSM:ERR188273_2\tPL:illumina_2" -o eg/ERR188273_chrX_rg2.bam eg/ERR188273_chrX_rg.bam
samtools head eg/ERR188273_chrX_rg2.bam
```

Popular alignment tools such as BWA MEM and STAR can add read groups; use the `-R` and `--outSAMattrRGline` parameters for the respective tool.

```
bwa mem \
  -M \
  -t ${thread} \
  -R "@RG\tID:${sample_name}\tSM:${sample}\tPL:${platform}" \
  ${fasta} \
  ${fastq1} \
  ${fastq2} |
  samtools sort -@ ${thread} -O BAM |\
  tee ${sample_name}.bam |\
  samtools index - ${sample_name}.bam.bai

STAR \
  --runMode alignReads \
  --genomeDir ${star_index} \
  --readFilesIn ${fastq1} ${fastq2} \
  --readFilesCommand "gunzip -c" \
  --outFileNamePrefix ${prefix}. \
  --outSAMtype BAM Unsorted \
  --twopassMode Basic \
  --outSAMattrRGline ID:${id} PL:Illumina PU:${pu} LB:${lb} PI:0 SM:${sm} \
  --outSAMattributes NH HI AS nM NM ch \
  --runThreadN ${num_threads}
```

## Interpreting the BAM flags

The second column in a SAM/BAM file is the flag column; use the `flags` subcommand to understand specific flags. They may seem confusing at first but the encoding allows details about a read to be stored by just using a few digits. The trick is to convert the numerical digit into binary, and then use the table to interpret the binary numbers, where 1 = true and 0 = false. I wrote a blog post on BAM flags at <http://davetang.org/muse/2014/03/06/understanding-bam-flags/>.

```{bash engine.opts='-l'}
samtools flags
```

Find out about a `73` flag.

```{bash engine.opts='-l'}
samtools flags 73
```

### Proper pair

Reads that are properly paired are mapped within an expected distance with each other and with one pair in the reverse complement orientation. The script `generate_random_seq.pl` can generate reads that originate from different references and are thus discordant and not properly paired (as well as properly paired reads). In the example below, 10% of reads are not properly paired (set with `-d 0.1`).

```{bash engine.opts='-l'}
script/generate_random_seq.pl 30 10000 1984 > test_ref.fa
script/random_paired_end.pl -f test_ref.fa -l 100 -n 10000 -m 300 -d 0.1
bwa index test_ref.fa 2> /dev/null
bwa mem test_ref.fa l100_n10000_d300_1984_1.fq.gz l100_n10000_d300_1984_2.fq.gz > aln.sam 2> /dev/null
```

`samtools flagstat` will indicate that some reads (about 10%) mapped to different chromosomes.

```{bash engine.opts='-l'}
samtools flagstat aln.sam
```

Flag of a proper pair.

```{bash engine.opts='-l'}
samtools flag $(samtools view -f 2 aln.sam | head -1 | cut -f2)
```

Flag of a pair (that is not a proper pair).

```{bash engine.opts='-l'}
samtools flag $(samtools view -F 2 aln.sam | head -1 | cut -f2)
```


## Filtering unmapped reads

Use `-F 4` to filter out unmapped reads.

```{bash engine.opts='-l'}
samtools view -F 4 -b eg/ERR188273_chrX.bam > eg/ERR188273_chrX.mapped.bam
```

Use `-f 4` to keep only unmapped reads.

```{bash engine.opts='-l'}
samtools view -f 4 -b eg/ERR188273_chrX.bam > eg/ERR188273_chrX.unmapped.bam
```

We can use the `flags` subcommand to confirm that a value of four represents an unmapped read.

```{bash engine.opts='-l'}
samtools flags 4
```

## Extracting entries mapping to a specific loci

Use `samtools view` and the `ref:start-end` syntax to extract reads mapping within a specific genomic loci; this requires a BAM index file.

```{bash engine.opts='-l'}
samtools view eg/ERR188273_chrX.bam chrX:20000-30000
```

Note that this takes into account the mapping of the entire read and not just the starting position. For example, if you specified chrX:20000-30000, a 75 bp long read that starts its mapping from position 19999 will also be returned. In addition, you can save the output as another BAM file if you want.

```{bash engine.opts='-l'}
samtools view -b eg/ERR188273_chrX.bam chrX:20000-30000 > eg/ERR188273_chrX_20000_30000.bam
```

If you want reads mapped to a single reference (e.g. chromosome), just specify the `ref` and leave out the start and end values.

```{bash engine.opts='-l'}
samtools view eg/ERR188273_chrX.bam chrX | head
```

You can also use a BED file, with several entries, to extract reads of interest.

```{bash engine.opts='-l'}
cat eg/my.bed 

samtools view -L eg/my.bed eg/ERR188273_chrX.bam
```

## Extracting only the first read from paired end BAM files

Sometimes you only want the first pair of a mate. 0x0040 is hexadecimal for 64 (i.e. 16 * 4), which is binary for 1000000, corresponding to the read in the first read pair.

```{bash engine.opts='-l'}
samtools view -b -f 0x0040 eg/ERR188273_chrX.bam > eg/first.bam
```

Once again, you can use `flags` to verify this (it also accepts hexadecimal input).

```{bash engine.opts='-l'}
samtools flags 0x0040
```

## Stats

For simple statistics use `samtools flagstat`.

```{bash engine.opts='-l'}
samtools flagstat eg/ERR188273_chrX.bam
```

For more stats, use `samtools stats`.

```{bash engine.opts='-l'}
samtools stats eg/ERR188273_chrX.bam | grep ^SN
```

## samtools calmd/fillmd

The `calmd` or `fillmd` tool is useful for visualising mismatches and insertions in an alignment of a read to a reference genome. The `-e` argument changes identical bases between the read and reference into `=`.

```{bash engine.opts='-l'}
samtools view -b eg/ERR188273_chrX.bam | samtools fillmd -e - genome/chrX.fa > eg/ERR188273_chrX_fillmd.bam

head eg/ERR188273_chrX_fillmd.bam
```

## Creating FASTQ files from a BAM file

Use the `fastq` tool to create FASTQ files from a BAM file. For paired-end reads, use `-1` and `-2` to create separate FASTA files.

```{bash engine.opts='-l'}
samtools fastq -1 eg/ERR188273_chrX_1.fq -2 eg/ERR188273_chrX_2.fq eg/ERR188273_chrX.bam
head eg/ERR188273_chrX_1.fq
```

## Random subsampling of BAM file

The SAMtools view `-s` parameter allows you to randomly sub-sample a BAM file. Using `-s 0.5` will create a new BAM file with a random half of all mapped reads; unmapped reads are not sampled.

```{bash engine.opts='-l'}
samtools view -s 0.5 -b eg/ERR188273_chrX.bam > eg/ERR188273_chrX_rand.bam
```

## Count number of reads

Use `samtools idxstats` to print stats on a BAM file; this requires an index file which is created by running `samtools index`. The output of idxstats is a file with four tab-delimited columns:

1. Reference name
2. Sequence length of reference
3. Number of mapped reads
4. Number of unmapped reads

```{bash engine.opts='-l'}
samtools idxstats eg/ERR188273_chrX.bam
```

We can use this with `awk` to calculate:

The number of mapped reads by summing the third column.

```{bash engine.opts='-l'}
samtools idxstats eg/ERR188273_chrX.bam  | awk '{s+=$3} END {print s}'
```

The number of reads, which is the sum of mapped and unmapped reads.

```{bash engine.opts='-l'}
samtools idxstats eg/ERR188273_chrX.bam | awk '{s+=$3+$4} END {print s}'
```

## Obtaining genomic sequence

Use `faidx` to fetch genomic sequence; coordinates are 1-based.

We need to first index the reference FASTA file that was used to map the reads.

```{bash engine.opts='-l'}
samtools faidx genome/chrX.fa
```

Now we can obtain the sequence.

```{bash engine.opts='-l'}
samtools faidx genome/chrX.fa chrX:300000-300100
```

## Comparing BAM files

The output from `mpileup` can be used to compare BAM files. The commands below generates alignments using `bwa` and `minimap2`.

```{bash engine.opts='-l'}
len=100
n=10000
m=300
script/generate_random_seq.pl 30 1000000 1984 > test_ref.fa
script/random_paired_end.pl -f test_ref.fa -l ${len} -n ${n} -m ${m}
bwa index test_ref.fa 2> /dev/null

bwa mem test_ref.fa l${len}_n${n}_d${m}_1984_1.fq.gz l${len}_n${n}_d${m}_1984_2.fq.gz 2> /dev/null | samtools sort - -o aln_bwa.bam
minimap2 -ax sr test_ref.fa l${len}_n${n}_d${m}_1984_1.fq.gz l${len}_n${n}_d${m}_1984_2.fq.gz 2> /dev/null | samtools sort - -o aln_mm.bam
```

The BAM files can be used with `mpileup` to compare the depths.

```{bash engine.opts='-l'}
samtools mpileup -s -f test_ref.fa aln_bwa.bam aln_mm.bam | head -20
```

Another approach is to use [deepTools](https://deeptools.readthedocs.io/en/develop/) and the [bamCompare](https://deeptools.readthedocs.io/en/develop/content/tools/bamCompare.html) command. The bigWig output file shows the ratio of reads between `b1` and `b2` in 50 bp (default) windows.

## Converting reference names

One of the most annoying bioinformatics problems is the use of different chromosome names, e.g. chr1 vs 1, in different references even when the sequences are identical. The GRCh38 reference downloaded from Ensembl has chromosome names without the `chr`:

    >1 dna:chromosome chromosome:GRCh38:1:1:248956422:1 REF

Whereas the reference names from UCSC has the `chr`:

    >chr1  AC:CM000663.2  gi:568336023  LN:248956422  rl:Chromosome  M5:6aef897c3d6ff0c78aff06ac189178dd  AS:GRCh38

Luckily you can change the reference names using `samtools reheader` but just make sure your reference sequences are actually identical.

```{bash engine.opts='-l'}
samtools view eg/ERR188273_chrX.bam | head -2
```

View header

```{bash engine.opts='-l'}
samtools view -H eg/ERR188273_chrX.bam
```

Substitute header with new name.

```{bash engine.opts='-l'}
samtools view -H eg/ERR188273_chrX.bam | sed 's/SN:chrX/SN:X/' > eg/my_header
```

Save bam file with new ref and check it out.

```{bash engine.opts='-l'}
samtools reheader eg/my_header eg/ERR188273_chrX.bam > eg/ERR188273_X.bam
samtools view eg/ERR188273_X.bam | head -2
```

## Coverage

Coverage can mean the:

1. average depth of each covered base
2. percentage of bases covered

`samtools depth` and `samtools mpileup` can be used to indicate the depth of
each covered base (and used to calculate the average depth. `samtools coverage`
will provide both the average depth and percentage of bases covered per
chromosome/reference sequence.

`samtools depth` will return three columns: reference, position, and coverage.

```{bash engine.opts='-l'}
samtools depth -@ 4 eg/ERR188273_chrX.bam > ERR188273_depth.tsv
head ERR188273_depth.tsv
```

The average depth can be calculated by summing the third column and dividing by the total number of bases (be sure to use `-a` with `samtools depth` as that will output all positions including zero depth).

```{bash engine.opts='-l'}
samtools depth -@ 4 -a eg/ERR188273_chrX.bam | perl -ane '$t += $F[2]; END {$cov = $t / $.; printf "Bases covered:\t%.3f\nCoverage:\t%.3f\n", $., $cov}'
```

The `samtools mpileup` command also provides depth information (but not for reads that have a mapping quality of 0, by default) with some additional information:

1. Sequence name
2. 1-based coordinate
3. Reference base (when used with `-f`)
4. Number of reads covering this position
5. Read bases
6. Base qualities
7. Alignment mapping qualities (when used with `-s`)

```{bash engine.opts='-l'}
samtools mpileup -f genome/chrX.fa -s eg/ERR188273_chrX.bam > ERR188273_mpileup.tsv
head ERR188273_mpileup.tsv
```

Note that the start of the `samtools mpileup` output differ from the start of the `samtools depth` output. This is because `mpileup` performs some filtering by default. In the case of this example, read pairs that are not both mapped will be ignored. To count these "orphan" reads, use the `--count-orphans` argument.

```{bash engine.opts='-l'}
samtools mpileup -f genome/chrX.fa --count-orphans -s eg/ERR188273_chrX.bam > ERR188273_mpileup_orphans.tsv
head ERR188273_mpileup_orphans.tsv
```

In addition `mpileup` performs "per-Base Alignment Quality" (BAQ) by default and will adjust base quality scores. The default behaviour to to skip bases with baseQ/BAQ smaller than 13. If you are finding discrepancies between `mpileup`'s coverage calculation with another coverage tool, you can either set `--min-BQ` to `0` or use `--no-BAQ` to disable BAQ.

I have an [old blog post](https://davetang.org/muse/2015/08/26/samtools-mpileup/) on using `mpileup`.

`samtools coverage` will provide the following coverage statistics:

1. `rname` - Reference name / chromosome
2. `startpos` - Start position
3. `endpos` - End position (or sequence length)
4. `numreads` - Number reads aligned to the region (after filtering)
5. `covbases` - Number of covered bases with depth >= 1
6. `coverage` - Proportion of covered bases [0..1]
7. `meandepth` - Mean depth of coverage
8. `meanbaseq` - Mean base quality in covered region
9. `meanmapq` - Mean mapping quality of selected reads

```{bash engine.opts='-l'}
samtools coverage eg/ERR188273_chrX.bam
```

The example BAM file only contains reads for `chrX` hence the statistics are only returned for `chrX`.

Returning to our coverage definition at the start of this section:

1. average depth of each covered base = `meandepth`
2. percentage of bases covered = `covbases`

The [mosdepth](https://github.com/brentp/mosdepth) tool can also calculate depth (and much faster than `samtools depth`) per base or within a given window. The output is given in a BED file, where the fourth column indicates the coverage.

```{bash engine.opts='-l'}
mosdepth ERR188273 eg/ERR188273_chrX.bam
gunzip -c ERR188273.per-base.bed.gz | head
```

`mosdepth` coverage.

```{bash engine.opts='-l'}
cat ERR188273.mosdepth.summary.txt
```

Coverage in using a 500 bp window.

```{bash engine.opts='-l'}
mosdepth -n --fast-mode --by 500 ERR188273_500 eg/ERR188273_chrX.bam
gunzip -c ERR188273_500.regions.bed.gz | head
```

## Stargazers over time

[![Stargazers over time](https://starchart.cc/davetang/learning_bam_file.svg)](https://starchart.cc/davetang/learning_bam_file)