Slides

Report
Previous Lecture: Next-Generation
DNA Sequencing Technology
This Lecture
NGS Alignment
Slides are from Stratos Efstathiadis, Cole Trapneli, Steven Salzberg, Ben Langmead,
Thomas Keane, Dennis Wall, Jianhua Ruan, J Fass
Learning Objectives
•
•
•
•
•
•
•
•
Challenge of NGS read alignment
Suffix Arrays
FM index
BWA algorithm
SAM/BAM format
CIGAR string for variants
Software to work directly with SAM/BAM files
Sequence viewers
Short Read Alignment
• Given a reference and a set of reads, report at least
one “good” local alignment for each read if one
exists
– Approximate answer to: where in genome did the read originate?
• What is “good”? For now, we concentrate on:
– Fewer mismatches is better
…TGATCATA…
GATCAA
– Failing to align a low-quality base
is better than failing to align a …TGATATTA…
GATcaT
high-quality base
– Match Uniqueness of the
alignment
better than
…TGATCATA…
GAGAAT
better than
…TGATcaTA…
GTACAT
Multiple mapping
• A single tag may occur more than once in the
reference genome.
• The user may choose to ignore tags that
appear more than n times.
• As n gets large, you get more data, but also
more noise in the data.
Inexact matching
?
• An observed tag may not exactly match any position in the reference
genome.
• Sometimes, the tag almost matches one or more positions.
• Such mismatches may represent a SNP (single-nucleotide polymorphism,
see wikipedia) or a bad read-out.
• The user can specify the maximum number of mismatches, or a phredstyle quality score threshold.
• As the number of allowed mismatches goes up, the number of mapped
tags increases, but so does the number of incorrectly mapped tags.
% of Paired K-mers with Uniquely
Assignable Location
Read Length is Not As Important For
Resequencing
100%
90%
80%
70%
60%
E.COLI
50%
HUMAN
40%
30%
20%
10%
0%
8
Jay Shendure
10
12
14 16
18
20
Length of K-mer Reads (bp)
New alignment algorithms must address the
requirements and characteristics of NGS reads
– Hundreds of Millions of reads per run (30x genome coverage)
– Short Reads (as short as 36bp)
– Different types of reads (single-end, paired-end, mate-pair,
etc.)
– Base-calling quality factors (should the aligner use them?)
– Sequencing errors ( ~ 1%)
– Repetitive regions
– Sequencing sample vs. reference genome
– Must adjust to evolving sequencing technologies and data
formats
Index: an auxiliary data structure
Two classes of indexing algorithms:
(1) Hash tables (the “old” way)
Hash of Reads (MAQ, ELAND, ZOOM, …)
Smaller but variable memory requirements (depends on
the amount of reads).
Hash of Reference (SOAP, MOSAIK, … )
Predictable memory requirements.
(2) Suffix arrays (the “new” way)
BWA, Bowtie, SOAP2, …
Indexing
• Genomes and reads are too large for direct
approaches like dynamic programming
• Indexing is required
Suffix tree
Suffix array
Seed hash tables
Many variants, incl. spaced seeds
• Choice of index is key to performance
Suffix Array
Find "ctat” in the reference
Is the process reversible ?
NGS Read Alignment
Burrows Wheeler Transformation (BWT)
• Invented by David Wheeler in 1983 (Bell Labs). Published in 1994.
“A Block Sorting Lossless Data Compression Algorithm”
Systems Research Center Technical Report No 124. Palo Alto, CA: Digital Equipment
Corporation, Burrows M, Wheeler DJ. 1994
• Originally developed for compressing large files (bzip2, etc.)
• Lossless, Fully Reversible
• Alignment Tools based on BWT: bowtie, BWA, SOAP2, etc.
• Approach:
– Align reads on the transformed reference genome, using an efficient index (FM index)
– Solve the simple problem first (align one character) and then build on that solution to solve a
slightly harder problem (two characters) etc.
• Results in great speed and efficiency gains (a few GigaByte of RAM for the entire H.
Genome). Other approaches require tens of GigaBytes of memory and are much
slower.
FM-index
Burrows-Wheeler Transform is the basis of the bzip2
file compression tool. B-W uses an FM-index (Ferragina
& Manzini) which allows efficient finding of substring
matches within compressed text – algorithm is
sub-linear with respect to time and storage space
required for a certain set of input data (reference
genome)
Reduced memory footprint, faster execution.
Burrows-Wheeler
• Store entire reference
genome.
• Align tag base by base from
the end.
• When tag is traversed, all
active locations are
reported.
• If no match is found, then
back up and try a
substitution.
Jianhua Ruan
The University of Texas at San Antonio
Burrows-Wheeler Transform (BWT)
BWT
acaacg$
$acaacg
aacg$ac
acaacg$
acg$aca
caacg$a
cg$acaa
g$acaac
gc$aaac
Burrows-Wheeler Matrix (BWM)
Burrows-Wheeler Matrix
$acaacg
aacg$ac
acaacg$
acg$aca
caacg$a
cg$acaa
g$acaac
Burrows-Wheeler Matrix
3
1
4
2
5
6
$acaacg
aacg$ac
acaacg$
acg$aca
caacg$a
cg$acaa
g$acaac
See the suffix array?
Key observation
a1c1a2a3c2g1$1
“last first (LF) mapping”
The i-th occurrence of character X in the
last column corresponds to
the same text character as the i-th
occurrence of X in the first column.
1$acaacg1
2aacg$ac1
1acaacg$1
3acg$aca2
1caacg$a1
2cg$acaa3
1g$acaac2
Exact match (another example)
BWT(agcagcagact) = tgcc$ggaaaac
gca
gca
Search for pattern: gca
gca
gca
$agcagcagact
$agcagcagact
$agcagcagact
$agcagcagact
act$agcagcag
act$agcagcag
act$agcagcag
act$agcagcag
agact$agcagc
agact$agcagc
agact$agcagc
agact$agcagc
agcagact$agc
agcagact$agc
agcagact$agc
agcagact$agc
agcagcagact$
agcagcagact$
agcagcagact$
agcagcagact$
cagact$agcag
cagact$agcag
cagact$agcag
cagact$agcag
cagcagact$ag
cagcagact$ag
cagcagact$ag
cagcagact$ag
ct$agcagcaga
ct$agcagcaga
ct$agcagcaga
ct$agcagcaga
gact$agcagca
gact$agcagca
gact$agcagca
gact$agcagca
gcagact$agca
gcagact$agca
gcagact$agca
gcagact$agca
gcagcagact$a
gcagcagact$a
gcagcagact$a
gcagcagact$a
t$agcagcagac
t$agcagcagac
t$agcagcagac
t$agcagcagac
Test with your own seq and pattern at: http://www.allisons.org/ll/AlgDS/Strings/BWT/
Exact Matching with FM Index
• To match Q in T using BWT(T), repeatedly apply
rule:
top = LF(top, qc); bot = LF(bot, qc)
– Where qc is the next character in Q (right-to-left) and
LF(i, qc) maps row i to the row whose first character
corresponds to i’s last character as if it were qc
FM Index is Small
• Entire FM Index on DNA reference consists of:
– BWT (same size as T)
– Checkpoints (~15% size of T)
– SA sample (~50% size of T)
Assuming 2-bit-per-base encoding and
no compression, as in Bowtie
Assuming a 16-byte checkpoint every
448 characters, as in Bowtie
Assuming Bowtie defaults for suffixarray sampling rate, etc
• Total: ~1.65x the size of T
~1.65x
>45x
>15x
>15x
FASTQ Format:
• The de-facto file format for sharing sequence read data
• Sequence and a per-base quality score
SAM (Sequence Alignment/Map) format:
• A unified format for storing read alignments to a reference
genome.
• Generally large files (a byte per bp)
• Very compact in size but computationally efficient to access.
BAM (Binary Alignment/Map) format:
• A Binary equivalent to SAM.
• Developed for fast processing and indexing
http://bioinformatics.oxfordjournals.org/cgi/reprint/btp352v1
Short-read mapping software
Software
Technique
Developer
License
Eland
Hashing reads
Illumnia
?
SOAP
Hashing refs
BGI
Academic
Maq
Hashing reads
Sanger (Li, Heng) GNUPL
Bowtie
BWT
Salzberg/UMD
BWA
BWT
Sanger (Li, Heng) GNUPL
SOAP2
BWT & hashing
BGI
GNUPL
Academic
http://www.oxfordjournals.org/our_journals/bioinformatics/nextgenerationsequencing.html
FASTQ Files
Sequence Id (Illumina)
36 bps read
@HWUSI-EAS610_0001:3:1:4:1405#0/1
GATAGTTCAATTCCAGAGATCAGAGAGAGGTGAGTG
+
B;30;<[email protected]/5@=?5?7?1>A2?0<6?<<80>79##
36 Quality scores
•
•
•
•
•
The de-facto file format for sharing DNA sequence read data
4 Lines per read
Sequence line and a per-base Phred quality score line per read
FASTQ Files are Text files
There is No file Header
An Introduction to Phred Quality Score
  10

