Transcriptome Analysis ofNon-Model Organisms
using Short-Read Sequences

Lesley Collins, Patrick Biggs
Claudia Voelckel, Simon Joly
Allan Wilson Centre,
Massey University
New ZealandTranscriptome Analysis
DNARNA

proteinncRNA

Phenotype

Genotype

Which genes are expressed in certain situations?
Which genes are expressed differently in different species?
ESTs (Expressed Sequence Tags) are collected from RNA
samples reflecting which genes are being expressed.

1 December 2008

GIW 2008 Collins and Biggs

1

Expression studies – Of chips and tags
Sample 1 vs Sample 2

Digital Tag Profiling

Microarrays

Tags frommRNA 1Tags frommRNA 2

Illumina
Genome Analyzer 2¾ count tags (tpm 1, tpm 2)¾ calculate tpm2/tpm1 ratios¾ map tags to reference genome

Data analysis
apply statistical tests to cy5/cy3 + tpm2/tpm1 ratios

1 December 2008

GIW 2008 Collins and Biggs

Benefits and constraints of…
Microarray studies&Tag profiling studies¾ microarray platform needed¾ open to any organism (but tag annotation
depends on available genome information)¾ queries limited to genes present on
the array¾ multi-species comparisons tricky¾ any transcript detectable (restriction site)¾ multi-species comparisons straightforward¾ 100 ug per sample¾ 1ug RNA per sample¾ low intensity data unreliable¾ intense data pre-processing¾ single copy transcripts detectable (3 tpm)
¾ little data pre-processing¾ reasonable costs¾ still expensive: cost likely to decrease
when multiplexing biological replicates

1 December 2008

GIW 2008 Collins and Biggs

2

The Pachycladon Transcriptome Project
Alpine Cress (Pachycladon, Brassicaceae)

P. enysii

P. fastigiata

2492 mobligate greywacke

1485 mcan grow on
schist/greywacke

What are the processes of ecological divergence?
Has there been an adaptive radiation to different geological substrates?
1 December 2008

GIW 2008 Collins and Biggs

But Pachycladon is a non-model plant…
There are 8 South Island species of recent origin (<1 mya):
P. fastigiata
P. enysii
P. stellata
P. novae-zealandiae
P. wallii
P. cheesemanii
P. exilis
P. latisiliqua

There is no Pachycladon genome available!

Can we look at the Pachycladon transcriptome without having tosequence the entire genome?

1 December 2008

GIW 2008 Collins and Biggs

3Pachycladon is polyploid!
Member of Brassicaceae familyAllopolyploid from very diverged
Brassicaceae parents

2 copies for most (all?) genes
One from Himalayas
Other from Australia



Rough divergence time
estimates from Arabidopsis
2 – 10 m year (S. Joly)

Simon Joly,submitted
GIW 2008 Collins and Biggs

1 December 2008

Sequencing of 3 Pachycladon species
Platform

TemplateSequence
lengthSolexa

cDNA

36 bp

454

cDNA

40-320 bp

SolexamRNA tags

18 bp

P. enysii

P. fastigiata

XP. cheesemanii

XXXX

1. Construct EST libraries with different sequencing depth using the same platform
2. Compare platforms for EST library construction
3. Establish EST databases/Reference transcriptomes
4. Use EST libraries for development of molecular markers
5. Digital Gene Expression pilot study to be compared with earlier microarray
results
1 December 2008

GIW 2008 Collins and Biggs

4

P. enysii and P. fastigiata – Sampling

Claudia Voelckel and Peter Heenan
1 December 2008

GIW 2008 Collins and Biggs

Some people had much more fun!
Pete Lockhart collecting samples from Molesworth station!!

1 December 2008

GIW 2008 Collins and Biggs

5Pachycladon transcriptome analysis
Collect reads and assemble into contigs

Pachy spp. ESTs

454 and Solexa

Compare contigs/reads against Arabidopsis

Which genes are easily identified?

Map tags against Arabidopsis

Preliminary expression analysis

Annotate contigs into ‘potential ESTs’ and map tags against thesePachycladon ‘ESTs’

Pachycladon Expression Analysis !!!!

1 December 2008

GIW 2008 Collins and Biggs

Pachycladon data collected

63,870 FLX-454 reads for P. fastigiata
78 contigs (length > 500 nt)
3825 contigs (length > 100 nt)

40 million Solexa 36nt reads for P. enysii
10 million Solexa 36nt reads for P.
cheesemanii14.9 million 18 nt Digital tags for P. enysii

15.8 million 18 nt Digital tags for P. fastigiata

