...
Tip | ||
---|---|---|
| ||
Use our summer school reservation (CoreNGSday5CoreNGS-Fri) when submitting batch jobs to get higher priority on the ls6 normal queue today:
|
Table of Contents |
---|
Overview
...
Code Block | ||||
---|---|---|---|---|
| ||||
idev -m 180 -N 1 -A OTH21164 -r CoreNGSday5 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 a renamed version of the yeast_pairedendpe.sort.bam file from our The Basic Alignment workflow Workflow exercises.
Code Block | ||||
---|---|---|---|---|
| ||||
mkdir -p $SCRATCH/core_ngs/samtools cd $SCRATCH/core_ngs/samtools cp $CORENGS/catchup/yeastfor_stuffsamtools/*.bam* . |
References
Handy links
- The SAM format specification
- especially section 1.4 - alignment section fields
- Manual for SAMTools
- especially the 1st section on samtools view.
...
Code Block | ||
---|---|---|
| ||
module load biocontainers # takes a while module load samtools cd $SCRATCH/core_ngs/samtools samtools view yeast_chip_pe.sort.bam | head |
With all the line wrapping, it looks pretty ugly. Still, you can pick out the CIGAR string in column 6. Let's select just that column with cut:
Code Block | ||
---|---|---|
| ||
samtools view yeast_chip_pe.sort.bam | cut -f 6 | head -20 |
...
Code Block | ||
---|---|---|
| ||
samtools view -F 0x4 yeast_chip_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) to be sure so that this character class syntax is recognized.
Code Block | ||
---|---|---|
| ||
samtools view -F 0x4 yeast_chip_pe.sort.bam | cut -f 6 | grep -P '[ID]' | head |
...
Code Block | ||||
---|---|---|---|---|
| ||||
samtools view -F 0x4 yeast_chip_pe.sort.bam | cut -f 6 | grep -P '[ID]' | wc -l |
...
Code Block | ||||
---|---|---|---|---|
| ||||
samtools view -c -F 0x4 yeast_chip_pe.sort.bam |
There should be 547664 mapped alignments.
...
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) 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).
Code Block | ||||
---|---|---|---|---|
| ||||
echo "`samtools view -F 0x4 -c yeast_chip_pe.sort.bam` `samtools \ $( samtools view -F 0x4 yeast_chip_pe.sort.bam | cut -f 6 | grep -P '[ID]' | wc -l`l )" \ | awk '{ print 100 * $2/$1,"%" }' |
...
Tip | |||||
---|---|---|---|---|---|
awk also has a printf function, which can take the standard formatting commands (see https://en.wikipedia.org/wiki/Printf_format_string#Type_field).
Notes:
|
...
Code Block | ||
---|---|---|
| ||
# 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}' |
...
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.
...
- Alternate locations for a mapped read are are 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)
...
- 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).
- its definition of a secondary alignment is different (and a bit non-standard)
- 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.
...
Expand | |||||||
---|---|---|---|---|---|---|---|
| |||||||
Note that we use the -b flag to tell samtools view to output BAM format (not the SAM format text we've been looking at).
|
...
Expand | |||||||
---|---|---|---|---|---|---|---|
| |||||||
There were 1,184,360 alignment records in the original BAM, and only 456,890 in the quality-filtered BAM, around 38% of our starting reads. Since we have only properly paired reads, the filtered BAM will contain equal numbers of both R1s and R2s. So the number of read pairs is 456890/2 or 228451. |
...