Page History

Overview

...

The most common samtools view filtering options are:

• -q N – only report alignment records with mapping quality of at least N (>= N).
• -f 0xXX – only report alignment records where the specified flags XX are all set (are all 1)
  • you can provide the flags in decimal, or as here as hexidecimalhexadecimal
• -F 0xXX – only report alignment records where the specified flags XX are all cleared (are all 0)

Setup

Login to login5.ls5ls6.tacc.utexas.edu, set up a directory for these exercises, and copy an indexed BAM file there. This is a renamed version of the yeast_pairedend.sort.bam file from our Alignment workflow exercises and start an idev session.

mkdir -p $SCRATCH/core_ngs/samtools
cd $SCRATCH/core_ngs/samtools
cp $CORENGS/catchup/for_samtools/* .

References

Handy links

• The SAM format specification
  
  • especially section 1.4 - alignment section fields
• Manual for SAMTools
  
  • especially the 1st section on samtools view.

SAM header fields

The 11 SAM alignment record required fields (Tab-delimited).

Image Removed

SAM flags field

idev -m 180 -N 1 -A OTH21164 -r CoreNGS-Fri
# or
idev -m 90 -N 1 -A OTH21164 -p development

Next set up a directory for these exercises, and copy an indexed BAM file there. This is the yeast_pe.sort.bam file from our The Basic Alignment Workflow exercises.

mkdir -p $SCRATCH/core_ngs/samtools
cd $SCRATCH/core_ngs/samtools
cp $CORENGS/catchup/for_samtools/* .

References

Handy links

• The SAM format specification
  
  • especially section 1.4 - alignment section fields
• Manual for SAMTools
  
  • especially the 1st section on samtools view.

SAM header fields

The 11 SAM alignment record required fields (Tab-delimited).

Image Added

SAM flags field

Meaning of each bit Meaning of each bit (flag) in the SAM alignment records flags field (column 2). The most commonly querried flags are denoted in red.

Image RemovedImage Added

SAM CIGAR string

...

We'll do this first exercise one step at a time to get comfortable with piping commands. Start by just looking at the first few alignment records of the BAM file:

cd $SCRATCH/core_ngs/alignment/samtools
module load biocontainers  # takes a while
module load samtools

cd $SCRATCH/core_ngs/samtools
samtools view yeast_pe.sort.bam | head

With all the line wrapping, it looks pretty ugly. Still, you can pick out the CIGAR string in colum column 6. Let's select just that column with cut:

samtools view yeast_pe.sort.bam | cut -f 6 | head -20

Next, make sure we're only looking at alignment records that represent mapped reads. The -F 0x4 option says to filter records where the 0x4 flag (read unmapped) is 0, resulting it only mapped reads being output.

samtools view -F 0x4 yeast_pe.sort.bam | cut -f 6 | head -20

Now we use grep with the pattern '[ID]' to select lines (which are now just CIGAR strings) that have Insert or Deletion operators. The brackets ( [ ] ) denote a character class pattern, which matches any of the characters inside the brackets. Be sure to ask for Perl regular expressions (-P) so that this character class syntax is recognized.

samtools view -F 0x4 yeast_pe.sort.bam | cut -f 6 | grep -P '[ID]' | head 

Ok, this looks good. We're ready to run this against all the alignment records, and count the results:

...

In addition to the rate, converted to a percentage, we also output the literal percent sign ( % ). The double quotes ( " ) denote a literal string in awk, and the comma between the number and the string says to separate the two fields with whatever the default Output Field Separator (OFS) is. By default, both awk's Input Field Separator (FS) and Output Field Separator are is whitespace (any number of space and Tab characters) and its Output Field Separator is a single space.

So what if we want to get fancy and do all this in a one-liner command? We can use echo and backtick evaluation to put both numbers in a string, then pipe that string to awk. This time we use awk body code, and refer to the two whitespace-separated fields being passed in: total count ($1) and indel count ($2).

echo "`samtools view -F 0x4 -c yeast_pe.sort.bam` \
 `samtools view   $( samtools view -F 0x4 yeast_pe.sort.bam | cut -f 6 | grep -P '[ID]' | wc -l`"l )" \
    | awk '{ print 100 * $2/$1,"%" }'

Filtering by location range

Sometimes you just want to examine a subset of reads in detail. Once you have a sorted and indexed BAM, you can use the coordinate filtering options of samtools view to do this. Here are some examples:

mkdir -p $SCRATCH/core_ngs/samtools
cd $SCRATCH/core_ngs/samtools
cp $CORENGS/catchup/for_samtools/* .
module load samtools

awk 'BEGIN{ printf("%.2f %%\n", 100*6697/547664) }'

Filtering by location range

Sometimes you just want to examine a subset of reads in detail. Once you have a sorted and indexed BAM, you can use the coordinate filtering options of samtools view to do this. Here are some examples:

mkdir -p $SCRATCH/core_ngs/samtools
cd $SCRATCH/core_ngs/samtools
cp $CORENGS/catchup/for_samtools/* .
module load biocontainers
module load samtools

cd $SCRATCH/core_ngs/samtools

# count the 

cd $SCRATCH/core_ngs/alignment/samtools

# count the number of reads mapped to chromosome 2 (chrII)
samtools view -c -F 0x4 yeast_pe.sort.bam chrII

# count the number of reads mapped to chromosomes 1 or M (chrI, chrM)
samtools view -c -F 0x4 yeast_pe.sort.bam chrI chrM

# count the number of reads mapped to chromosomeschromosome 1 that overlap coordinates 1000-20002 (chrII)
samtools view -c -F 0x4 yeast_pe.sort.bam chrI:1000-2000
chrII

# sincecount therethe arenumber onlyof 20reads readsmapped into the chrI:1000-2000 region, examine them individuallychromosomes 1 or M (chrI, chrM)
samtools view -c -F 0x4 yeast_pe.sort.bam chrI:1000-2000 chrM

# lookcount atthe anumber subsetof ofreads fieldmapped for chrI:1000-2000 reads
#   2=flags, 3=contig, 4=start, 5=mapping quality, 6=CIGAR, 9=insert sizeto chromosomes 1 that overlap coordinates 1000-2000
samtools view -c -F 0x4 yeast_pe.sort.bam chrI:1000-2000

# since |there cutare -f 2-6,9

Since you will find yourself wanting to interpret the flags field, and it's easier to do that when it is represented as hexadecimal, let's use awk to help do that. 

# look at a subset of field for chrI:1000-only 20 reads in the chrI:1000-2000 region, examine them individually
samtools view -F 0x4 yeast_pe.sort.bam chrI:1000-2000

# look at a subset of field for chrI:1000-2000 reads
#   2=flags, 3=contig, 4=start, 5=mapping quality, 6=CIGAR, 9=insert size
samtools view -F 0x4 yeast_pe.sort.bam chrI:1000-2000 | cut -f 2-6,9 | awk '
  BEGIN{FS=OFS="\t"}{$1 = sprintf("0x%x", $1); print}'

Notes:

• Usually, if a command line is continued on more than one line, the line continuation character ( \ ) is used.
  • here we don't need a line continuation character because the linefeed after the first single quote is part of the script.
• There is a sprintf (string print formatted) function in awk, just as there is in many modern programming languages.
• We use sprintf to re-format the 1st input field (here the cut flags) in hex
  
  • using the format directive (%x)
  • also adding a literal "0x" prefix to denote the base is hexadecimal
• The results of sprintf are stored back into the 1st field ($1), replacing the original value.
• awk's print statement is then used with no other arguments to print all its input fields.

Exercise: How many of the chrI:1000-2000 alignments are from properly paired mapped reads?

samtools view -F 0x4 yeast_pe.sort.bam chrI:1000-2000 | awk '
  BEGIN{FS=OFS="\t"}
  {$2 = sprintf("0x%x", $2); print $2}' | grep -c '3$'

Since you will find yourself wanting to interpret the flags field, and it's easier to do that when it is represented as hexadecimal, let's use awk to help do that. 

# look at a subset of field for chrI:1000-2000 reads
#   2=flags, 3=contig, 4=start, 5=mapping quality, 6=CIGAR, 9=insert size
samtools view -F 0x4 yeast_pe.sort.bam chrI:1000-2000 | cut -f 2-6,9 | \
  awk 'BEGIN{FS=OFS="\t"}
       {$1 = sprintf("0x%x", $1); print}'

Notes:

• If a command line is continued on more than one line, the line continuation character ( \ ) is used.
• There is a sprintf (string print formatted) function in awk, just as there is in many modern programming languages.
• We use sprintf to re-format the 1st input field (here the cut flags) in hex
  
  • using the format directive (%x)
  • also adding a literal "0x" prefix to denote the numeric base is hexadecimal
• The results of sprintf are stored back into the 1st field ($1), replacing the original value.
• awk's print statement is then used with no other arguments to print all its input fields, including the changed $1.

Exercise: How many of the chrI:1000-2000 alignments are from properly paired mapped reads?

samtools view -F 0x4

samtools view -c -F 0x4 -f 0x2 yeast_pe.sort.bam chrI:1000-2000

About mapping quality

Mapping qualities are a measure of how likely a given sequence alignment to its reported location is correct. If a read's mapping quality is low (especially if it is zero, or mapQ 0 for short) the read maps to multiple locations on the genome (they are multi-hit reads), and we can't be sure whether the reported location is the correct one.

Aligners also differ in whether they report alternate alignments for multi-hit r reads. Some things to keep in mind:

• Alternate locations for a mapped read are are flagged as secondary (flag 0x100)
• While they often provide valuable information, secondary reads must be filtered for some downstream applications
  • e.g., ChIP-seq peak calling and variant analysis with GATK.
• When secondary reads are reported, the total number of alignment records in the BAM file is greater than the number of reads in the input FASTQ files
  • this affects how the true mapping rate must be calculated
  • true mapping rate = (pirmary mapped reads) / (total BAM file sequences - secondary mapped reads)

Here are some examples of how different aligners handle reporting of multi-hit reads and their mapping qualities:

• bwa aln (global alignment) and bowtie2 with default parameters (both --local and --global) report at most one location for a read that maps
  
  • this will be the location with the best mapping quality and alignment
  • if a read maps equally well to multiple locations, these aligners pick one location at random
    • bwa aln will always report a 0 mapping quality for these multi-hit reads
    • bowtie2 will report a low mapping quality (< 10), based on the complexity (information content) of the sequence
• bwa mem (local alignment) can always report more than one location for a mapped read
  • for example, if one part of a read maps to one location and other part maps somewhere else (e.g. because of RNA splicing)
  • there is no way to disable reporting of secondary reads with bwa mem.
• bowtie2 can be configured to report more than one location for a mapped read
  
  • the -k N option says to report up to N alignments for a read
• most transcriptome-aware RNA-seq aligners by default can report more than one location for a mapped read
  • e.g. tophat2, hisat2, star.

Filtering high-quality reads

Using our yeast_pe.sort.bam file, let's do some some quality filtering.

mkdir -p $SCRATCH/core_ngs/samtools
cd $SCRATCH/core_ngs/samtools
cp $CORENGS/catchup/for_samtools/* .
module load samtools

Exercise:  Use samtools view with -F, -f and -q options to create a BAM containing only mapped, properly paired, high-quality (mapQ 20+) reads. Call the output file yeast_pe.sort.filt.bam.

# if needed... cd $SCRATCH/core_ngs/alignment/samtools module load samtools samtools view -F 0x04 -f 0x2 -q 20 -o yeast_pe.sort.filt.bam yeast_pe.sort.bam

Exercise: How many records are in the filtered BAM compared to the original? How many read pairs does this represent?

 | awk '
  BEGIN{FS=OFS="\t"}
  {$2 = sprintf("0x%x", $2); print $2}' | grep -c '3$'

samtools view -c -F 0x4 -f 0x2 yeast_pe.sort.bam chrI:1000-2000

About mapping quality

Mapping qualities are a measure of how likely a given sequence alignment to its reported location is correct. If a read's mapping quality is low (especially if it is zero, or mapQ 0 for short) the read maps to multiple locations on the genome (they are multi-hit or multi-mapping reads), and we can't be sure whether the reported location is the correct one.

Aligners also differ in whether they report alternate alignments for multi-hit reads. Some things to keep in mind:

• Alternate locations for a mapped read will be flagged as secondary (flag 0x100)
• While they often provide valuable information, secondary reads must be filtered for some downstream applications
  • e.g., ChIP-seq peak calling and variant analysis with GATK.
• When secondary reads are reported, the total number of alignment records in the BAM file is greater than the number of reads in the input FASTQ files!
  • this affects how the true mapping rate must be calculated
  • true mapping rate = (pirmary mapped reads) / (total BAM file sequences - secondary mapped reads)

Here are some examples of how different aligners handle reporting of multi-hit reads and their mapping qualities:

• bwa aln (global alignment) and bowtie2 with default parameters (both --local default end-to-end mode) report at most one location for a read that maps
  
  • this will be the location with the best mapping quality and alignment
  • if a given read maps equally well to multiple locations, these aligners pick one location at random
    • bwa aln will always report a 0 mapping quality for these multi-hit reads
    • bowtie2 will report a low mapping quality (< 10), based on the complexity (information content) of the sequence
• bwa mem (local alignment) can always report more than one location for a mapped read
  • its definition of a secondary alignment is different (and a bit non-standard)
    
    • if one part of a read maps to one location and another part maps somewhere else (e.g. because of RNA splicing), the longer alignment is marked as primary and the shorter as secondary.
  • there is no way to disable reporting of secondary alignments with bwa mem.
    • but they can be filtered from the sorted BAM with -F 0x100 (secondary alignment flag = 0).
• bowtie2 can be configured to report more than one location for a mapped read
  
  • the -k N option says to report up to N alignments for a read
• most transcriptome-aware RNA-seq aligners by default report more than one location for a mapped read
  • e.g. hisat2, star, tophat2.
  • when reads are quantified (counted with respect to genes), multiply-mapped reads can be counted fractionally
    • e.g. if a read maps to 5 genes, it can be counted as 1/5 for each of the genes

Filtering for high-quality reads

Using our yeast_pe.sort.bam file, let's do some some quality filtering.

mkdir -p $SCRATCH/core_ngs/samtools
cd $SCRATCH/core_ngs/samtools
cp $CORENGS/catchup/for_samtools/* .
module load biocontainers
module load samtools

Exercise:  Use samtools view with -F, -f and -q options to create a BAM containing only mapped, properly paired, high-quality (mapQ 20+) reads. Call the output file yeast_pe.sort.filt.bam.

samtools view -b -F 0x4 -f 0x2 -q 20 yeast_pe.sort.bam > yeast_pe.sort.filt.bam 
# or
samtools view -b -F 0x4 -f 0x2 -q 20 -o yeast_pe.sort.filt.bam  yeast_pe.sort.bam 

Exercise: How many records are in the filtered BAM compared to the original? How many read pairs does this represent?

samtools view -c yeast_pe.sort.bam
samtools view -c yeast_pe.sort.filt.bam

# or to be really fancy...
echo "`samtools view -c yeast_pe.sort.bam` \
      `samtools view -c yeast_pe.sort.filt.bam`" | \
       awk '{printf "%d original reads\n", $1;
             printf "%d Q20 reads (%d read pairs)\n", $2, $2/2;
             printf "%0.2f%% high-quality reads\n", 100*$2/$1}'

samtools view -c yeast_pe.sort.bam
samtools view -c yeast_pe.sort.filt.bam

Exercise: If our original BAM contained secondary reads, How many primary aligned reads (0x100 = 1) how would we exclude those also0) are in the bwa_local.sort.dup.bam file?

samtools view -b -F 0x104 -f 0x2 -q 20 -o yeast_pe.sort.filt.bam yeast_pe.sort.bam
#or
samtools view -b -F 0x104 -f 0x2 -q 20 yeast_pe.sort.bam > yeast_pe.sort.filt.bam -c bwa_local.sort.dup.bam

samtools view -F 0x4 -f 0x100 -c bwa_local.sort.dup.bam