Short Read Mapping

Sam, Bam File Formats

SCP

Can I use Blast or Diamond?

• Alignment methods must have tradeoffs for speed vs accuracy
• Depending on the application, may want to make different tradeoffs
• Different types of alignment objectives lead to different categories of aligners

Short Read Mappers

• BLAST is much faster than original algorithms (Smith Waterman) but still too slow for the amount of data produced by NGS technology
• Resequencing usually involves comparing very similar sequences (>90% identity in residues) to a reference genome
• Nucleotides only
• Software that leverages this high percent identity can be faster

ACGTTGCTAGCGGTCAGTAATCGGCCCGTATCGATATATGGCG

Reference Genome

Read 1 Read 2

CAATCGCCCC ACGTTCTA

​

? ?

Heinz 1706

Cherokee purple

Where in the genome did each of the short reads originate?

ACGTTGCTAGCGGTCAGTAATCGGCCCGTATCGATATATGGCG

Reference Genome

Read 2 Read 1

ACGTT CTA CAATCGCCCC

||||| ||| ||||| |||

​

Short Read Mappers

• Orders of magnitude faster than BLAST
• tens to hundreds of millions of reads mapped per hour per CPU
• Only matches of 95% identity or greater are found
• Usually only output the best hit or the set of all equivalently good hits
• The point is usually to find the origin in the reference genome
• Other genomic regions of lower identity are not considered useful

Bwa-mem2

• Use for any sequences greater than 70bp up to 1Mb
• Faster (1.3-3.1X) than bwa mem but needs MUCH more RAM
• Will work well with reads with
• 2% error/differences for a 100bp read
• 3% error/differences for a 200bp read
• Individual chromosomes cannot be longer than 2GB

​

Ambystoma mexicanum Paris japonica Pinus taeda

32Gb 150Gb 22Gb

Indexing Strategies

• The first step is to transform part of the data into a more suitable form for fast searching
• Indexing – creating a glossary or look up table

Algorithms

BWA mem

1. BWT = Burrows-Wheeler transform
1. Excellent compression - clusters repeated characters together
2. Can be searched without decompressing
3. Very fast, very efficient for highly similar strings
1. FM-index
1. Maximal exact matches only
2. Indels and mismatches are handled during extension of the alignment

​

Unique Read Mapping

• Some reads can be mapped uniquely to the reference
• Some map to multiple locations
• How to deal with multi mapped reads?

​

? ? ?

The read matches to each location equally well.

​

Clever Ways to find “Best” Alignments

• Use the quality values
• Penalize mismatches at high quality bases more than mismatches at low quality bases
• Paired End information
• If one read does not map uniquely, but the other does, use that information to place the non-unique one

How BWA-MEM2 Processes Multi-Mapping Reads

​

Scenario 1. The full read maps to one or more locations.

​

SECONDARY

Alignments

Scenario 2. Part of the read maps to one location, part of the read maps to another.

​

SUPPLEMENTARY

Alignments

Ref genome

Ref genome

read

read

How BWA-MEM2 Processes Multi-Mapping Reads

​

Scenario 1. The full read maps to one or more locations.

​

SECONDARY Alignments

• BWA-MEM picks one location at random to be the primary alignment.
• It assigns a MAPQ score of 0 to indicate it is not uniquely mapped.
• Up to 5 alternative equal-scoring locations are listed in the XA:Z tag.

​

Bwa-mem2 outputs NO secondary alignments by default; you only get flag-256 records with -a.

Scenario 2. Part of the read maps to one location, part of the read maps to another.

​

SUPPLEMENTARY Alignments

• one part becomes the primary alignment.
• The other parts are flagged as supplementary alignments (flag 2048).

​

Did It Actually Work? Mapping QC

samtools flagstat aln.bam # % mapped, % properly paired,

duplicates

samtools stats aln.bam # insert size, error rate, coverage

(plot-bamstats)

​

Red flags:

• <90% mapped → wrong reference, contamination, or untrimmed adapters
• Low properly-paired % or a bimodal insert size → library or assembly problem
• High duplicate rate → over-amplified, low-complexity library

Read Groups: the @RG Line You Need for GATK and Picard Tools

• A read group stamps every read with where it came from

-R '@RG\tID:run1\tSM:sample1\tPL:ILLUMINA\tLB:lib1'

• ID = unique per run/lane SM = sample name (what variant callers group on)
• PL = platform LB = library (used by duplicate marking)
• GATK requires read groups
• Bcftools and deep variant do not
• Adding them later is possible but slow - must rewriting the entire BAM (samtools addreplacerg)

When to use bwa-mem2

• Local alignment: ends that do not match are soft-clipped (S), not forced to align
• Contrast: bowtie2 defaults to end-to-end — the whole read must align
• Soft clipping is a feature: adapter remnants, chimeras, SV breakpoints

