THE RNA-SEQ ANALYSIS PIPELINE Alicia Oshlack Murdoch - - PowerPoint PPT Presentation

THE ¡RNA-­‑SEQ ¡ANALYSIS ¡PIPELINE ¡

Alicia ¡Oshlack ¡ Murdoch ¡Childrens ¡Research ¡Ins5tute ¡

Two ¡ways ¡to ¡look ¡at ¡sequencing ¡data ¡

Sequence ¡of (mapped) ¡read ¡

• genome ¡sequencing ¡ ¡
• variant ¡detec5on ¡
• Muta5on ¡detec5on ¡
• genomic ¡rearrangements ¡
• Bisulfite-­‑seq ¡(methyla5on) ¡
• RNA ¡edi5ng ¡etc. ¡

Posi5on ¡of ¡mapped ¡ read ¡

• RNA-­‑seq ¡
• ChIP-­‑seq ¡ ¡
• MeDIP-­‑seq ¡for ¡DNA ¡methyla5on ¡
• etc. ¡

2 ¡

Two ¡ways ¡to ¡look ¡at ¡RNA-­‑seq ¡data ¡

Sequence ¡of (mapped) ¡read ¡

• ¡Determining ¡genes/

transcripts ¡

Posi5on ¡of ¡mapped ¡ read ¡

• Expression ¡levels ¡
• Differen5al ¡expression ¡

3 ¡

Benefits ¡and ¡opportuni5es ¡of ¡RNA-­‑seq ¡

• All ¡transcripts ¡are ¡sequenced ¡not ¡just ¡ones ¡for ¡

which ¡probes ¡are ¡designed ¡(cf ¡microarrays) ¡

• Annota5on ¡of ¡new ¡exons, ¡transcribed ¡regions, ¡

genes ¡or ¡non-­‑coding ¡RNAs ¡

• Whole ¡transcriptome ¡sequencing ¡

– The ¡ability ¡to ¡look ¡at ¡alterna5ve ¡splicing ¡ – Allele ¡specific ¡expression ¡ – RNA ¡edi5ng ¡

This ¡talk ¡

• Analysis ¡of ¡RNA-­‑seq ¡data ¡for ¡the ¡purpose ¡of ¡

determining ¡differen5al ¡expression ¡

• How ¡much ¡are ¡expression ¡levels ¡changing ¡

between ¡samples? ¡

RNA-­‑seq ¡

Pepke ¡et ¡al, ¡Nature ¡Methods, ¡2009 ¡

Raw ¡data ¡

• Short ¡sequence ¡reads ¡
• Quality ¡scores ¡

@SEQ_ID GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT + !''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65 ¡

RNA-­‑seq ¡analysis ¡steps ¡

Raw ¡sequence ¡reads ¡ Map ¡onto ¡genome ¡ Summarize ¡reads ¡to ¡transcripts ¡ Sta5s5cal ¡tes5ng: ¡Determine ¡ differen5ally ¡expressed ¡genes ¡ Systems ¡biology ¡ De ¡novo ¡assembly ¡ Annota5on ¡based ¡ Genome ¡guided ¡assembly ¡ Which ¡transcriptome? ¡

Map ¡reads ¡to ¡the ¡genome ¡

• Accuracy ¡vs ¡speed ¡
• Many ¡aligners ¡available ¡(BWA, ¡Bow5e, ¡

Novoalign,…) ¡

CDS ¡ CDS ¡ CDS ¡ CDS ¡ CDS ¡ CDS ¡ CDS ¡ CDS ¡

Gene ¡ transcript ¡

Sequencing ¡transcripts ¡not ¡the ¡genome ¡

Splice ¡site ¡mapping ¡

• Build ¡a ¡junc5on ¡library ¡from ¡

all ¡combina5ons ¡of ¡known ¡ exon ¡boundaries ¡

• Determine ¡where ¡splice ¡

junc5ons ¡occur ¡using ¡the ¡ data ¡itself ¡-­‑ ¡unbiased ¡by ¡

