Slide 1

Slide 1 text

Fast differential expression analysis annotation-agnostic across groups with biological replicates Leonardo Collado-Torres tweet: @fellgernon blog: tinyurl.com/FellBit

Slide 2

Slide 2 text

@fellgernon #biostatJC2013 Field overview Ultimate Goal What is the biological (genomic) cause, if any, of X disease? Currently What are the most likely genomic difference(s) between two+ groups?

Slide 3

Slide 3 text

@fellgernon #biostatJC2013 Tools • Molecular biology: reverse transcriptase • High-throughput sequencing • $$ and – > Large number of biological replicates • Computers • Biostatistics Image: http://bit.ly/15MVhSU

Slide 4

Slide 4 text

@fellgernon #biostatJC2013 High-throughput sequencing: DNA Mardis. Ann Rev of Genomics and Human Genetics 2008

Slide 5

Slide 5 text

@fellgernon #biostatJC2013 High-throughput sequencing Mardis. Ann Rev of Genomics and Human Genetics 2008

Slide 6

Slide 6 text

@fellgernon #biostatJC2013 High-throughput sequencing Mardis. Ann Rev of Genomics and Human Genetics 2008

Slide 7

Slide 7 text

@fellgernon #biostatJC2013 Not as pretty in real life Whiteford et al. Bioinformatics 2009

Slide 8

Slide 8 text

@fellgernon #biostatJC2013 Many things can go wrong, like phasing Whiteford et al. Bioinformatics 2009

Slide 9

Slide 9 text

@fellgernon #biostatJC2013 So why do we use HTS? http://cdn.memegenerator.net/ instances/400x/33773502.jpg

Slide 10

Slide 10 text

@fellgernon #biostatJC2013 Transcription -> Translation Alberts et al. Molecular Biology of the Cell, 4th edition 2002 Fig 6-21

Slide 11

Slide 11 text

@fellgernon #biostatJC2013 Mapping reads from mRNAs Trapnell et al. Nat. Biotech 2009

Slide 12

Slide 12 text

@fellgernon #biostatJC2013 Mapping result: our initial data n samples à 3 billion nt Adapted from @leekgroup

Slide 13

Slide 13 text

@fellgernon #biostatJC2013 Split by chromosome and filter n samples à ~760 million nt Rows with at least 1 sample with coverage > 5

Slide 14

Slide 14 text

@fellgernon #biostatJC2013 Test for base-level DE Adapted from @andrewejaffe

Slide 15

Slide 15 text

@fellgernon #biostatJC2013 F-statistic at each base-pair • Null model • Alternative Model

Slide 16

Slide 16 text

@fellgernon #biostatJC2013 Threshold on F-statistics F-statistic corresponding to p-value < 10-8 (F5,30 ) Adapted from @andrewejaffe

Slide 17

Slide 17 text

@fellgernon #biostatJC2013 Q-values: qvalue::qvalue Permute model matrices and find null regions for all chromosomes.

Slide 18

Slide 18 text

@fellgernon #biostatJC2013 How can we make it fast? • Avoid Input/Output as much as possible • Work by chromosome • Reduce memory – Run Length Encoding (IRanges::Rle) 0000111111222 = (0, 1, 2) (4, 6, 3) • Use multiple cores (parallel::mclapply) – Split data to use cores efficiently • Calculate F-stats using Rcpp (Has + and -)

Slide 19

Slide 19 text

@fellgernon #biostatJC2013 Finding candidate DERs: example dataRegions 450 500 550 600 segs 450 500 550 600 pieces 450 500 550 600 ders 450 500 550 600 450 500 550 600 0.0 1.0 2.0 Index f

Slide 20

Slide 20 text

@fellgernon #biostatJC2013 Remember to map back to genomic coords dataRegions 450 500 550 600

Slide 21

Slide 21 text

@fellgernon #biostatJC2013 Calculate F-stats and threshold 450 500 550 600 0.0 0.5 1.0 1.5 2.0 2.5 Index f