​

​

​

​

Choosing a mapper:

• minimap2 — long/noisy reads (PacBio, Oxford Nanopore) and whole-genome alignment
• HISAT2 / STAR — RNA-seq against a genome
• Do not use bwa-mem2 for RNA-seq or long reads

Caution: Reference Mismatches

• chr1 (UCSC) vs 1 (Ensembl) — same genome, different contig names
• Best case the tool errors; worst case it silently reports zero overlaps
• GRCh37 / hg19 / GRCh38 / hg38 coordinates are NOT interchangeable
• Your BAM, your annotation (GFF/GTF), and your VCF must all use the same assembly

​

PAUSE

Go get fastqc started….

III. SAM, BAM AND CRAM FORMAT

SAM Format

• SAM = Sequence Alignment/Map format
• Tab delimited plain text
• Store large nucleotide sequence alignments
• Alignment of every read
• Including gaps, SNPs and structural variants
• Pairing of reads
• Can record more than one alignment location in the genome
• Stores quality values

​

Formats we’ve already learned:

• fasta
• fastq

SAM Format

Strengths

• Flexible
• Useful for operations on very large sequences
• Extremely detailed documentation
• https://samtools.github.io/hts-specs/
• Manipulations can be done with the software samtools

​

SAM

Example sam with one read:

​

@SQ SN:gi|254160123|ref|NC_012967.1| LN:4629812

@PG ID:bwa PN:bwa VN:0.7.12-r1039 CL:/lustre/projects/rnaseq_ws/apps/bwa-0.7.12/bwa sampe ../raw_data/NC_012967.1.fasta aln_SRR030257_1.sai aln_SRR030257_2.sai ../raw_data/SRR030257_1.fastq ../raw_data/SRR030257_2.fastq

SRR030257.1 99 gi|254160123|ref|NC_012967.1| 950180 60 36M = 950295 151 TTACACTCCTGTTAATCCATACAGCAACAGTATTGG AAA;A;AA?A?AAAAA?;?A?1A;;????566)=*1 XT:A:U NM:i:1 SM:i:37 AM:i:25 X0:i:1 X1:i:0 XM:i:1 XO:i:0 XG:i:0 MD:Z:32C3

​

SAM - Header

• Structure
• Optional Header at top of file
• Alignment information

@

@

@

@

Read Info

Read Info

Read Info

​

SAM - Header

• Header lines start with @ symbol
• Always at top of file
• Contain lots of information about what was mapped, what it was mapped to, and how (metadata)
• the version information for the SAM/BAM file
• whether or not and how the file is sorted
• information about the reference sequences
• any processing that was used to generate the various reads in the file
• software version

​

Simple Header

@HD = first line

VN = version of SAM format

SO = sort order (this is sorted by coordinates)

​

@HD VN:1.5 SO:coordinate

@SQ SN:ref LN:45

​

@SQ = reference sequence

SN = Sequence reference Name

LN = seqeunce reference length

Decipher header information: https://samtools.github.io/hts-specs/SAMv1.pdf

Alignment Line

Below the headers are the alignment records

Tab-delimited fields

​

1 QNAME Query template/pair NAME

2 FLAG bitwise FLAG

3 RNAME Reference sequence NAME

4 POS 1-based leftmost POSition/coordinate of clipped sequence

5 MAPQ MAPping Quality (Phred-scaled)

6 CIGAR extended CIGAR string