1 December 2008

Otago Genomics
Sequencer,
Otago University,
New Zealand

AWC Genome
Sequencing Service,
Massey University,
New Zealand

GIW 2008 Collins and Biggs

6

Short-Read36nt
sequences

Sequencing

Construction of in silico genomeof concatenated ESTs fromclosest reference genome18 nt
Digital Tags

De novo assembly

Mapping

Mapping ofcomplete data against
in silico genome

De novo assembly oncomplete data

Digital Taganalysis

Sequence length
optimisation using
partial data

Analysis onComplete data

BLAST contigs against
in silico genome

Database sorting of results
and visualisation

Expression
Analysis

SNP
analysis

Duplicate genesortingComparisons againstreference genome ESTs

Genomic Analysis

Optimised mapping against Arabidopsis TAIR7 ESTs
100%

Evolutionary distance between
Pachycladon and Arabidopsis

90%
80%
70%

% total results

Data was mapped with Eland (Illumina)
 Allows only 2 mismatches per sequence Not all data will map

R0R1R2NMU0U1U260%
50%
40%
30%20%
10%0%
17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32

Sequence length

ELAND
resultELAND
Result type

Length = 19(1 lane)% Total(1 lane)Number(7 lanes)

Av % Total(7 lanes)

Std Dev(7 lanes)

U0

Unique 0 mismatch

341,804

8.99

3,577,678

8.86

1.2

U1

Unique 1 mismatch

543,248

14.29

5,697,913

14.23

0.34

U2

Unique 2 mismatch

520,442

13.695,547,612

13.85

0.29

R0

Repeat 0 mismatch

960,11925.26

1,943,961

24.84

0.93

R1

Repeat 0 mismatch

191,895

5.05

2,914,570

4.88

0.22

R2

Repeat 0 mismatch

275,036

7.24

10,455,270

7.27

0.16

QC

Quality filtered

14

0.00

121

0.00

0.00

NMNo match

968,84225.49

9916563

26.07

0.79

Total

-

3,801,400

100

40,053,567

100-

1 December 2008

GIW 2008 Collins and Biggs

7

Optimised Contig assembly parameters
Edena

1.00.9Hernandez et al. 2008 GenomeRes 18:802

Velvet

Zerbino and Birney 2008
Genome Res 18:821

0.8
Cumulatve Fraction

0.7
0.6

all conditions (n = 464642)

0.5

kmer = 15 (n = 108528)

0.4

kmer = 17 (n = 144532)
kmer = 19 (n = 90736)

0.3

kmer = 21 (n = 59058)

0.2

kmer = 23 (n = 38353)
kmer = 25 (n = 23435)

0.1
0.0
10100

1000

Assembled Contig Size (nt)

K-mer
Size

Number ofContigs

BLAST
hits

Number ofContigs with
BLAST hits

% Contigs
with BLAST
hits

Number of genes hit
All aligns

Aligns>40 ntAligns>85 nt

15231

50.00

3201712,638

18,844

8,780

69.50

10,579

9,312

6,786

19

22,631

38,535

16,41372.50

14,906

13,105

10,30421

20,531

39,554

15,22774.20

14,594

12,2679,504

23

16,873

34,486

12,73975.50

13,20910,732

8,222

25

12,75027,360

9751

76.50

10,871

8,852

6,595

Total

85,42562,91173.65

1 December 2008

GIW 2008 Collins and Biggs

Optimised contig mapping parameters

Number of Arabidopsis ESTs to which Pachycladon contigs mappedwith any overlap (bin size = 20 nt)

E.g. These contigs mapped to 617 Arabidopsis ESTs with anoverlap of 1-20 nt

1 December 2008

GIW 2008 Collins and Biggs

8Arabidopsis TAIR7 vrs TAIR8This project began using TAIR7…
…then TAIR8 was released
 TAIR7


27,029 protein coding genes (3866 have splice variants - 12%)TAIR827,235 protein coding genes (4330 have splice variants - 13%)1291 new genes
2009 new gene models
23% of existing TAIR7 genes modified

GIW 2008 Collins and Biggs

1 December 2008

Mapping against the Arabidopsis TAIR8 Genome• Arabidopsis transcriptome: 38,963 ESTs (27,235 protein coding genes)• With Eland reads map uniquely (U), repeatedly (R) or not at all (NM & QC)
100%% total

90%
80%

R2count70%

R1count60%

R0count
QCcount

50%

NMcount

40%

Mapping results for1 lane (3.8 million reads)
with 20mers:

U2count

30%

U1count

20%

U0count

10%0%
15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32