• annota5on. ¡
• Several ¡so]ware ¡packages ¡

to ¡do ¡this ¡such ¡as ¡TopHat, ¡ SplitSeek, ¡SpliceMap… ¡

Gerber ¡et ¡al, ¡Nat ¡Methods, ¡2011 ¡

Which ¡transcriptome ¡to ¡use? ¡

Raw ¡sequence ¡reads ¡ Map ¡onto ¡genome ¡ Summarize ¡reads ¡to ¡transcripts ¡ Sta5s5cal ¡tes5ng: ¡Determine ¡ differen5ally ¡expressed ¡genes ¡ Systems ¡biology ¡

RNA-­‑seq ¡analysis ¡steps ¡

De ¡novo ¡assembly ¡ Annota5on ¡based ¡ Genome ¡guided ¡assembly ¡

Which ¡transcriptome? ¡

Op5on ¡1 ¡

• Use ¡annota5on ¡

– Works ¡well ¡for ¡human ¡and ¡mouse ¡but ¡it’s ¡only ¡as ¡ good ¡as ¡your ¡annota5on ¡ – No ¡novel ¡transcripts ¡are ¡analysed ¡

Op5on ¡2: ¡Genome ¡guided ¡transcript ¡ assembly ¡

• Uses ¡the ¡loca5on ¡and ¡density ¡of ¡reads ¡along ¡

the ¡genome ¡to ¡assemble ¡transcripts ¡

• E.g. ¡Cufflinks ¡
• Can’t ¡assemble ¡across ¡breaks ¡in ¡the ¡genome ¡

– Cancer, ¡poor ¡genomes ¡

Op5on ¡3: ¡De ¡novo ¡transcriptome ¡ assembly ¡

• Assemble ¡transcripts ¡from ¡the ¡data ¡without ¡using ¡a ¡

reference ¡genome ¡

• “Harder” ¡than ¡genome ¡assembly ¡

– Orders ¡of ¡magnitude ¡varia5on ¡in ¡coverage ¡ – Con5gs ¡are ¡short ¡ – Alterna5ve ¡isoforms/transcripts ¡have ¡overlapping ¡sequences ¡ – *Very* ¡computa5onally ¡intensive ¡

• So]ware ¡includes ¡ ¡

– Oases/velvet ¡ – TransAbyss ¡ – Trinity ¡ ¡ – … ¡

Example: ¡Annota5ng ¡the ¡chicken ¡W ¡ chromosome ¡

¡ ¡Z ¡ ¡ ¡ ¡Z ¡ ¡ ¡Z ¡ ¡ ¡W ¡

Male ¡ Female ¡

Two ¡hypotheses ¡for ¡mechanisms ¡of ¡avian ¡sex ¡determina5on: ¡

• 1. Dominant ¡ovary ¡determining ¡gene ¡on ¡W ¡(cf ¡mammals) ¡
• 2. Dosage ¡of ¡Z-­‑linked ¡genes ¡

There ¡is ¡an ¡annotated ¡chicken ¡genome ¡

• Chicken ¡W ¡chromosome ¡is ¡poorly ¡assembled ¡
• Are ¡genes ¡on ¡other ¡chromosomes ¡really ¡on ¡

the ¡W, ¡in ¡par5cular ¡the ¡random ¡chromosome? ¡

Chromosome Assembled Size (Mb) Size inc. random (Mb) Estimated Size (Mb) Estimated Genes (Ensembl) Z 69 70 80 796 W 0.24 0.89 18-54 46 Un_random 56

• 1287

¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡Experimental ¡design ¡

PCR ¡Sexing ¡

+12hour ¡Blastoderms ¡

Hand ¡plate ¡for ¡PCR ¡Sexing ¡

RNA ¡

Stage ¡26 ¡paired ¡gonads ¡(day ¡4.5) ¡