Slide 22

Slide 22 text

@fellgernon #biostatJC2013 Thresholding identifies segments segs 450 500 550 600

Slide 23

Slide 23 text

@fellgernon #biostatJC2013 Disjoin data regions & segments pieces 450 500 550 600

Slide 24

Slide 24 text

@fellgernon #biostatJC2013 Identify pieces that overlap segments ders 450 500 550 600

Slide 25

Slide 25 text

@fellgernon #biostatJC2013 Example: re-cap dataRegions 450 500 550 600 segs 450 500 550 600 pieces 450 500 550 600 ders 450 500 550 600 450 500 550 600 0.0 1.0 2.0 Index f

Slide 26

Slide 26 text

@fellgernon #biostatJC2013 Example: result Cluster for region with name COL6A1 and q-value 0.8256 chr21 chr21 Coverage 1 2 group CEU YRI Mean coverage 0.125 0.500 group CEU YRI Regions significantQval TRUE FALSE tx_name (gene_id) tx_name(gene_id) 47411000 47411200 47411400 47411600

Slide 27

Slide 27 text

@fellgernon #biostatJC2013 Public datasets • derfinderExample: – Blood CEU vs YRI non-related individuals • derHippo: – Brain hippocampus from cocaine addicts, alcohol addicts, and controls • derSnyder: – Michael Snyder time course (~1 year): 2 x diseases, recovery & healthy periods • derStem: – 5 stem cell types, 2 replicates per group

Slide 28

Slide 28 text

@fellgernon #biostatJC2013 P-values: derfinderExample

Slide 29

Slide 29 text

@fellgernon #biostatJC2013 P-values: derStem

Slide 30

Slide 30 text

@fellgernon #biostatJC2013 MA-style plots: derfinderExample

Slide 31

Slide 31 text

@fellgernon #biostatJC2013 derHippo: CO vs CT

Slide 32

Slide 32 text

@fellgernon #biostatJC2013 derHippo: ETOH vs CT

Slide 33

Slide 33 text

@fellgernon #biostatJC2013 derHippo Genome overview of regions with q-value < 0.1

Slide 34

Slide 34 text

@fellgernon #biostatJC2013 derHippo: cluster #2

Slide 35

Slide 35 text

@fellgernon #biostatJC2013 derfinderExample: cluster #6

Slide 36

Slide 36 text

@fellgernon #biostatJC2013 derSnyder: cluster #3

Slide 37

Slide 37 text

@fellgernon #biostatJC2013 Coverage adjustment? ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 6.0e+07 8.0e+07 1.0e+08 1.2e+08 1.4e+08 1.6e+08 1.8e+08 6.0e+08 8.0e+08 1.0e+09 1.2e+09 1.4e+09 1.6e+09 chr 1 total Cov vs Cov < quantile 0.9 Coverage at bases < quantile 0.9 Total coverage 0to186 186to294 294to322 322to400 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 1.5e+07 2.0e+07 2.5e+07 3.0e+07 3.5e+07 4.0e+07 6.0e+08 8.0e+08 1.0e+09 1.2e+09 1.4e+09 1.6e+09 chr 1 total Cov vs # bases with data Number bases with Cov > 0 Total coverage 0to186 186to294 294to322 322to400 derSnyder Similar to metagenomeSeq::cumNorm

Slide 38

Slide 38 text

@fellgernon #biostatJC2013 derfinder R package https://github.com/lcolladotor/derfinder https://github.com/lcolladotor/derfinderReport https://github.com/lcolladotor/derfinderExample

Slide 39

Slide 39 text

@fellgernon #biostatJC2013 Acknowledgements Leek Group Jeffrey Leek Alyssa Frazee Hopkins Sarven Sabunciyan Ben Langmead Lieber Institute (LIBD) Andrew Jaffe Harvard Rafa Irizarry Funding NIH (Aug 2012- July 2013) LIBD (Aug 2013 - now) CONACyT México