NM27%

R

29%

U

45%

eland

Optimal k-mer: 20
Conclusion: TAIR8 is slightly better but TAIR7 data still very useful
1 December 2008

GIW 2008 Collins and Biggs

9Now onto the full analysis
Mapping ofcomplete data toTAIR8

De novo assemblyof Pachycladon
ESTs

Analysis of DigitalTags

GIW 2008 Collins and Biggs

1 December 2008

Contigs assembled using Velvet and Edena


Criteria
BLASTN HSP alignment must cover at least 90% of the contig length
Alignment must contain at least 90% identical nucleotides
P. enysii
Edena

P. enysii
VelvetP. cheesemanii
Edena

P. cheesemanii
VelvetContigs >500nt

542

1042

365

573

Hits inArabidopsis

971256872

Contigs from 454 for P. fastigiata reanalysedEST database Version 0.1 under construction

Many thanks to Oliver Deusch and Nicole Grünheit

1 December 2008

GIW 2008 Collins and Biggs

10

1 December 2008

GIW 2008 Collins and Biggs

Gbrowse against Arab Tair8

1 December 2008

GIW 2008 Collins and Biggs

11

Finding Pachycladon gene pairs
GGA

Most genes will have two copies when compared to Arabidopsis

GA
A
E.g. Two overlapping Pachycladon contigs (114bp) that map to At5g64740.1, cellulose synthase 6

GIW 2008 Collins and Biggs

1 December 2008

Digital Tag Analysis:
P. enysii and P. fastigiata against ArabidopsisPe - 5,591,951 sequences (37.64%)
mapped uniquely
Pf - 5,767,796 sequences (36.48%)
mapped uniquely
1/3 of all sequences mapped repeatedly
1/3 of all sequences did not map

P. enysii
Genes

U0U012Mapped

15296

31375

tpm > 3

1311

19451

Occur in both species

1068

18263

Occur only in P. enysii

243

1188

P. fastigiata
243, 15%

283, 18%38, 2%

48, 3%

987, 62%

1 December 2008

Only in Pe

Genes

Up in Pe

Mapped

16204

31511

tpm > 31351

20753

Not different

Occur in both species

1068

18263

Up in PfOccur only in P.
fastigiata

283

2490

U0U012

Only in Pf
Lydia Hopp and Janina Mothes

GIW 2008 Collins and Biggs

12

P. enysii vs. P. fastigiata – Previous Chip study
P. fastigiata (1485 m, glabrous)

310 genes
P. enysii (1885 m, hairy)

324 genesGene expression differences predict
differences in secondary metabolites(e.g. type of mustard oils produced)
Follow up studies confirm: P. enysii
leaves less toxic
herbivory important for species
differentiation !

In press: Molecular Ecology!
Transcriptionalandbiochemicalsignatures of divergence in natural
populations of two species of New
Zealand alpine Pachycladon

C. Voelckel1*, P. B. Heenan2, B.
Janssen3, M. Reichelt4, K. Ford2, R.
Hofmann5, P. J. Lockhart1

We can now compare our Digital Tag data with this chip study!
GIW 2008 Collins and Biggs

1 December 2008

SummaryPachycladon provides many challenges for genomics:
Non-model plant

Using Arabidopsis as a reference genome worked to a large extent
But we need to construct a Pachycladon EST database for moremeaningful evolutionary studiesPolyploidMultiple copies of the same geneCauses problems with contig assembly

Alternative splicing

1 December 2008

Multiple transcripts of each gene copy
Can cause problems with mapping

GIW 2008 Collins and Biggs

13

And still more to come…
Preliminary Pachycladon EST libraries in assembly
Pachycladon enysii
Pachycladon fastigiata
Pachycladon cheesemanii

Comparison of Digital Tag data to Arabidopsis Chip study

Analysis of Digital Tag data with Pachycladon ESTs

SNP analysis of Pachycladon data
Between Pachycladon species
Between Pachycladon and Arabidopsis

1 December 2008

GIW 2008 Collins and Biggs

AcknowledgementsPatrick Biggs
Claudia Voelckel
Simon Joly
Pete Lockhart

Oliver Deusch and Nicole Grünheit (Bill Martin’s lab,Duesseldorf))
Lydia Hopp and Janina Mothes (Greifswald University)



AWC Genome Sequencing Service

Lorraine Berry
Maurice Collins





Allan Wilson Centre (David Penny)
Institute of Molecular BioSciences
Marsden FundRoyal Society of New Zealand
Humboldt Foundation

1 December 2008

GIW 2008 Collins and Biggs

14