12 ¡Female ¡ ¡ ¡ RNA ¡ 16 ¡Female ¡gonads ¡ Pooled ¡Samples ¡ 12 ¡Female ¡ 16 ¡Female ¡gonads ¡ 12 ¡Male ¡ 12 ¡Male ¡ 16 ¡Male ¡gonads ¡ 16 ¡Male ¡gonads ¡

RNA-­‑seq ¡

• Illumina ¡HiSeq2000 ¡
• Paired-­‑end ¡100bp ¡
• 4 ¡lanes ¡
• >80million ¡reads/

sample ¡

Defining ¡the ¡transcriptome ¡

• Annota5on ¡~20,000 ¡genes ¡
• Genome ¡guided ¡assembly ¡(Cufflinks) ¡~45,000 ¡

genes ¡

• De ¡novo ¡transcriptome ¡assembly ¡~2.5 ¡million ¡

transcripts ¡(Abyss ¡with ¡filtering)! ¡ A ¡combined ¡approach ¡

• Assemble ¡cufflink ¡genes ¡using ¡transcripts ¡from ¡
• ur ¡de ¡novo ¡assembly ¡

Annota5on ¡of ¡the ¡chicken ¡W ¡ combined ¡all ¡three ¡approaches ¡

W/W_random Chromsome Un_random Chromosome Autosomes Blastoderm Coverage Gonads Coverage Abyss Transcripts Cufflinks Transcripts Ensembl Transcripts

1000 1500 2000 2500 3000 3500

base position

391

Coverage Genome Ensembl Cufflinks Abyss

RASA1−W

Full ¡list ¡of ¡W ¡genes/transcripts ¡for ¡differen5al ¡expression ¡

Summariza5on ¡

Take ¡your ¡“transcriptome” ¡and ¡add ¡ up ¡the ¡reads ¡

Coding ¡Sequence ¡ Exons ¡ Introns ¡ Splice ¡Junc5ons ¡

CDS ¡ CDS ¡ CDS ¡ CDS ¡

• ¡Reads ¡in ¡exons ¡
• ¡Exons ¡+ ¡junc5ons ¡
• ¡All ¡reads ¡start ¡to ¡end ¡of ¡transcript ¡
• ¡De ¡novo ¡methods ¡

CDS ¡ CDS ¡ CDS ¡ CDS ¡ CDS ¡ CDS ¡ CDS ¡ CDS ¡ CDS ¡ CDS ¡ CDS ¡

Even ¡when ¡all ¡transcripts ¡are ¡“known” ¡summariza5on ¡or ¡expression ¡ quan5fica5on ¡is ¡difficult. ¡How ¡do ¡you ¡assign ¡reads ¡to ¡transcripts? ¡

Transcript ¡1 ¡ Transcript ¡2 ¡

CDS ¡ CDS ¡ CDS ¡ CDS ¡

Transcript ¡3 ¡

Summariza5on ¡turns ¡mapped ¡reads ¡ into ¡a ¡table ¡of ¡counts ¡

Tag ¡ID A1 A2 B1 B2

ENSG00000124208 478 619 4830 7165 ENSG00000182463 27 20 48 55 ENSG00000125835 132 200 560 408 ENSG00000125834 42 60 131 99 ENSG00000197818 21 29 52 44 ENSG00000125831 ENSG00000215443 4 4 9 7 ENSG00000222008 30 23 ENSG00000101444 46 63 54 53 ENSG00000101333 2256 2793 2702 2976

… … ¡tens ¡of ¡thousands ¡more ¡tags ¡…

** ¡very ¡high ¡dimensional ¡data ¡** ¡

RNA-­‑seq ¡analysis ¡steps ¡

Raw ¡sequence ¡reads ¡ Map ¡onto ¡genome ¡ Summarize ¡reads ¡to ¡transcripts ¡ Sta5s5cal ¡tes5ng: ¡Determine ¡ differen5ally ¡expressed ¡genes ¡ Systems ¡biology ¡ De ¡novo ¡assembly ¡ Annota5on ¡based ¡ Genome ¡guided ¡assembly ¡ Which ¡transcriptome? ¡

Assessing ¡differen5al ¡expression ¡(DE) ¡