Q Phred
 is the Error Probability: The
probability that a base call is wrong.
10
Q Phred  10  log 10 ( )
Q: Phred Quality Score
Q

Probability the base call
in wrong (confidence)
40
0.0001
0.01% (99.99%)
30
0.001
0.1% (99.9%)
20
0.01
1% (99%)
10
0.1
10% (90%)
• Phred Quality Score encoding in FASTQ/SAM files: ASCII Character = Q + 33
• FASTQ Files: Q represents Base Call Quality: Probability the base call is wrong.
• SAM Files: Q represents Mapping Quality: Probability the mapping position of the read is incorrect.
$perl –e ‘print chr(33);’
http://en.wikipedia.org/wiki/FASTQ_format
The SAM file
SAM data fields
Mapping Quality (MAPQ) in BWA
Mapping Quality is a function of Edit Distance and the Uniqueness of the alignment.
BWA Mapping Quality
0
A read aligns equally well to multiple positions
(hits). BWA picks randomly one of the positions
and assigns MAPQ=0
1 – 23
25
Only 1 Best hit (with no suboptimal hits) with
more than 2 mismatches.
Or
Only 1 Best hit, with 1 suboptimal hit.
37
Only 1 Best hit (no suboptimal hits), with up to 2
mismatches
(edit distance could be more than 2)
SAM/BAM format
@HD
@SQ
@RG
@RG
Header section
VN:1.0
SN:chr20 AS:HG18 LN:62435964
ID:L1 PU:SC_1_10 LB:SC_1 SM:NA12891
ID:L2 PU:SC_2_12 LB:SC_2 SM:NA12891
position of alignment
Alignment section
Query Name
V00-HWI-EAS132:3:38:959:2035#0
V00-HWI-EAS132:4:99:122:772#0
V00-HWI-EAS132:4:44:473:970#0
V00-HWI-EAS132:4:29:113:1934#0
Ref sequence
147
177
25
99
chr1
chr1
chr1
chr1
query sequence (same strand as ref)
28
32
40
44
255
255
255
255
36M
36M
36M
36M
=
=
*
=
79
9127
0
107
0
0
0
0
GATCTGATGGCAGAAAACCCCTCTCAGTCCGTCGTG
AAAGGATCTGATGGCAGAAAACCCCTCTCAGTCCGT
GTCGTGGTGAAGGATCTGATGGCAGAAAACACCTCT
GGGTTTTCTGCCATCAGATCCTTTACCACGACAGAC
query quality
aaX`[\`Y^Y^]ZX``\EV_BBBBBBBBBBBBBBBB
aaaaaa\OWaI_\WL\aa`Xa^]\ZUaa[XWT\^XR
__YaZ`W[aZNUZ[U[_TL[KVVX^QURUTDRVZBB
aaaQaa__``]\\_^``^a^`a`_^^^_XQ[ZS\XX
NM:i:1
NM:i:1
NM:i:2
NM:i:1
Post-processing: Tools and programming APIs for
parsing and manipulating alignments:
Samtools: http://samtools.sourceforge.net/
Convert SAM to BAM and vice versa
Sort and Index BAM files
Merge multiple BAM files
Show alignments in text viewer
Remove Duplicates from PCR amplification step
Picard Tools: (Java-based)
http://picard.sourceforge.net/index.shtml
BAM is Indexed & Binary Compressed SAM
• BAM is indexed by genome location.
• Software toolkit allows other software to
extract sequence data from BAM at specific
genome – do not need to store entire data file
in memory during operation of program!
SAMTools/Picard
http://samtools.sourceforge.net/
• SAMTools is a simple toolkit to transform SAM
to BAM, merge, sort, index
• Can also calculate statistics (mean quality,
depth of coverage, etc.)
• filter duplicate reads
• create multiple alignments of all reads over a
genomic interval
• call variants
SAMTools commands
• Download and install SAMTools
http://sourceforge.net/projects/samtools/files/samtools/
Make the BAM file:
samtools view –bt ref.fa –o data.bam data.sam
Sort the BAM file:
samtools sort data.bam data.sorted.bam
Index the BAM file:
samtools index data.sorted.bam data.st_index.bam
http://www.htslib.org/doc/sam.html
http://samtools.sourceforge.net/samtools-c.shtml
Manual Reference Pages
- samtools (1)
DESCRIPTION
• Samtools is a set of utilities that manipulate alignments in the BAM
format. It imports from and exports to the SAM (Sequence Alignment/Map)
format, does sorting, merging and indexing, and allows to retrieve
reads in any regions swiftly.
SYNOPSIS
samtools
samtools
samtools
samtools
samtools
samtools
samtools
samtools
samtools
samtools
bcftools
bcftools
bcftools
view -bt ref_list.txt -o aln.bam aln.sam.gz
sort aln.bam aln.sorted
index aln.sorted.bam
idxstats aln.sorted.bam
view aln.sorted.bam chr2:20,100,000-20,200,000
merge out.bam in1.bam in2.bam in3.bam
faidx ref.fasta
pileup -vcf ref.fasta aln.sorted.bam
mpileup -C50 -gf ref.fasta -r chr3:1,000-2,000 in1.bam in2.bam
tview aln.sorted.bam ref.fasta
index in.bcf
view in.bcf chr2:100-200 > out.vcf
view -vc in.bcf > out.vcf 2> out.afs
Heng Li from the Sanger Institute wrote the C version of samtools.
Picard (Java) Toolkit
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
AddCommentsToBam
AddOrReplaceReadGroups
BamToBfq
BamIndexStats
BedToIntervalList
BuildBamIndex
CalculateHsMetrics
CleanSam
CollectAlignmentSummaryMetrics
CollectBaseDistributionByCycle
CollectGcBiasMetrics
CollectInsertSizeMetrics
CollectMultipleMetrics
CollectTargetedPcrMetrics
CollectRnaSeqMetrics
CollectWgsMetrics
CompareSAMs
CreateSequenceDictionary
DownsampleSam
ExtractIlluminaBarcodes
EstimateLibraryComplexity
FastqToSam
FifoBuffer
FilterSamReads
FilterVcf
FixMateInformation
GatherBamFiles
GenotypeConcordance
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
IlluminaBasecallsToFastq
IlluminaBasecallsToSam
CheckIlluminaDirectory
IntervalListTools
MakeSitesOnlyVcf
MarkDuplicates
MarkDuplicatesWithMateCigar
MeanQualityByCycle
MergeBamAlignment
MergeSamFiles
MergeVcfs
NormalizeFasta
ExtractSequences
QualityScoreDistribution
ReorderSam
ReplaceSamHeader
RevertSam
RevertOriginalBaseQualitiesAndAddMateCigar
SamFormatConverter
SamToFastq
SortSam
SortVcf
UpdateVcfSequenceDictionary
VcfFormatConverter
MarkIlluminaAdapters
SplitVcfs
ValidateSamFile
ViewSam
VcfToIntervalList
Visualization
• BAM file format also simplifies visualization of NGS
data.
• Must make a .BAI index for each BAM file using
samtools index command.
• BAM and .BAI files must be located on your own
computer
• Index allows viewer to quickly find and load only read
data for a specific genomic interval.
• Integrative Genomics Viewer (IGV) from Broad Institute
is our current favorite.
• Use the Java Webstart to download and run IGV (use
the 1.2 GB version):
http://www.broadinstitute.org/igv/
Summary
•
•
•
•
•
•
•
•
Challenge of NGS read allignment
Suffix Array
FM index
BWA algorithm
SAM/BAM format
CIGAR string for variants
Software to work directly with SAM/BAM files
Sequence viewers
Next Lecture: ChIP-seq

similar documents