Transcript Yan Guo_08.05.2014

Vanderbilt Center for Quantitative
Sciences
Summer Institute
Sequencing Analysis
Yan Guo
What is Sequencing?
• Sequencing is the process of determining the precise order of
nucleotides.
• Non high throughput sequencing: Sanger Sequencing: The
basic chain termination method, developed by Frederick
Sanger in 1974. Generates all possible single-stranded DNA
molecules complementary to a given template, and beginning
at a common 5' base.
The Pros and Cons of Sanger
Sequencing
• Pros: Highly accurate
• targetable
• Cons:
• Cost $15 per /1000 base pairs, to sequencing the whole
genome will cost roughly: 30bil/1000x$15=$15m
• Low detection rate of alternative allele
Current Generation Sequencing
Illumina
ABI Solid
454 Life
Science
Price
Low
medium
High
Read Length
50-100
50-100
400-1000
Read Depth
High
High
Low
Difficulty
Easy
High
Easy
Sequencing Type By Source
• RNA: mRNA, Small RNA, Total RNA
• DNA: Whole Genome or targeted (Exome,
mitochondrial, genes of interest, etc)
Sequencing Data
• Raw Image data is more than 2TB per sample
• Raw data is about 5-15GB per single end sample or 10-30GB per pair end
sample for RNAseq or Exome Sequencing. Whole genome data can easily
exceed 200GB per sample.
• In general 5x raw data size is needed to finish processing
• Raw data is usually in FASTQ format, the base quality is in Phred scale
• Older Illumina pipeline uses Phred 64 scale, newer CASAVA 1.8 pipeline
uses Sanger scale.
Single end vs Paired end
• Paired end data has double
amount of data than single
end.
• Paired end is more
expensive than single end.
• Paired end data is easier to
do quality control (insert
size, removing duplicate)
• Paired end data provides
more opportunities to
detect structural variance.
What can you obtain from DNAseq
•
•
•
•
•
SNPs (require only normal or tumor)
Somatic Mutations (require tumor and normal pair)
Copy Number Variation (work best with whole genome sequencing)
Small Structural Variance: Insertion, deletion
Large Structure Variance: (Translocation, Inversion)
What can you obtain from RNAseq
•
•
•
•
Gene Expression
SNP (only for expressed genes)
Novel Splicing Variants
Genes Fusion
RNAseq has been
used primarily as a
replacement of
microarray
How does RNAseq compare to
Microarray?
• Since 2008, people has been saying that
RNAseq will replace microarray for gene
expression profiling.
• VANTAGE stopped offering microarray service
earlier this year.
• Wang, Z., M. Gerstein, and M. Snyder, RNA-Seq: a revolutionary tool for
transcriptomics. Nat Rev Genet, 2009. 10(1): p. 57-63.
• 2. Shendure, J., The beginning of the end for microarrays? Nat Methods,
2008. 5(7): p. 585-7.
Data Distribution
Guo, Y., et al., Large Scale Comparison of Gene Expression Levels by Microarrays and RNAseq Using TCGA Data. PLoS One, 2013. 8(8): p.
e71462.
Result Consistency
Guo, Y., et al., Large Scale Comparison of Gene Expression Levels by Microarrays and RNAseq Using TCGA Data. PLoS One, 2013.
8(8): p. e71462.
RNAseq vs Microarray - advantages
RNAseq
Miroarray
Result Type
Rich, not limited to expression
Limited to expression only
Expression
Can quantify expression on exon
and gene level
Can quantify expression on exon or
gene level
Novel Discovery
Can be used for novel discovery
Can only detect what is on the chip
Analysis
Difficult
Easy
Interpretation
Difficult
Easy
Price for assay
Price has become comparable to
microarray, however the analysis
Price is stable
hardware and analysis time may
increase the final cost
Processing RNA
Raw data
•
•
•
•
@HWI-ST508:203:D078GACXX:8:1101:1296:1011 1:N:0:ATCACG
NTGGAGTCCTAGGCACAGCTCTAAGCCTCCTTATTCGAGCCGAGCTGGGCC
+
#4=DDDDDDDDDDE<DAEEEIDFEIEIEIEIIIIIIDEDDDDA@DDDDII@
@EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG
EAS139
the unique instrument name
136
the run id
FC706VJ
the flowcell id
2
flowcell lane
2104
tile number within the flowcell lane
15343
'x'-coordinate of the cluster within the
tile
197393
'y'-coordinate of the cluster within the
tile
1
the member of a pair, 1 or 2 (paired-end
or mate-pair reads only)
Y
Y if the read fails filter (read is bad), N
otherwise
18
0 when none of the control bits are on,
otherwise it is an even number
ATCACG
index sequence
Phred Score
Phred Quality Score
Probability of incorrect
base call
Base call accuracy
10
1 in 10
90 %
20
1 in 100
99 %
30
1 in 1000
99.9 %
40
1 in 10000
99.99 %
50
1 in 100000
99.999 %
Quality Control
•
Quality control should be
conducted at multiple steps
during sequencing data
processing
– Raw data
– Alignment
– Results (Expression for RNA,
and SNP/mutation for DNA)
Guo, Y., et al., Three-stage quality control strategies for DNA re-sequencing data. Brief Bioinform, 2013.
Raw Data QC - Tools
• FAST QC
http://www.bioinformatics.babraham.ac.uk/p
rojects/fastqc/
• FASTX-Toolkit
http://hannonlab.cshl.edu/fastx_toolkit/
• QC3 https://github.com/slzhao/QC3
• NGS QC Toolkit
http://59.163.192.90:8080/ngsqctoolkit/
Raw Data QC - What to Look For
Alignment QC - Tools
• QC3 https://github.com/slzhao/QC3
• Qqplot
http://genome.sph.umich.edu/wiki/QPLOT
• SAMStat http://samstat.sourceforge.net/
Alignment QC - What to Look For
Expression QC - Tools
• MultiRankSeq
https://github.com/slzhao/MultiRankSeq
Clustering Algorithms
• Start with a collection of n objects each
represented by a p–dimensional feature
vector xi , i=1, …n.
• The goal is to divide these n objects into k
clusters so that objects within a clusters are
more “similar” than objects between
clusters. k is usually unknown.
• Popular methods: hierarchical, k-means,
SOM, mixture models, etc.
Distance Calculation in Sequencing
Smith-Waterman algorithm
Sequence 1 = ACACACTA
Sequence 2 = AGCACACA
w(gap) = 0
w(match) = +2
w(a, − ) = w( − ,b) = w(mismatch) = − 1
Distance Calculation in Microarray
• Pearson Correlation
 x1 