• Which ¡genes ¡are ¡changing ¡in ¡their ¡abundance ¡

between ¡samples? ¡

• Sta5s5cal ¡tests ¡for ¡DE ¡(edgeR) ¡

Normaliza5on ¡

Accoun5ng ¡for/removing ¡technical ¡ sources ¡of ¡varia5on ¡

Simple ¡thought ¡experiment ¡

• Two ¡samples ¡A ¡and ¡B, ¡sequenced ¡to ¡the ¡same ¡depth ¡

(same ¡library ¡size) ¡

• Every ¡gene ¡that ¡is ¡expressed ¡in ¡A, ¡is ¡expressed ¡in ¡B ¡

at ¡the ¡same ¡level ¡

• Say ¡there ¡are ¡a ¡small ¡number ¡of ¡genes ¡that ¡are ¡

expressed ¡uniquely ¡to ¡sample ¡B, ¡but ¡they ¡are ¡quite ¡ highly ¡expressed ¡(lots ¡of ¡reads) ¡

• Many ¡genes ¡within ¡the ¡common ¡set ¡will ¡appear ¡

differen5ally ¡expressed ¡(B ¡< ¡A) ¡

Another ¡way ¡to ¡view ¡it ¡

• Hypothe5cal ¡example: ¡

Sequence ¡6 ¡libraries ¡to ¡ the ¡same ¡depth, ¡with ¡ varying ¡levels ¡of ¡unique-­‑ to-­‑sample ¡expression ¡

• Differences ¡in ¡observed ¡

counts ¡among ¡the ¡ common ¡genes ¡

Red=low, ¡goldenyellow=high ¡

Marioni ¡et ¡al. ¡2008 ¡RNA-­‑seq ¡data ¡

Distribu5on ¡of ¡log-­‑ra5os ¡of ¡ counts ¡ 0 ¡counts ¡are ¡omiued ¡ Red ¡line ¡= ¡trimmed ¡mean ¡ trimmed ¡ mean ¡

• 5 ¡lanes ¡of ¡liver ¡RNA, ¡5 ¡lanes ¡of ¡kidney ¡RNA ¡
• Compare ¡two ¡single ¡kidney ¡libraries ¡(technical ¡

replicates), ¡a]er ¡adjus5ng ¡for ¡library ¡size ¡

(a)

log2(Kidney1 NK1) ! log2(Kidney2 NK2) Density

• 6
• 4
• 2

2 4 6 0.0 0.4 0.8 log2(Liver NL) - log2(Kidney NK) Density

• 6
• 4
• 2

2 4 6 0.0 0.2 0.4

(b)

Same ¡sample. ¡ Technical ¡replicates ¡ Different ¡samples. ¡ Liver ¡vs ¡Kidney ¡

Another ¡way ¡to ¡view ¡the ¡data ¡

.

• 20
• 15
• 10
• 5

5 A = log2( Liver NL Kidney NK) M = log2(Liver NL) - log2(Kidney NK)

! !

Housekeeping genes Unique to a sample

(c)

• Sequencing ¡“real ¡

estate” ¡is ¡fixed. ¡

• Underlying ¡RNA ¡

composi5on ¡can ¡be ¡ very ¡different ¡

– e.g. ¡several ¡liver-­‑ specific ¡genes ¡

• An ¡adjustment ¡at ¡

the ¡analysis ¡stage ¡ should ¡be ¡made ¡

The ¡adjustment ¡to ¡data ¡analysis ¡is ¡a ¡ scaling ¡constant ¡

• Assump5on: ¡core ¡set ¡of ¡genes ¡that ¡do ¡not ¡change ¡in ¡
• expression. ¡
• Pick ¡a ¡reference ¡sample, ¡compute ¡TMM ¡rela5ve ¡to ¡reference ¡
• TMM ¡(Trimmed ¡Mean ¡of ¡M ¡values) ¡M=log ¡ra5o ¡
• Adjustment ¡to ¡sta5s5cal ¡analysis: ¡