7 MRNM Mate Reference sequence NaMe (`=' if same as RNAME)

8 MPOS 1-based Mate POSistion

9 TLEN inferred Template LENgth (insert size)

10 SEQ query SEQuence on the same strand as the reference

11 QUAL query QUALity (ASCII-33 gives the Phred base quality)

12+ OPT variable OPTional fields in the format TAG:VTYPE:VALUE

​

Lets unpack this alignment line, taken from a SAM file. Tab delimited:

SRR030257.2000020 83 gi|254160123|ref|NC_012967.1| 3295752 60 36M = 3295706 -82 TGCTGGCGGCGATATCGTCCGTGGTTCCGATCTGGT ?%<91<?>>??AAAAAAAAAAAAAAAAAAAAAAAAA XT:A: NM:i:0 SM:i:37 AM:i:37 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:36

SAM Field 1

Query name = Read name

​

SRR030257.2000020

​

Field 2:

Flag

​

83

​

​

64 + 16 + 2 + 1

​

1 = Read is paired

2 = Read mapped in proper pair

16 = Read mapped to reverse strand

64 = First in pair

Look up a SAM flag: https://broadinstitute.github.io/picard/explain-flags.html

Online tools can help decipher the flag:

https://broadinstitute.github.io/picard/explain-flags.html

SAM Field 3

Reference sequence name (useful especially if you have multiple chromosomes)

​

gi|254160123|ref|NC_012967.1|

​

SAM Field 4

Position- 1-based leftmost mapping position of the first matching base

​

3295752

SAM Field 5

Mapping Quality, Phred-scaled. bwa caps at 60; the SAM spec allows 0–254 (255 = unavailable)

• equals −10 log₁₀ Pr{mapping position is wrong}, rounded to the nearest integer
• Same as phred
• MAPQ 30 = 99.9% chance the position is right (1 in 1,000 wrong)
• MAPQ of 0 = Maps equally well to multiple locations
• MAPQ of 1-3 = read matches a small number of alternative positions

​

60 0.0001% probability wrong (1 in 1,000,000)

​

•  MAPQ doesn't measure how well a read aligned — it measures how much better its best alignment is than its best alternative.

SAM Field 6

CIGAR String

​

​

​

36M 36 bases consume query AND reference

8S28M 8 bases soft-clipped, 28 aligned

bwa mem only ever emits M — never = or X — so mismatches hide inside M

M does NOT mean "perfect match" — it means alignment match and INCLUDES mismatches.

• A read with 5 SNPs still reports 36M.

More CIGAR

Position: 5

CIGAR: 3M1I3M1D5M

The Full CIGAR Operation Table

Op Meaning Consumes

M alignment match (match OR mismatch) query + ref

I insertion relative to the reference query

D deletion from the reference ref

N skipped region (intron in RNA-seq) ref

S soft clip — bases REMAIN in SEQ query

H hard clip — bases REMOVED from SEQ neither

P padding

= true match

X mismatch

• bwa mem only ever emits M — never = or X — so mismatches hide inside M
• S vs H matters: supplementary records often use H, so SEQ looks truncated
• M + I + S + = + X must sum to the read length

SAM Field 7

Reference sequence for the next read in the template

• For a forward read, this is the reference where the reverse read maps
• For a reverse read, this is the reference where the forward read maps

​

= reverse read maps on the same reference

SAM Field 8

Position where the next read maps

​

3295706

​

​

SAM Field 9

Observed template length

From leftmost base to rightmost base

Negative if this read is the rightmost read

-82

​

R1

R2

Template length

SAM Field 10

Sequence of the read

​

TGCTGGCGGCGATATCGTCCGTGGTTCCGATCTGGT

SAM Field 11

Quality of the read

​

?%<91<?>>??AAAAAAAAAAAAAAAAAAAAAAAAA

SAM Field 12

Optional MORE information

​

TAG:TYPE:VALUE format

​

XT:A:U NM:i:0 SM:i:37 AM:i:37 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:36

​

Anything with an X is specified by the user or by the mapping software, and is not part of the SAM spec.

Tags from bwa-mem2

NM:i: edit distance to the reference (mismatches + indel bases)

MD:Z: which reference bases mismatched — reconstruct ref without loading it

AS:i: alignment score of THIS alignment

XS:i: score of the best OTHER alignment

XA:Z: up to 5 alternative hits, as chr,pos,CIGAR,NM;

SA:Z: pointer to supplementary alignments (the split reads, flag 2048)

MC:Z / MQ:i mate’s CIGAR and MAPQ RG:Z read group

• MAPQ relies on the difference between AS and XS
• AS ≈ XS drives MAPQ toward 0 — the read fits two places equally well

BAM Format

• Sister format to SAM
• BAM – Binary version of SAM
• compressed BGZF (Blocked GNU Zip Format) - a variant of GZIP (GNU ZIP),
• files are bigger than GZIP files, but they are much faster for random access
• Can index and then look up information embedded in the file without decompressing the whole file
• Up to 75% smaller in size - Always store read alignments as bam not sam!
• Not readable by people

CRAM

• Introduced in 2011 by EMBL/EBI
• Even smaller and more efficient than BAM files
• Reference-based compression: CRAM stores only the DIFFERENCES from the reference
• You MUST have the exact reference FASTA to decode a CRAM

Fritz, Markus Hsi-Yang, et al. "Efficient storage of high throughput DNA sequencing data using reference-based compression." Genome research 21.5 (2011): 734-740.

Moving Data

• The best way to move data to/from your laptop or lab computer to ISAAC/Sphinx is with secure copy

SCP

• Secure copy
• Command line interface based on SSH

​

Remember our previous cp command:

cp <location of what we want> <location where it should go>

​

cp from to

​

​

SCP

• Secure copy
• Command line interface based on SSH

From a local directory to Sphinx:

scp file.txt mstaton1@sphinx.ag.utk.edu:~

​

​

​

From Sphinx to your local directory:

scp mstaton1@sphinx.ag.utk.edu:~/file.txt .

​

Local file

Username

Remote Machine

Remote location to place file

Username

Remote Machine

Remote Location

Local Location