x     and
 x N 
Two profiles (vectors)

 
C pearson ( x , y ) 
 y1 

y    
 y N 
N
( xi  mx )( yi  m y )
i 1
[i 1 ( xi  mx ) ][ i 1 ( yi  m y ) 2 ]
N

x

y
2
N
+1  Pearson Correlation  – 1
mx 
1
N

my 
1
N

N
x
n 1 n
N
n 1
yn

x

y
Similarity Measurements
• Euclidean Distance
 
d ( x, y) 
2
(
x

y
)
n 1 n n
N
 x1 
 y1 


x     y    
 
 x N 
 y N 
Linkage
• Single Linkage: D(X, Y) = min(d(x, y)), x ϵ X,
yϵY
• Complete Linkage: D(X, Y) = max(d(x, y)),
x ϵ X, y ϵ Y
• Average Linkage:
Experssion QC - What to Look For
Batch Effect
Correction of Batch Effect
Guo, Y., et al., Statistical strategies for microRNAseq batch effect reduction. Translational Cancer Research, 2014. 3(3): p. 260-265.
Normalization of RNAseq
Reads Per
Kilo base per Million reads (RPKM)
RNAseq Data Alignment
• TopHat2
http://ccb.jhu.edu/software/tophat/index.sht
ml
• MapSplice
http://www.netlab.uky.edu/p/bioinfo/MapSpl
ice
Gene Quantification
• CufflInks for RPKM
http://cufflinks.cbcb.umd.edu/
• HTSeq for read count http://wwwhuber.embl.de/users/anders/HTSeq/doc/over
view.html
Data
Gene Symbol
DDR1
RFC2
HSPA6
PAX8
GUCA1A
UBE1L
THRA
PTPN21
CCL5
CYP2E1
EPHB3
ESRRA
CYP2A6
SCARB1
TTLL12
C2orf59
WFDC2
MAPK1
MAPK1
ADAM32
1
9.376298
7.950475
5.584798
6.355186
2.961001
7.437969
6.606546
7.392678
2.710744
3.871231
4.289411
7.151026
4.568492
6.134823
9.346916
4.42666
4.706794
4.777312
7.875045
4.629726
2
8.961996
7.795976
5.12491
6.245788
3.226968
7.422707
6.687768
6.772702
2.479818
4.085553
3.771091
7.219117
4.33565
6.440855
8.955574
5.219388
4.974295
4.797072
7.902457
5.27395
3
9.271935
7.124782
5.77907
6.388794
3.092915
8.298944
6.910623
6.834253
2.51898
5.031865
3.798425
6.900173
4.5123
5.739945
8.868433
4.799542
5.149892
4.249238
7.572943
4.351249
4
8.968211
8.156603
5.849914
6.737545
3.187618
6.124551
7.166293
6.840313
2.61285
5.053069
3.893421
7.841436
4.672211
6.269867
9.825905
5.204245
4.417064
4.252584
8.10576
5.249061
5
8.663588
7.821047
5.593596
6.662428
3.067353
6.263097
6.711748
6.813115
2.885117
5.080394
4.01667
7.254173
4.587597
5.534482
9.387397
4.846079
4.273504
3.687591
7.793828
5.204216
6
9.214028
6.613421
5.042853
6.279758
3.159364
7.548323
6.632955
6.68312
2.668616
5.557095
4.200385
7.119073
4.561608
5.281546
9.1008
3.934838
4.638822
4.412024
7.635768
5.412291
Example of Quantile Normalization
Red = G1; Green = G2; Blue = G3; Yellow = G4; Black = G5
Original
Sort S1
Sort S2
Sort S3
S1 S2
S3
S1
S2
S3
G1
2
4
4
2
3
4
G2
5
4
14
3
4
G3
4
6
8
3
G4
3
5
8
G5
3
3
9
Sorted
S1 S2
S3
G1
2
3
4
8
G2
3
4
8
4
8
G3
3
4
8
4
5
9
G4
4
5
9
5
6
14
G5
5
6
14
Take Average for Each Row
Sorted
S1 S2
S3
S1 S2
S3
S1 S2
S3
S1 S2
S3
S1 S2
S3
2
3
4
3
3
3
3
3
3
3
3
3
3
3
3
4
8
5
5
5
5
5
5
5
5
5
3
4
8
5
5
5
5
5
5
4
5
9
6
6
6
5
6
14
3
Averaged
S1 S2
S3
3
3
3
5
5
5
5
5
5
6
6
6
8
8
8
Reorder
Red = G1; Green = G2; Blue = G3; Yellow = G4; Black = G5
Averaged
S1 S2
S3
S1 S2
S3
S1 S2
S3
S1 S2
S3
S1 S2
S3
3
3
3
3
3
3
5
3
3
5
3
3
5
3
5
5
5
8
5
8
8
5
8
8
5
8
5
5
5
6
8
5
6
8
5
6
6
6
5
6
5
8
8
8
5
S1 S2
S3
3
5
3
8
5
8
6
8
5
5
6
5
5
3
6
Differential Expression Analysis
• Cuffdiff from Cufflinks package Trapnell, C., et al., Differential gene and transcript expression analysis of RNA-seq
experiments with TopHat and Cufflinks. Nat Protoc, 2012. 7(3): p. 562-78.
• DESeq
http://bioconductor.org/packages/release/bioc/html/DESeq.html
• EdgeR
• http://www.bioconductor.org/packages/release/bioc/html/edgeR.h
tml
• NBPSeq http://cran.rproject.org/web/packages/NBPSeq/index.html
• TSPM http://omictools.com/sequencing/rna-seq/normalizationde/tspm-r-s2496.html
• baySeq
http://www.bioconductor.org/packages/release/bioc/html/baySeq.
html
Which Method Is the Best?
Guo, Y., et al., Evaluation of read count based RNAseq analysis methods. BMC Genomics, 2013. 14 Suppl 8: p. S2.
Consistency
Consistency
Inconsistency
Method
DESeq
edgeR
baySeq
Cuffdiff
Adj pvalue
0.278
0.047
0.907
<0.001
Disease1
Read count (IGHG2) 391
Total Read count
49870084
Adjusted Read Count 78
Log2FC
3.00
2.92
NA
5.83
Rank
2572
712
24962
13
Disease2
2038
65550902
311
Disease3
338
71454121
47
Control1
634
35641084
178
Control2
10282
44863975
2292
Control3
1764
49052840
360
Combined Approach
1log2FoldChan pValue(DESe
log2FoldChan pValue(edge
log2FoldChan Likelihood(ba AdjLikelihood
ge(DESeq2) q2)
pAdj(DESeq2) ge(edgeR)
R)
pAdj(edgeR) ge(raw)
ySeq)
(baySeq)
rank(DESeq) rank(edgeR) rank(baySeq) rankMethod1
ENSMUSG000000
90862_Rps13
-5.86335
7.02E-210
1.07E-205
-6.19904
1.80E-109
4.01E-105
-6.02447
4.21E-07
1.81E-07
1
1
4
6
ENSMUSG000000
58546_Rpl23a
-3.67515
3.27E-140
2.49E-136
-3.75807
2.14E-58
2.38E-54
-3.57503
2.53E-05
5.21E-06
2
2
5
9
ENSMUSG000000
91957_Gm8841
-4.68658
1.91E-72
7.27E-69
-5.33723
1.60E-50
8.90E-47
-5.14873
4.71E-05
1.70E-05
4
4
8
16
ENSMUSG000000
62683_Atp5g2
-4.62274
4.86E-69
1.48E-65
-5.27352
4.06E-50
1.81E-46
-5.06178
7.54E-05
2.83E-05
5
5
10
20
ENSMUSG000000
82697_Gm12913
-3.94956
8.13E-80
4.13E-76
-4.21738
4.01E-52
2.98E-48
-4.10532
0.000296
9.17E-05
3
3
14
20
ENSMUSG000000
60128_Gm10075
-4.59774
7.59E-64
1.65E-60
-5.34137
2.30E-42
6.41E-39
-5.13809
2.68E-05
8.81E-06
7
8
6
21
ENSMUSG000000
58558_Rpl5
-3.06317
8.73E-68
2.22E-64
-3.1805
1.29E-42
4.12E-39
-2.99193
0.000313
0.000119
6
7
16
29
ENSMUSG000000
63316_Rpl27
-3.26559
2.86E-62
5.45E-59
-3.4434
2.48E-43
9.20E-40
-3.27313
0.000409
0.000151
8
6
18
32
ENSMUSG000000
73702_Rpl31
-2.71836
4.03E-41
4.72E-38
-2.8753
7.47E-33
1.67E-29
-2.70413
0.000453
0.000167
13
10
19
42
ENSMUSG000000
85279_Gm15965
-4.21114
1.41E-33
1.34E-30
-6.49343
3.47E-24
4.29E-21
-6.53916
7.18E-05
2.31E-05
16
18
9
43
ENSMUSG000000
78686_Mup9
-3.77538
8.13E-32
6.19E-29
-4.7123
5.65E-22
5.47E-19
-4.58514
1.71E-13
1.71E-13
20
23
1
44
ENSMUSG000000
93337_Mir5109
2.810771
9.06E-36
9.20E-33
3.112215
3.04E-32
6.16E-29
3.311894
0.000751
0.000244
15
11
22
48
ENSMUSG000000
49517_Rps23
-2.38227
3.15E-44
4.37E-41
-2.45472
9.73E-29
1.67E-25
-2.25997
0.001089
0.000337
11
13
25
49
Guo, Y., et al., MultiRankSeq: Multiperspective Approach for RNAseq Differential Expression Analysis and Quality Control. BioMed Research International,
2014. 2014: p. 8.
Presentation Using Heatmap and
Cluster
Zhao, S., et al., Advanced Heat Map and Clustering Analysis Using Heatmap3. BioMed Research International, 2014. 2014: p. 6.
Difference Between Heatmaps
Questions We Can Answer with Cluster
• Microarray data quality checking
– Does replicates cluster together?
– Does similar conditions, time points, tissue types cluster together?
Presentation Using Volcano Plot
Presentation Using Circos Plot
Test Your Hypothesis Without
Performing Any Analysis
• GEO http://www.ncbi.nlm.nih.gov/geo/
Test Your Hypothesis Without
Performing Any Analysis
Functional Analysis
Samples Space n
F
M
Suppose in a study, we are
trying to find out if the
proportion of smoking
individual is significantly
different between men and
women.
Smoking
c
a
b
d
Fisher’s Exact Test
Male
Female
Total
a
b
a+b
Nonsmoking c
d
c+d
Total
b+d
a+b+c+d=n
Smoking
a+c
H0 : The proportion of smoking in male == the proportion of smoking in female
H1 : The proportion of smoking in male != the proportion of smoking in female
http://www.graphpad.com/quickcalcs/contingency1.cfm
Fisher’s Exact Test – in Functional
Analysis
All Genes
d
Breast
Cancer
Genes
b
a
Winner
Genes
c
Breast
Cancer
Genes
Non
Brest
Cancer
Genes
Winner
Genes
Non
Winner
Genes
a
b
c
d
Analogy
• There are 18000 Balls: 200 + 17800 in a
box.
• Blindfolded, you randomly draw 100 balls.
• What is the probability that you draw less
than 50
WebGestalt
• http://bioinfo.vanderbilt.edu/webgestalt/
Gene Set Enrichment Analysis
• KS test based analysis (Ref)
• GSEA does not need a winner list first
• http://www.broadinstitute.org/gsea/index.jsp
SNV and Indel
• Difficulty due to high false positive rate
• RNAMapper (Miller, et al. Genome Research,
2013)
• SNVQ (Duitama, et al. (BMC Genomics, 2013)
• FX (Hong, et al. Bioinformatics, 2012)
• OSA (Hu, et al. Binformatics, 2012)
Microsatellite instability
Examples:
• Yoon, et al. Genome Research 2013
• Zheng, et al. BMC Genomics, 2013
RNA Editing and Allele-specific
expression
RNA editing tools and database
• DARNED, REDidb, dbRES, RADAR
Allele-specific expression
• asSeq (Sun, et al. Biometrics, 2012)
• AlleleSeq (Rozowsky, et al. Molecular Systems
Biology, 2011)
Exogenous RNA
• Virus (Same as DNA)
• Food RNA (you are what you eat)
Wang, et al. PLOS ONE, 2012
nonCoding RNA