– Use ¡addi5onal ¡offset ¡(GLM) ¡

Robinson ¡& ¡Oshlack, ¡Genome ¡Biology, ¡2010 ¡

DE ¡analysis ¡happens ¡in ¡R ¡

• Several ¡packages ¡currently ¡available ¡for ¡DE ¡analysis ¡
• f ¡RNA-­‑seq ¡data ¡in ¡R: ¡

– baySeq, ¡ ¡ – DEGSeq, ¡ ¡ – DESeq ¡and ¡ ¡ – edgeR ¡ – NBPSeq ¡ ¡ (+ ¡more ¡methods ¡with ¡less ¡refined ¡code, ¡e.g. ¡TSPM) ¡

• edgeR ¡developed ¡at ¡WEHI, ¡so ¡this ¡is ¡my ¡favourite ¡
• edgeR ¡is ¡only ¡package ¡so ¡far ¡with ¡full ¡GLM ¡

capabili5es ¡

• Great ¡support ¡on ¡Bioconductor ¡mailing ¡list! ¡

Output ¡of ¡analysis ¡

• A ¡list ¡of ¡genes ¡with ¡p-­‑values ¡for ¡tes5ng ¡

differen5al ¡expression ¡between ¡samples ¡

• A ¡ranked ¡list ¡of ¡genes ¡

RNA-­‑seq ¡analysis ¡steps ¡

Raw ¡sequence ¡reads ¡ Map ¡onto ¡genome ¡ Summarize ¡reads ¡to ¡transcripts ¡ Sta5s5cal ¡tes5ng: ¡Determine ¡ differen5ally ¡expressed ¡genes ¡ Systems ¡biology ¡ De ¡novo ¡assembly ¡ Annota5on ¡based ¡ Genome ¡guided ¡assembly ¡ Which ¡transcriptome? ¡

Going ¡beyond ¡gene ¡lists ¡

Gene ¡length ¡bias ¡and ¡ GO ¡analysis ¡

Length ¡bias ¡in ¡RNA-­‑seq ¡

Equal ¡number ¡of ¡transcripts ¡6X ¡length ¡

6X ¡number ¡of ¡fragments ¡ More ¡power ¡to ¡detect ¡DE ¡at ¡a ¡given ¡ threshold ¡

Propor5on ¡of ¡DE ¡genes ¡

Gene length (bases) 2000 4000 6000 8000 Proportion DE 0.00 0.05 0.10 0.15 0.20 0.25 0.30

(a)

Gene ¡Ontology ¡analysis ¡

• Gene ¡Ontology ¡categorises ¡genes ¡into ¡

func5onal ¡groups ¡

• We ¡wish ¡to ¡know ¡if ¡certain ¡gene ¡ontologies ¡

have ¡more ¡DE ¡genes ¡than ¡expected ¡

• If ¡a ¡category ¡has ¡lots ¡of ¡long ¡genes ¡we ¡expect ¡

it ¡to ¡have ¡more ¡DE ¡genes ¡

• Can’t ¡use ¡simple ¡sta5s5cal ¡tests ¡

GO ¡analysis ¡for ¡RNA-­‑seq ¡

¡Develop ¡a ¡computa5onal ¡method ¡for ¡gene ¡ category ¡tes5ng ¡that ¡can ¡account ¡for ¡gene ¡ level ¡biases ¡in ¡differen5al ¡expression ¡detec5on ¡ ¡Category ¡tes5ng ¡refers ¡to ¡tes5ng ¡if ¡a ¡set ¡of ¡genes ¡has ¡an ¡over ¡

representa5on ¡of ¡DE ¡genes ¡

Young ¡et ¡al, ¡Genome ¡Biology, ¡2010 ¡

Three ¡step ¡procedure ¡

• 1. Determine ¡which ¡genes ¡are ¡differen5ally ¡

expressed ¡

• 2. Define ¡a ¡probability ¡weigh5ng ¡func5on ¡
• 3. Generate ¡many ¡random ¡samples ¡to ¡

