SLIDE 1
RNA-seq read mapping Pr Engstrm SciLifeLab - - PowerPoint PPT Presentation
RNA-seq read mapping Pr Engstrm SciLifeLab - - PowerPoint PPT Presentation
RNA-seq read mapping Pr Engstrm SciLifeLab RNA-seq workshop April 2016 Input: sequence reads (FASTQ format) @D2BJQQN1:142:C03FNACXX:5:1101:12935:2174
SLIDE 2
SLIDE 3
Goal: ¡reads ¡mapped ¡to ¡genome ¡
Scale chr2: Common SNPs(144) RepeatMasker 2 kb hg19 136,872,000 136,873,000 136,874,000 136,875,000 136,876,000 Gm12878 RNA-seq reads Gm12878 RNA-seq coverage, minus strand RefSeq Genes 100 vertebrates Basewise Conservation by PhyloP Simple Nucleotide Polymorphisms (dbSNP 144) Found in >= 1% of Samples Repeating Elements by RepeatMasker CXCR4 CXCR4 Gm12878_minus 40 _ 0 _ 100 Vert. Cons 4.88 _
- 4.5 _
0 -
SLIDE 4
Spliced ¡alignment ¡
- k
- Garber ¡et ¡al. ¡Nature ¡Methods ¡2011 ¡
SLIDE 5
Introns ¡can ¡be ¡very ¡large! ¡
10 100 1,000 10,000 100,000 1M 10M 0.0 0.2 0.4 0.6 0.8 1.0
Human introns (Ensembl) Intron size (bp) Cumulative proportion
SLIDE 6
Limited ¡sequence ¡signals ¡at ¡splice ¡sites ¡
Iwata ¡and ¡Gotoh ¡BMC ¡Genomics ¡2011 ¡
- H. sapiens
- D. rerio
- D. melanogaster
- P. chrysosporium
- S. cerevisiae
- P. tricornutum
- A. thaliana
5’ss BP 3’ss
- 3
- 2
- 1
1 2 5 4 3
- 3
- 2
- 1
1 2 5 4 3
- 3
- 2
- 1
1 2 5 4 3
- 3
- 2
- 1
1 2 5 4 3
- 3
- 2
- 1
1 2 5 4 3
- 3
- 2
- 1
1 2 5 4 3
- 3
- 2
- 1
1 2 5 4 3
- 3
- 2
- 1
1 2 5 4 3
- 12
- 11
- 10
- 9
- 8
- 7
- 6
- 5
- 4
- 3
- 2
- 1
- 12
- 11
- 10
- 9
- 8
- 7
- 6
- 5
- 4
- 3
- 2
- 1
- 12
- 11
- 10
- 9
- 8
- 7
- 6
- 5
- 4
- 3
- 2
- 1
- 12
- 11
- 10
- 9
- 8
- 7
- 6
- 5
- 4
- 3
- 2
- 1
- 12
- 11
- 10
- 9
- 8
- 7
- 6
- 5
- 4
- 3
- 2
- 1
- 12
- 11
- 10
- 9
- 8
- 7
- 6
- 5
- 4
- 3
- 2
- 1
- 12
- 11
- 10
- 9
- 8
- 7
- 6
- 5
- 4
- 3
- 2
- 1
GT…AG ¡ 98.6% ¡ GC…AG ¡ 1.3% ¡ AT…AC ¡ 0.1% ¡
SLIDE 7
MulX-‑mapping ¡reads ¡and ¡pseudogenes ¡
FuncXonal ¡gene ¡ Processed ¡pseudogene ¡ Correct ¡read ¡alignment ¡ IdenXcal, ¡spliced ¡ Incorrect ¡read ¡alignment ¡ Mismatches, ¡not ¡spliced ¡ Note: ¡
- An ¡aligner ¡may ¡report ¡both ¡alignments ¡or ¡either ¡
- Some ¡search ¡strategies ¡and ¡scoring ¡schemes ¡give ¡preference ¡to ¡unspliced ¡alignments ¡
SLIDE 8
Current ¡RNA-‑seq ¡aligners ¡
TopHat2 ¡ Kim ¡et ¡al. ¡Genome ¡Biology ¡2013 ¡ HISAT ¡& ¡HISAT2 ¡ Kim ¡et ¡al. ¡Nature ¡Methods ¡2015 ¡ STAR ¡ Dobin ¡et ¡al. ¡Bioinforma8cs ¡2013 ¡ GSNAP ¡ Wu ¡and ¡Nacu ¡Bioinforma8cs ¡2010 ¡ OLego ¡ Wu ¡et ¡al. ¡Nucleic ¡Acids ¡Research ¡2013 ¡ HPG ¡aligner ¡ Medina ¡et ¡al. ¡DNA ¡Research ¡2016 ¡ MapSplice2 ¡ hap://www.netlab.uky.edu/p/bioinfo/MapSplice2 ¡
SLIDE 9
The ¡predecessor: ¡BLAT ¡
“In ¡the ¡process ¡of ¡assembling ¡and ¡annotaXng ¡the ¡human ¡genome, ¡I ¡was ¡faced ¡ with ¡two ¡very ¡large-‑scale ¡alignment ¡problems: ¡aligning ¡three ¡million ¡ESTs ¡and ¡ aligning ¡13 ¡million ¡mouse ¡whole-‑genome ¡random ¡reads ¡against ¡the ¡human ¡
- genome. ¡These ¡alignments ¡needed ¡to ¡be ¡done ¡in ¡less ¡than ¡two ¡weeks’ ¡Xme ¡
- n ¡a ¡moderate-‑sized ¡(90 ¡CPU) ¡Linux ¡cluster ¡in ¡order ¡to ¡have ¡Xme ¡to ¡process ¡
an ¡updated ¡genome ¡every ¡month ¡or ¡two. ¡To ¡achieve ¡this ¡I ¡developed ¡a ¡very-‑ high-‑speed ¡mRNA/DNA ¡and ¡translated ¡protein ¡alignment ¡algorithm. ¡“ ¡ ¡ (Kent ¡Genome ¡Research ¡2002) ¡
¡
SLIDE 10
InnovaXons ¡in ¡RNA-‑seq ¡alignment ¡soiware ¡
- Read ¡pair ¡alignment ¡
- Consider ¡base ¡call ¡quality ¡scores ¡
- SophisXcated ¡indexing ¡to ¡decrease ¡CPU ¡and ¡memory ¡usage ¡ ¡
- Map ¡to ¡geneXc ¡variants ¡
- Resolve ¡mulX-‑mappers ¡using ¡regional ¡read ¡coverage ¡
- Consider ¡juncXon ¡annotaXon ¡
- Two-‑step ¡approach ¡(juncXon ¡discovery ¡& ¡final ¡alignment) ¡
¡
SLIDE 11
Two-‑step ¡RNA-‑seq ¡read ¡mapping ¡
1st)run)of)HISAT)to)discover)splice)sites) 2nd)run)of)HISAT)to)align)reads)by)making)use)of)the)list)of)splice)sites)collected)above) ) ) ) ) ) )
e1# e3# e2#
) ) ) ) ) ) ) )
e1# e3# e2#
) )
mapped) unmapped) Read) Exon) Intron) Global)Search) Local)Search) Extension) Junc:on)extension)
Thi
Kim ¡et ¡al. ¡Nature ¡Methods ¡2015 ¡
SLIDE 12
Mapping ¡accuracy ¡
Accuracy ¡for ¡20 ¡million ¡simulated ¡human ¡100 ¡bp ¡reads ¡with ¡0.5% ¡mismatch ¡rate ¡ ¡
92.3 (93.8) 97.6 (99.2) 92.1 (93.5) 97.6 (99.3) 91.6 (91.6) 88.9 (90.5) 96.7 (99) 94.4 (97.4) 25 50 75 100 HISAT×1 OLego Percentage of reads GSNAP STAR STAR×2 HISAT HISAT×2 TopHat2 Correctly and uniquely mapped Correctly mapped (multimapped) Incorrectly mapped % % + %) ( Unmapped
Kim ¡et ¡al. ¡Nature ¡Methods ¡2015 ¡
SLIDE 13
Mapping ¡accuracy ¡for ¡reads ¡with ¡small ¡anchors ¡
Kim ¡et ¡al. ¡Nature ¡Methods ¡2015 ¡
90.4 (91.7) 97.2 (98.9) 88.2 (89.4) 97.2 (99) 88.5 (88.5) 51.1 (52.2) 94.4 (96.4) 93.5 (95.5) 9.2 (9.4) 94.4 (97.3) 7.6 (7.7) 95.4 (98.3) (0) (0) 79.8 (92.6) 77.8 (96) 25 50 75 100 25 50 75 100 Percentage of reads, 2M_8_15 Percentage of reads, 2M_1_7 HISAT×1 OLego GSNAP STAR STAR×2 HISAT HISAT×2 TopHat2 Correctly and uniquely mapped Correctly mapped (multimapped) Incorrectly mapped Unmapped % % + %) ( % % + %) ( 62.4% (M) 3.1% (gt_2M) 4.2% (2M_1_7) 5.1% (2M_8_15) 25.1% (2M_gt_15)
b
SLIDE 14
Mapping ¡accuracy ¡for ¡spliced ¡RNA-‑seq ¡reads ¡
TopHat2 ann TopHat2 TopHat1 ann TopHat1 STAR 2−pass ann STAR 2−pass STAR 1−pass ann STAR 1−pass SMALT ReadsMap PASS cons PASS PALMapper cons ann PALMapper cons PALMapper ann PALMapper MapSplice ann MapSplice GSTRUCT ann GSTRUCT GSNAP ann GSNAP GEM cons ann GEM cons GEM ann BAGET ann
Simulation 1 Simulation 2
Percent of simulated spliced reads
20 40 60 80 100
Percent of simulated spliced reads
5 10 15
Percent of simulated spliced reads
20 40 60 80 100
Percent of simulated spliced reads
5 10 15
Perfectly mapped Part correctly mapped Mapped, no base correct No base correcly mapped but intersecting correct location
High ¡accuracy ¡at ¡mapping ¡to ¡correct ¡locus: ¡GSNAP, ¡GSTRUCT, ¡MapSplice, ¡STAR ¡ ¡ High ¡rate ¡of ¡perfect ¡spliced ¡alignments: ¡ReadsMap, ¡TopHat2 ¡ann ¡ ¡ ¡ Engström ¡et ¡al. ¡Nature ¡Methods ¡2013 ¡
SLIDE 15
Major ¡differences ¡in ¡indel ¡frequencies ¡
TopHat2 ann TopHat2 TopHat1 ann TopHat1 STAR 2−pass ann STAR 2−pass STAR 1−pass ann STAR 1−pass SMALT ReadsMap PASS cons PASS PALMapper cons PALMapper MapSplice ann MapSplice GSTRUCT ann GSTRUCT GSNAP ann GSNAP GEM cons ann GEM cons GEM ann BAGET ann Insertions (%)
20 40 60 80 100
5.86 6.71 2.05 2.05 2.02 2.02 2.03 2.00 8.91 2.70 2.38 2.44 0.68 31.54 1.65 1.65 4.90 4.94 4.84 5.80 84.91 85.76 83.32 14.46 Deletions (%)
20 40 60 80 100
1 2 3 4 5–7 8+ Indel size (bases):
a
6.94 6.09 7.33 7.29 4.50 4.37 4.50 4.14 9.92 4.48 4.77 4.95 0.30 61.71 5.00 4.98 9.12 9.16 8.97 10.25 29.39 29.51 29.12 13.07 Indel ¡frequencies ¡are ¡tabulated ¡(number ¡of ¡indels ¡per ¡thousand ¡sequenced ¡reads). ¡Data ¡set: ¡K562 ¡(mean). ¡ ¡
Engström ¡et ¡al. ¡Nature ¡Methods ¡2013 ¡
SLIDE 16
Indel ¡accuracy ¡on ¡simulated ¡data ¡
- GSNAP ¡and ¡GSTRUCT ¡exhibit ¡high ¡sensiXvity ¡for ¡both ¡long ¡and ¡short ¡deleXons ¡
- TopHat2 ¡ann ¡is ¡most ¡sensiXve ¡for ¡long ¡inserXons ¡ ¡
¡
1 2 3 4 5 6 7 8+ 1 2 3 4 5 6 7 8+ 1 2 3 4 5 6 7 8+
TopHat2 ann TopHat2 TopHat1 ann TopHat1 STAR 2−pass ann STAR 2−pass STAR 1−pass ann STAR 1−pass SMALT ReadsMap PASS cons PASS PALMapper cons ann PALMapper cons PALMapper ann PALMapper MapSplice ann MapSplice GSTRUCT ann GSTRUCT GSNAP ann GSNAP GEM cons ann GEM cons GEM ann BAGET ann
1 2 3 4 5 6 7 8+
b
Precision Precision Recall Recall Insertions Deletions Indel size (bases)
0% 100% n.a.
Engström ¡et ¡al. ¡Nature ¡Methods ¡2013 ¡
SLIDE 17
Novel ¡juncXons ¡are ¡typically ¡supported ¡by ¡few ¡alignments ¡
50,000 100,000 150,000
Minimum number of supporting mappings Annotated junctions
1 2 3 4 5 6 7 8 9 10 100,000 200,000 300,000
Minimum number of supporting mappings Novel junctions
1 2 3 4 5 6 7 8 9 10
Annotated junctions K562 whole cell K562 whole cell
b
Engström ¡et ¡al. ¡Nature ¡Methods ¡2013 ¡
SLIDE 18
Improved ¡juncXon ¡accuracy ¡by ¡filtering ¡on ¡coverage ¡
5,000 10,000 15,000 20,000 25,000 100,000 110,000 120,000
True junctions
- Simulation 1
False junctions
- BAGET ann
GEM ann GEM cons GEM cons ann GSNAP GSNAP ann GSTRUCT GSTRUCT ann MapSplice MapSplice ann PALMapper
- PALM. ann
- PALM. cons
- PALM. cons ann
PASS PASS cons ReadsMap SMALT STAR 1p STAR 1p ann STAR 2p STAR 2p ann TopHat1 TopHat1 ann TopHat2 TopHat2 ann
Engström ¡et ¡al. ¡Nature ¡Methods ¡2013 ¡
SLIDE 19
Improved ¡juncXon ¡accuracy ¡by ¡filtering ¡on ¡coverage ¡
- BAGET ann
GEM ann GEM cons GEM cons ann GSNAP GSNAP ann GSTRUCT GSTRUCT ann MapSplice MapSplice ann PALMapper
- PALM. ann
- PALM. cons
- PALM. cons ann
PASS PASS cons ReadsMap SMALT STAR 1p STAR 1p ann STAR 2p STAR 2p ann TopHat1 TopHat1 ann TopHat2 TopHat2 ann
10,000 20,000 30,000 40,000 100,000 110,000 120,000
False junctions True junctions
- Simulation 2
Engström ¡et ¡al. ¡Nature ¡Methods ¡2013 ¡
SLIDE 20
Several ¡methods ¡show ¡over-‑confidence ¡in ¡annotaXon ¡
- BAGET ann
GEM ann GEM cons GEM cons ann GSNAP GSNAP ann GSTRUCT GSTRUCT ann MapSplice MapSplice ann PALMapper
- PALM. ann
- PALM. cons
- PALM. cons ann
PASS PASS cons ReadsMap SMALT STAR 1p STAR 1p ann STAR 2p STAR 2p ann TopHat1 TopHat1 ann TopHat2 TopHat2 ann
Annotated junctions
c d e
25,000 5,000 10,000 15,000 20,000 68,000 72,000 76,000 80,000
- 80,000
Simulation 1 False junctions True junctions
Engström ¡et ¡al. ¡Nature ¡Methods ¡2013 ¡
SLIDE 21
Top ¡performers ¡(RGASP) ¡
In ¡general, ¡GSNAP, ¡GSTRUCT, ¡MapSplice ¡and ¡STAR ¡compared ¡ favorably ¡to ¡the ¡other ¡methods, ¡but ¡also ¡displayed ¡certain ¡ weaknesses: ¡ ¡
- MapSplice ¡is ¡a ¡conservaXve ¡aligner, ¡both ¡with ¡respect ¡to ¡
mismatch ¡frequency, ¡indel ¡calls ¡and ¡exon ¡juncXon ¡calls. ¡ ¡
- The ¡largest ¡issue ¡with ¡GSNAP, ¡GSTRUCT ¡and ¡STAR ¡is ¡the ¡
presence ¡of ¡many ¡false ¡exon ¡juncXons ¡in ¡the ¡output. ¡
Engström ¡et ¡al. ¡Nature ¡Methods ¡2013 ¡
SLIDE 22
Compute ¡requirements ¡
- Program
Run time (min) Memory usage (GB) HISATx1 22.7 4.3 HISATx2 47.7 4.3 HISAT 26.7 4.3 STAR 25 28 STARx2 50.5 28 GSNAP 291.9 20.2 OLego 989.5 3.7 TopHat2 1,170 4.3
Run times and memory usage for HISAT and other spliced aligners to align 109 million 101-bp RNA-seq reads from a lung fibroblast data set. We used three CPU cores to run the programs on a Mac Pro with a 3.7 GHz Quad-Core Intel Xeon E5 processor and 64 GB of RAM.
Kim ¡et ¡al. ¡Nature ¡Methods ¡2015 ¡
SLIDE 23
RecommendaXons ¡
- Use ¡a ¡two-‑pass ¡workflow ¡
- STAR ¡and ¡GSNAP ¡generally ¡perform ¡well ¡
- HISAT ¡also ¡seems ¡to ¡do ¡well ¡
- HISAT ¡and ¡STAR ¡are ¡the ¡fastest ¡
- HISAT2 ¡has ¡not ¡been ¡evaluated ¡but ¡the ¡authors ¡recommend ¡it ¡over ¡HISAT ¡
- HISAT2 ¡and ¡GSNAP ¡can ¡use ¡SNP ¡data, ¡which ¡may ¡give ¡higher ¡sensiXvity ¡
- If ¡you ¡want ¡to ¡run ¡Cufflinks, ¡use ¡TopHat2 ¡or ¡HISAT2 ¡
- For ¡long ¡(PacBio) ¡reads, ¡STAR, ¡BLAT ¡or ¡GMAP ¡can ¡be ¡used ¡
SLIDE 24
Important ¡SAM ¡fields ¡
Command: samtools view –X file.bam Perfectly ¡and ¡uniquely ¡aligned ¡read ¡pair: HWI-ST1018:3:1305:21090:45397#0 pPR1 chr1 4426 255 101M = 4435 110 GT… C@… NH:i:1 HI:i:1 AS:i:200 nM:i:0 HWI-ST1018:3:1305:21090:45397#0 pPr2 chr1 4435 255 101M = 4426 -110 CG… 5<… NH:i:1 HI:i:1 AS:i:200 nM:i:0 Problema<c ¡read ¡pair: HWI-ST1018:3:2109:6170:66353#0 pPR2s chr1 5058 3 65M36S = 5058 95 CA… B@… NH:i:2 HI:i:2 AS:i:135 nM:i:9 HWI-ST1018:3:2109:6170:66353#0 pPr1s chr1 5058 3 7S73M1D21M = 5058 -95 CC… ##… NH:i:2 HI:i:2 AS:i:135 nM:i:9
SLIDE 25