produce ¡a ¡null ¡distribu5on ¡in ¡order ¡to ¡ calculate ¡significance ¡of ¡a ¡category ¡

Probability ¡weigh5ng ¡func5on ¡

Fit ¡a ¡func5on ¡to ¡the ¡binary ¡data ¡ series ¡1=DE, ¡0=!DE ¡ ¡ We ¡chose ¡a ¡spline ¡with ¡a ¡ monotonicity ¡requirement ¡

Gene length (bases) 2000 4000 6000 8000 Proportion DE 0.00 0.05 0.10 0.15 0.20 0.25 0.30

(a)

Random ¡sampling ¡

• Select ¡a ¡random ¡set ¡of ¡genes ¡the ¡same ¡size ¡

as ¡the ¡set ¡of ¡DE ¡genes ¡

• However ¡the ¡probability ¡of ¡selec5ng ¡a ¡gene ¡

is ¡weighted ¡by ¡the ¡value ¡of ¡the ¡probability ¡ weigh5ng ¡func5on ¡(based ¡on ¡the ¡length ¡or ¡ read ¡count ¡of ¡the ¡gene) ¡

• Then ¡count ¡how ¡many ¡genes ¡in ¡the ¡DE ¡set ¡

have ¡the ¡GO ¡category ¡of ¡interest ¡

• Repeat ¡this ¡many ¡5mes ¡
• Calculate ¡a ¡p-­‑value ¡for ¡the ¡category ¡

GOseq ¡results ¡

Categories ¡with ¡short ¡genes ¡get ¡a ¡ higher ¡rank ¡in ¡GOseq ¡ Categories ¡with ¡long ¡genes ¡get ¡a ¡ lower ¡rank ¡

• 1500
• 1000
• 500

500 1000 Increase in rank in GOseq 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 log10(Average GO category gene length)

(a)

GOseq ¡finds ¡a ¡different ¡set ¡of ¡ significant ¡categories ¡

Even ¡for ¡a ¡small ¡number ¡of ¡ categories ¡there ¡are ¡a ¡significant ¡ number ¡of ¡discrepancies ¡between ¡ the ¡standard ¡method ¡and ¡our ¡gene ¡ length ¡corrected ¡method. ¡ Over ¡all ¡20% ¡different ¡ ¡

andom sampling (200k reps)) 100 200 300 400 Length of lists Number of discrepencies 20 40 60 80

(b)

Hypergeometric Wallenius

RNA-­‑seq ¡analysis ¡steps ¡

Raw ¡sequence ¡reads ¡ Map ¡onto ¡genome ¡ Summarize ¡reads ¡to ¡transcripts ¡ Sta5s5cal ¡tes5ng: ¡Determine ¡ differen5ally ¡expressed ¡genes ¡ Systems ¡biology ¡ De ¡novo ¡assembly ¡ Annota5on ¡based ¡ Genome ¡guided ¡assembly ¡ Which ¡transcriptome? ¡ Learn ¡ something! ¡

The ¡future ¡

• Effec5vely ¡dealing ¡with ¡no ¡genome ¡or ¡poor ¡

quality ¡genomes ¡

• Integra5ng ¡RNA-­‑seq ¡data ¡with ¡other ¡genomics ¡

data ¡

• Analysis ¡methodology ¡is ¡cri5cal ¡and ¡s5ll ¡

developing ¡for ¡specific ¡purposes ¡

• Opportuni5es ¡to ¡use ¡this ¡data ¡in ¡new ¡and ¡

imagina5ve ¡ways ¡– ¡requires ¡new ¡analysis ¡ methodology ¡

Acknowledgements ¡

• Mark ¡Robinson ¡(Uni ¡Zurich) ¡
• Mauhew ¡Young ¡(Cambridge) ¡
• Nadia ¡Davidson ¡(MCRI) ¡
• Mauhew ¡Wakefield ¡(WEHI) ¡
• Gordon ¡Smyth ¡(WEHI) ¡
• Davis ¡McCarthy ¡(Oxford) ¡