全 文 ：基于 Illumina RNA鄄Seq短序列的转录组从头组装软件比较与优化*
赵摇 磊, Zachary LARSON鄄RABIN, 陈斯云, 郭振华**
(中国科学院昆明植物研究所中国西南野生生物种质资源库, 云南 昆明摇 650201)
摘要: 利用公共数据库中果蝇 F1 代和栽培水稻基于高通量 Illumina 测序平台的 RNA鄄Seq 短序列数据, 比
较了 8 个 (ABySS, Velvet, SOAPdenovo, Oases, Trinity, Multiple鄄k, T鄄IDBA and Trans鄄ABySS) 转录组从头组
装软件。 结果显示, 在基于单一 k鄄mer和多重 k鄄mer方法的两类软件中, Trinity和 Trans鄄ABySS分别表现出
最好的组装性能, 而其它软件性能比较接近。 我们还发现基于多重 k鄄mer 比单一 k鄄mer 可以组装获得更多
的总碱基数目, 但是即使利用最好的多重 k鄄mer 组装软件, 所获得的数据质量也比研究人员所期望的要
低。 鉴于此, 我们提出了 “ETM冶 优化方法, 将多重 k鄄mer方法组合到 Trinity中, 使其在具有最好的组装
性能的基础上兼具了多重 k鄄mer的优势, 测试结果显示了该方法具有一定的优越性。 我们的研究结果为用
户选择合适的软件提供了依据, 对推动基于高通量 Illumina测序的转录组研究具有重要意义。
关键词: 高通量; 二代测序; 转录组; 从头组装; 优化
中图分类号: Q 75 摇 摇 摇 摇 摇 摇 摇 文献标识码: A摇 摇 摇 摇 摇 摇 摇 文章编号: 2095-0845(2012)05-487-15
Comparing De Novo Transcriptome Assemblers
Using Illumina RNA鄄Seq Reads
ZHAO Lei, Zachary LARSON鄄RABIN, CHEN Si鄄Yun, GUO Zhen鄄Hua**
(Plant Germplasm and Genomics Center, Germplasm Bank of Wild Species, Kunming Institute of Botany,
Chinese Academy of Sciences, Kunming 650201, China)
Abstract: In this study, we carried out a systematic comparison of the de novo transcriptome assembly performance
of eight assemblers (ABySS, Velvet, SOAPdenovo, Oases, Trinity, Multiple鄄k, T鄄IDBA and Trans鄄ABySS), pro鄄
cessing Illumina RNA鄄Seq reads from F1 hybrids (Drosophila MS) of Drosophila melanogaster and Drosophila sech鄄
ellia and cultivated rice. Our study showed that Trinity and Trans鄄ABySS were the most effective for producing tran鄄
scriptomes from our trial datasets using single k鄄mer and multiple k鄄mer methods, respectively, although the per鄄
formance levels of the other tested assemblers were comparable. We found that using single k鄄mer assemblers gener鄄
ally produced fewer total numbers of bases than multiple k鄄mer assemblers, although even the best assembler爷s re鄄
sults showed lower quality than some researchers may desire. Therefore, we developed and tested a novel de novo
transcriptome assembly method, ETM, which employs a combination of multiple k鄄mer tools with Trinity assembler.
The ETM method yielded superior results from our trial datasets. Our results will assist the growing number of tran鄄
scriptome projects using Illumina RNA鄄Seq reads and provide guidelines for choosing appropriate software.
Key words: High鄄throughput; NGS; Transcriptome; de novo assembly; Optimized
植 物 分 类 与 资 源 学 报摇 2012, 34 (5): 487 ~ 501
Plant Diversity and Resources摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 DOI: 10. 3724 / SP. J. 1143. 2012. 12084
*
**
Foundation items: The National Natural Science Foundation of China (30990244), the Knowledge Innovation Project of the Chinese Acade鄄
my of Sciences (KSCX2鄄YW鄄N鄄067), the Scientific Research Foundation for the Returned Overseas Chinese Scholars,
State Education Ministry and the Young Academic and Technical Leader Raising Foundation of Yunnan Province
(2008PY065), a Chinese Academy of Sciences Young International Scientist Fellowship ( awarded to Zachary Larson鄄
Rabin), Yunnan Provincial Government through an innovation team program
Author for correspondence; E鄄mail: guozhenhua@ mail. kib. ac. cn
Received date: 2012-06-01, Accepted date: 2012-06-19
作者简介: 赵摇 磊 (1978-) 男, 硕士, 从事生物信息与比较基因组学。 E鄄mail: zhaolei@ mail. kib. ac. cn
摇 A transcriptome is the set of all RNA molecules
present in a given cell or tissue at a given time, and
the sequencing of a transcriptome is vital to fully un鄄
derstand the genetics and biology of that cell, tissue,
or organism. Although transcriptome sequencing and
analysis have been possible for more than 20 years,
recent technological advances have made the process
much faster and cheaper. It has become feasible to
quickly assemble transcriptomes and quantify tran鄄
script abundance in species for which no genome se鄄
quence has yet been constructed (i. e., de novo tran鄄
scriptome assembly), using so鄄called next鄄generation
sequencing technologies (NGS). There are currently
three types of NGS platforms, including Roche 454
FLX, Illumina, and Applied Biosystems SOLiD. The
Roche 454 platform generates reads of about 330 bp,
and such relatively long read lengths are useful for de
novo assembly since they improve the accuracy of
mapping repetitive regions. Unfortunately, it is more
expensive than the other two platforms and produces
higher error rates in homopolymer repeats (Metzker,
2010; Shendure and Ji, 2008). The ABI SOLiD plat鄄
form is economical, but generates reads that are often
too short to be well鄄suited for de novo assembly (An鄄
sorge, 2009). The Illumina system is also relatively
inexpensive, but generates moderate read lengths
(75 - 100 bp) that suffice for de novo assembly.
Thus, it is the most widely used platform for NGS
and has been recommended for de novo assembly
(Fullwood et al., 2009; Simon et al., 2009).
The “RNA鄄Seq冶 method of transcriptome ana鄄
lysis makes use of the high numbers of short reads
produced by the NGS platforms. RNA鄄Seq employs
a “ shotgun 冶 鄄style approach both to sequence a
transcriptome and to quantify the abundance of its
transcripts, obviating the need for BAC and cDNA
cloning ( Wang et al., 2009; Haas and Zody,
2010). The RNA鄄Seq method starts with the high鄄
throughput sequencing of double鄄stranded cDNA
fragments, followed either by direct mapping of the
sequences to a reference genome, or by their de novo
assembly. This results in a comprehensive transcrip鄄
tome analysis that includes transcripts爷 untranslated
regions ( UTRs) and alternative鄄splicing isoforms,
as well as non鄄coding RNAs (ncRNAs). Since the
abundance of each transcript is quantified, RNA鄄Seq
can be used to assess the differential expression of
genes in both physiological and pathological condi鄄
tions, or to compare the transcriptomes of different
tissue types, developmental stages, or genotypes
(Nagalakshmi et al., 2008; Mortazavi et al., 2008).
The use of a reference genome to reconstruct a
transcriptome is termed ‘genome鄄guided assembly爷,
while ‘de novo assembly爷 attempts transcriptome as鄄
sembly without such aid (Garber et al., 2011). Ge鄄
nome鄄guided approaches, such as Cufflinks (Trap鄄
nell et al., 2010) and Scripture (Guttman et al.,
2010), first map all of the reads to a reference ge鄄
nome, and then merge the overlapping reads into
transcripts. However, the increasing interest in re鄄
searching non鄄model organisms renders de novo as鄄
sembly approaches increasingly important. In recent
years, many programs have been developed specific鄄
ally for de novo assembly of RNA鄄Seq reads genera鄄
ted by Illumina (Table 1).
In this study, we compared eight free de novo
assemblers to assist the growing number of transcrip鄄
tome projects using Illumina RNA鄄Seq reads (Table
2): ABySS ( Simpson et al., 2009; Birol et al.,
2009), Velvet (Zerbino and Birney, 2008), SOAP鄄
denovo (Li et al., 2009), Oases (D. R. Zerbino,
European Bioinformatics Institute), Trinity ( Grab鄄
herr et al., 2011), Multiple鄄k ( Surget鄄Groba and
Montoya鄄Burgos, 2010 ), T鄄IDBA ( Peng et al.,
2011) and Trans鄄ABySS (Robertson et al., 2010),
each of which is based on De Bruijn graph algorithms
(Miller et al., 2010; Zerbino and Birney, 2008).
Our study did not compare commercial software such
as SeqMan (Santure et al., 2011), CLC (Wickra鄄
masinghe et al., 2011; Kumar and Blaxter, 2010),
or Rnnotator (Martin et al., 2010), although these
programs can also process Illumina reads. We uti鄄
lized the publicly available Illumina RNA鄄Seq reads
from F1 hybrids (Drosophila MS) of D. melanogaster
884摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 植 物 分 类 与 资 源 学 报摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 第 34 卷
Table 1摇 Software programs previously used for de novo transcriptome assembly using Illumina RNA鄄Seq reads
Tool Organism
ABySS Tumor (Birol et al., 2009); Conus bullatus (Hu et al., 2011); Burmese pythons (Wall et al., 2011)
Velvet
Chinese hamster ovary (Birzele et al., 2010); Anopheles funestus (Crawford et al., 2010); Eucalyptus tree (Mizrachi et al.,
2010); Pachycladon enysii (Collins et al., 2008); Aedes aegypti and Anopheles gambiae (Gibbons et al., 2009); Medicago
sativa (Yang et al., 2011); Potato (Hamilton et al., 2011); Ruditapes philippinarum (Ghiselli et al., 2011)
SOAPdenovo
Whitefly (Wang et al., 2010b); Sweetpotato (Wang et al., 2010c); Nilaparvata lugens (Xue et al., 2010); Migratory locust
(Chen et al., 2010); Lateolabrax japonicus (Xiang et al., 2010); Camellia sinensis (Shi et al., 2011); Salvia miltiorrhiza
(Wenping et al., 2011); Populus euphratica (Qiu et al., 2011); Dugesia japonica (Qin et al., 2011); Siraitia grosvenorii
(Tang et al., 2011); Acacia auriculiformis and Acacia mangium (Wong et al., 2011); Taxus (Hao da et al., 2011)
Oases
Chickpea (Garg et al., 2011); Sheep (Jager et al., 2011); Eichhornia paniculata (Ness et al., 2011); Radix balthica (Feld鄄
meyer et al., 2011); Xiphophorus maculatus (Garcia et al., 2011); Pandalus latirostris (Kawahara鄄Miki et al., 2011)
Trinity Schizosaccharomyces pombe (Grabherr et al., 2011)
Multiple鄄k Loricaria (Surget鄄Groba and Montoya鄄Burgos, 2010)
T鄄IDBA Embryonic stem cells (Peng et al., 2011)
Trans鄄ABySS Liver (Robertson et al., 2010)
Table 2摇 Characteristics of de novo RNA鄄Seq assembly programs
Assembler Strategy Cost k鄄mer parameter URL
ABySS DBG Free Single k鄄mer www. bcgsc. ca / platform / bioinfo / software / abyss
Velvet DBG Free Single k鄄mer www. ebi. ac. uk / ~ zerbino / velvet /
SOAPdenovo DBG Free Single k鄄mer http: / / soap. genomics. org. cn / soapdenovo. html
Oases DBG Free Single k鄄mer www. ebi. ac. uk / ~ zerbino / oases /
Trinity DBG Free Single k鄄mer http: / / trinityrnaseq. sourceforge. net
Multiple鄄k DBG Free Multiple k鄄mer www. surget鄄groba. ch / download / stm. tar. gz
T鄄IDBA DBG Free Multiple k鄄mer http: / / i. cs. hku. hk / ~ alse / hkubrg / projects / tidba /
Trans鄄ABySS DBG Free Multiple k鄄mer www. bcgsc. ca / platform / bioinfo / software / trans鄄abyss
Abbreviations: DBG=De Bruijn graph; K鄄mer=a sequence “word冶 of k nucleotides in length
and D. sechellia, and cultivated rice transcriptomes
(McManus et al., 2010; Zhang et al., 2010) to ex鄄
ploit the comprehensiveness of the genomic and pro鄄
teomic annotation of these model species. Our goals
were to gather and compare assembler performance
data, including requirements for time and RAM
(Random Access Memory) of the computers; the
maximum, mean and median sizes of the tentative
consensus transcript sequences ( TCTSs) and their
length distributions; as well as the N50 and N90
summary statistics of each dataset. We also compared
the assembly accuracy, integrity, coverage and mis鄄
matches ratios resulting from each program, including
the percentages of reads successfully mapped back to
the assemblies, and the efficiency of de novo tran鄄
scriptome assemblers in dealing with alternative spli鄄
cing. We compiled the various factors and used these
data formulate suggestions for researchers to consider
when choosing appropriate de novo transcriptome as鄄
semblers. Finally, we designed a new workflow for
optimized assembly using multiple programs.
Materials and methods
Transcriptome reads and reference datasets
RNA鄄Seq paired鄄end reads of the Drosophila MS
and cultivated rice transcriptomes produced on the Il鄄
lumina sequencing platform were obtained from the
NCBI Gene Expression Omnibus (http:/ / www. ncbi.
nlm.nih.gov / geo) under accession numbers GSE20421,
GSE16631 and GSE16507 and used for de novo tran鄄
scriptome assembly (McManus et al., 2010; Zhang
et al., 2010). There were 4 019 544 876 (4G) 37鄄
bp paired鄄end sequence reads, 5 982 700 350 (6G)
75鄄bp paired鄄end sequence reads, with sequence
depths of 25伊 and 15 伊, in the Drosophila MS and
cultivated rice datasets, respectively. Genome, pro鄄
9845 期摇 摇 摇 摇 ZHAO Lei et al. : Comparing De Novo Transcriptome Assemblers Using Illumina RNA鄄Seq Reads摇 摇 摇 摇 摇
tein, and gene datasets for use as reference se鄄
quences were obtained from FlyBase (www. flybase.
org) for D. melanogaster and from MSU Rice and
RAP ( http: / / rapdb. dna. affrc. go. jp / download / in鄄
dex. html) for Oryza sativa.
Assembly parameters and extraction
For each single k鄄mer assembler, different series
of k鄄mer parameter were set (k=23, 25, 27, 29, 31)
(see Additional File 1, Table S1, http: / / journal.
kib.ac.cn / UserFiles / File / 2012鄄084_editing(1).rar).
For each multiple k鄄mer assembler, incremental mul鄄
tiple k鄄mer parameter combinations were assigned (k
=(23, 25), (23, 25, 27), (23, 25, 27, 29),
(23, 25, 27, 29, 31)) (see Additional File 1, Ta鄄
ble S1, http: / / journal. kib. ac. cn / UserFiles / File /
2012鄄084_editing (1). rar). Each assembled TCTS
(逸100 bp in length) from each program and dataset
was extracted and analyzed by using the mathematical
statistical method contained in our Perl script targe鄄
ting TCTS length distributions (Tables 3 and 4; Ad鄄
ditional File 2, http: / / journal. kib. ac. cn / UserFiles /
File / 2012鄄084_editing(1). rar). The best assembly
produced by each assembler was then chosen for com鄄
parison.
All programs were run on a supercomputer com鄄
prising 16 Quad鄄Core AMD 2. 2 GHz CPUs with 128
GB memory. The X_86 64 bit operating system used
SUSE Linux Enterprise Server 10.
Mapping and alignment
To obtain the percentage of reads utilized in
each assembly, Bowtie was used to map all reads
back to each transcriptome (Tables 3 and 4). All pa鄄
rameters were set to 鄄n 2 鄄l 28 鄄e 70 (cf. Bowtie in鄄
struction manual) (Langmead et al., 2009). To de鄄
termine accuracy and integrity statistical estimates,
the NCBI Standalone BLAST 2. 2. 25 was used to a鄄
lign each assembly to the D. melanogaster and O. sa鄄
tiva protein sequence sets, using an E鄄value thresh鄄
old of 1e鄄5 (Altschul et al., 1990; Mount, 2007).
For the alignment of each query sequence, only the
best鄄scoring hit was used, then mismatches rate ( to鄄
tal mismatched bases / total de novo assembled bases)
was computed using Perl script ( see Additional file
6, http: / / journal. kib. ac. cn / UserFiles / File / 2012鄄
084_ editing (1). rar) . Assembly accuracy (% of
TCTSs hit) and integrity (% of proteins hit) (Ta鄄
bles 5 and 6) were then calculated using Perl script
( Additional file 3, http: / / journal. kib. ac. cn /
UserFiles / File / 2012鄄084_editing(1). rar) . For the
percentages of TCTSs achieving full coverage, the
UCSC Standalone BLAT v. 34 (Kent, 2002) was
used to align all TCTSs to the D. melanogaster and
O. sativa gene sequences, and the maximum intron
size was set to 3 000 bp (Deutsch and Long, 1999;
Atambayeva et al., 2008). In the alignment of each
query sequence, only the best鄄scoring hit was used,
Table 3摇 Statistical characteristics for TCTSs assembled by each program using Drosophila MS RNA鄄Seq reads
AB Ve SO Oa Tr Mk TI TA
TCTS逸100 bp 69 477 60 112 66 407 46 929 21 662 126 212 90 202 164 798
TCTS臆200 bp 44 005 38 241 36 358 13 916 0 72 539 69 971 48 407
TCTS逸1 kb 801 432 7 218 17 358 4 748 1 294 125 11 079
Max 3 800 2 873 10 974 22 442 8 519 2 873 2 734 4 326
Mean 232 221 432 1 065 773 243 174 425
Median 160 163 177 604 556 178 137 318
N50 273 246 932 2 081 943 285 173 583
N90 117 118 140 551 379 124 108 200
% of reads mapped 28. 9 22. 3 31. 3 56. 6 78. 6 46. 3 27. 5 64. 3
Multi鄄hit reads 逸5 0 0 26 399 8 557 360 0 0 0 6 902 154
Abbreviations: AB: ABySS; Ve: Velvet; SO: SOAPdenovo; Oa: Oases; Tr: Trinity; Mk: Multiple鄄k; TI: T鄄IDBA; TA: Trans鄄ABySS. First row
reports assembler names; Rows 2 to 9 report the TCTS size statistics, specifically: the number of TCTSs with size 逸100 bp, the number of TCTSs
with size 臆200 bp, and the number of TCTSs with size 逸1 kb. The remaining statistics represent the following TCTS characteristics: maximum sizes,
mean sizes, median sizes, N50 and N90 sizes, the percentages of reads mapped back to the assemblies, and the numbers of multi鄄hit reads
094摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 植 物 分 类 与 资 源 学 报摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 第 34 卷
and mismatches rates were then calculated using
Perl script ( see Additional file 6, http: / / journal.
kib. ac. cn / UserFiles / File / 2012鄄084 _ editing ( 1 ).
rar) . The score of each alignment was computed u鄄
sing the formula: score = matches鄄mismatches. The
same method was used when aligning gene models to
Table 4摇 Statistical characteristics for TCTSs assembled by each program using cultivated rice RNA鄄Seq reads
AB Ve SO Oa Tr Mk TI TA
TCTS逸100 bp 101 295 143 289 112 246 120 504 86 492 220 546 155 501 366 812
TCTS臆200 bp 47 354 86 649 42 312 31 453 0 129 214 91 962 144 556
TCTS逸1 kb 6 264 1 305 17 679 29 346 26 221 2 037 3 592 34 867
Max 9 616 5 184 12 950 14 730 11 946 5 184 5 539 12 347
Mean 358 230 562 735 914 235 255 437
Median 213 173 270 387 626 178 172 249
N50 518 254 1 038 1 400 1 247 261 304 695
N90 151 123 207 293 410 125 119 177
% of reads mapped 38. 7 22. 2 36. 8 59. 8 75. 4 25. 0 23. 2 54. 4
Multi鄄hit reads逸5 0 0 3 195 5 968 148 0 0 0 5 740 958
Abbreviations: AB: ABySS; Ve: Velvet; SO: SOAPdenovo; Oa: Oases; Tr: Trinity; Mk: Multiple鄄k; TI: T鄄IDBA; TA: Trans鄄ABySS. First row
reports assembler name; Rows 2 to 9 report the TCTS size statistics, specifically: the number of TCTSs with size 逸100 bp, the number of TCTSs with
size 臆200 bp, and the number of TCTSs with size 逸1 kb. The remaining statistics represent the following TCTS characteristics: maximum sizes,
mean sizes, median sizes, N50 and N90 sizes, the percentages of reads mapped back to the assemblies, and the numbers of multi鄄hit reads mapped
back to the assemblies
Table 5摇 Accuracy, integrity and mismatches ratio of TCTSs obtained by each assembler using Drosophila MS RNA鄄Seq reads
AB Ve SO Oa Tr Mk TI TA
Accuracy: % of TCTSs hit to a known
D. melanogaster protein sequence 75. 8 75. 3 75. 4 76. 5 92. 5 78. 6 78. 7 81. 6
% of TCTSs showing 80% identity to a
known D. melanogaster protein sequence 72. 9 74. 5 73. 3 73. 8 91. 8 78. 0 77. 8 80. 8
Integrity: % of D. melanogaster protein
sequences with a TCTS hit 84. 6 83. 3 85. 2 87. 1 76. 9 86. 5 83. 1 86. 1
% of D. melanogaster protein sequences
showing 80% identity with a TCTS 70. 7 69. 6 71. 2 73. 6 58. 3 73. 6 69. 1 72. 6
Mismatches ratio / % 0. 30 0. 31 0. 31 0. 53 0. 35 0. 35 0. 34 0. 31
Notes: E鄄value cutoff =1e鄄5. To measure the accuracy of each assembly program, 23 711 D. melanogaster proteins were compared to the TCTSs assem鄄
bled by each program using BLASTX. To measure integrity, the same D. melanogaster protein sequences were compared to de novo assembled TCTSs
using TBLASTN. Mismatches ratio is computed by total mismatched bases / total de novo assembled bases
Table 6摇 Accuracy, integrity and mismatches ratio of TCTSs obtained by each assembler using cultivated rice RNA鄄Seq reads
AB Ve SO Oa Tr Mk TI TA
Accuracy: % of TCTSs hit to a known
O. sativa protein sequence 54. 8 49. 8 44. 7 62. 9 82. 1 46. 9 46. 5 68. 5
% of TCTSs showing 80% identity to
a known O. sativa protein sequence 50. 4 46. 2 39. 9 56. 5 74. 1 43. 2 43. 3 63. 8
Integrity: % of O. sativa protein
sequences with a TCTS hit 81. 7 81. 4 84. 5 84. 0 87. 6 83. 4 81. 7 83. 5
% of O. sativa protein sequences
showing 80% identity with a TCTS 44. 6 44. 3 46. 8 49. 0 46. 6 47. 0 44. 0 48. 9
Mismatches ratio / % 0. 65 0. 64 0. 51 0. 76 0. 78 0. 67 0. 64 0. 68
Notes: E鄄value cutoff =1e鄄5. To measure the accuracy of each assembly program, 67 393 O. sativa proteins were compared to the TCTSs assembled by
each program using BLASTX for accuracy. To measure integrity, the same O. sativa protein sequences were compared to de novo assembled TCTSs u鄄
sing TBLASTN. Mismatches ratio is computed by total mismatched bases / total de novo assembled bases
1945 期摇 摇 摇 摇 ZHAO Lei et al. : Comparing De Novo Transcriptome Assemblers Using Illumina RNA鄄Seq Reads摇 摇 摇 摇 摇
all TCTSs to determine the number of genes showing
full coverage. Percentages of TCTSs and genes that a鄄
chieved full and partial coverage were calculated
(Tables 7 and 8) with Perl script (see Additional file
4, http: / / journal. kib. ac. cn / UserFiles / File / 2012鄄
084_editing(1).rar).
Analysis of alternative splice sites
Current GTF annotations for D. melanogaster
and O. sativa were downloaded from EnsemblGe鄄
nomes (ftp: / / ftp. ensembl. org / pub / release鄄64 / gtf /
drosophila _ melanogaster / Drosophila _ melanogaster.
BDGP5. 25. 64. gtf. gz and Oryza_sativa. MSU6. gtf.
gz) . De novo assembly transcripts were mapped to ge鄄
nome sequences using GMAP (http:/ / research鄄pub.
gene. com/ gmap / ) (Wu and Watanabe, 2005), show鄄
ing only the single best scoring alignment, and splice
levels were detected as recommended for GSNAP
(GMAP manual). Novel and known splice sites were
reported from each of the alignments (Fig. 5).
Results
Time and RAM consumptions
Each de novo RNA鄄Seq assembly program was
initially scored for required RAM consumption and
for the time needed to complete assembly calcula鄄
tions. Both factors depend on the complexity of the
dataset and the efficiency of the assembly algorithm.
For both the Drosophila MS and the cultivated rice
datasets, SOAPdenovo proved most efficient both in
terms of processing speed and memory usage. In
contrast, Trinity consumed the most computational
resources, while the other assemblers displayed in鄄
termediate requirements (Figures 1 and 2). Howev鄄
er, it should be noted that analysis efficiency was
not linked to assembly accuracy.
Comparisons of TCTS contiguity produced by
the eight programs
摇 摇 The eight de novo assembly programs produced
very different TCTS populations, and a variety of
summary statistics and graphical representations have
been used to describe their size characteristics, or
contiguity. The assembly programs varied widely in
the size distributions of the TCTSs themselves, while
also varying between results for the Drosophila MS
and cultivated rice datasets (Figures 3 and 4). For
instance, the long and relatively thick tails of TCTS
Table 7摇 Summaries of D. melanogaster gene coverage, the percentages of TCTSs covering 95% of TCTS
matched gene sequences and mismatches ratio for each de novo assembly program
AB Ve SO Oa Tr Mk TI TA
% of TCTSs achieving TCTS coverage 31. 2 36. 9 18. 5 40. 0 72. 9 42. 4 18. 4 55. 8
Number of genes achieving total gene coverage 31 9 127 1 097 549 38 6 146
Mismatches ratio / % 0. 77 0. 78 0. 21 0. 73 0. 80 0. 86 0. 80 0. 85
Notes: TCTS coverage is at least 95% of each TCTS length matched the gene, and is computed by match length / TCTS length; gene coverage is at
least 80% of each gene length matched the TCTS, and is computed by match length / gene length. The assembly of each assembler was aligned to
15 222 D. melanogaster gene sequences using BLAT for TCTS coverage. For gene coverage, the same D. melanogaster gene sequences were aligned
to the assembly of each assembler using BLAT. Mismatches ratio is computed by total mismatched bases / total de novo assembled bases
Table 8摇 Summaries of O. sativa gene coverage, the percentages of TCTSs covering 95% of TCTS
matched gene sequences and mismatches ratio for each de novo assembly program
AB Ve SO Oa Tr Mk TI TA
% of TCTSs achieving TCTS coverage 37. 6 37. 2 20. 3 40. 0 56. 2 36. 4 29. 7 36. 7
Number of genes achieving total gene coverage 127 36 625 1 172 906 81 56 627
Mismatches ratio / % 0. 22 0. 28 0. 10 0. 21 0. 19 0. 30 0. 31 0. 25
Notes: TCTS coverage is at least 95% of each TCTS length matched the gene, and is computed by match length / TCTS length; gene coverage is at
least 80% of each gene length matched the TCTS, and is computed by match length / gene length. The assembly of each assembler was aligned to
56 797 O. sativa gene sequences using BLAT for TCTS coverage. For gene coverage, the same O. sativa gene sequences were aligned to the assem鄄
bly of each assembler using BLAT. Mismatches ratio is computed by total mismatched bases / total de novo assembled bases
294摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 植 物 分 类 与 资 源 学 报摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 第 34 卷
lengths produced by Oases and Trinity for both data鄄
sets indicated a more desirable TCTS length distribu鄄
tion ( i. e., more long TCTSs). Trans鄄ABySS yiel鄄
ded fewer and shorter TCTSs for Drosophila MS than
for cultivated rice, while SOAPdenovo showed the
opposite trend with fewer and shorter TCTSs for the
Drosophila MS dataset. Other summary statistics in鄄
cluded measurements of the total numbers of TCTSs
(categorized by size: 逸100 bp, 臆200 bp, and 逸
1 000 bp) as well as their maximum, mean, and
median lengths (Tables 3 and 4). Trans鄄ABySS as鄄
sembled the most total bases (approximately 70 MB
for Drosophila MS and 160 MB for cultivated rice),
although it also produced more of the shorter TCTSs
Fig. 1摇 Processing requirements in hours and RAM for each de novo
transcriptome assembler using Drosophila MS RNA鄄Seq reads
1. ABySS; 2. Velvet; 3. SOAPdenovo; 4. Oases; 5. Trinity;
6. Multiple鄄k; 7. T鄄IDBA; 8. Trans鄄ABySS
Fig. 2摇 Processing requirements in hours and RAM for each de novo
transcriptome assembler using cultivated rice RNA鄄Seq reads
1. ABySS; 2. Velvet; 3. SOAPdenovo; 4. Oases; 5. Trinity;
6. Multiple鄄k; 7. T鄄IDBA; 8. Trans鄄ABySS
3945 期摇 摇 摇 摇 ZHAO Lei et al. : Comparing De Novo Transcriptome Assemblers Using Illumina RNA鄄Seq Reads摇 摇 摇 摇 摇
with nearly 47% and 58% of total TCTSs for both
datasets (Figures 3 and 4). Trinity produced the fe鄄
west total TCTSs and shorter TCTSs; yet combined
TCTS lengths were 16 752 840 bp and 79 027 325 bp
for Drosophila MS and cultivated rice, respectively,
supporting the conclusion that the length of each
TCTS was much longer. The largest and second lar鄄
gest maximum TCTS sizes were produced by Oases
and SOAPdenovo, respectively, for both Drosophila
MS and cultivated rice (Tables 3 and 4). The lar鄄
gest mean and median TCTS sizes were produced by
Oases for Drosophila MS, but by Trinity for cultivat鄄
ed rice, and the second鄄largest mean and median si鄄
zes were produced by Trinity for Drosophila MS and
Oases for cultivated rice (Tables 3 and 4).
N50 and N90 have also been commonly de鄄
scribed summary statistics for RNA鄄Seq transcrip鄄
tome assembly, each defined as the TCTS length for
which 50% or 90% of the TCTSs, respectively, are
of the same length or larger ( Robertson et al.,
2010). They were determined by sorting all TCTSs
by length from longest to shortest, and then selecting
the length of the shortest TCTS such that the sum of
all TCTSs of equal length or longer was at least 50%
and 90% of the total length of all TCTSs. N50 and
N90 estimate the connectivity of the assembly and
can be used to compare estimated assembly quality,
with higher values indicating better performance.
For Drosophila MS, Oases produced the largest N50
and N90 values at 2 081 bp and 551 bp, respective鄄
ly, while Trinity came in a distant second with 943
bp and 379 bp (Table 3). The same N50 trend was
detected for cultivated rice transcriptome assembly,
with Oases scoring best with 1 400 bp and Trinity
scoring second鄄best with 1 247 bp (Table 4). The
biggest N90 TCTS for the cultivated rice dataset,
however, was produced by Trinity with 410 bp, while
Oases was second biggest with 293 bp (Table 4).
494摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 植 物 分 类 与 资 源 学 报摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 第 34 卷
Taken together, the contiguity鄄based summary
statistics of these eight assemblers showed that Oa鄄
ses, Trinity, Trans鄄ABySS and SOAPdenovo per鄄
formed better than the others, although other factors
besides TCTS contiguity should also be considered.
Measured by total numbers of bases, the overall TCTS
yields produced by the multiple k鄄mer assemblers
(Multiple鄄k, T鄄IDBA and Trans鄄ABySS) ranged from
20% to 60% larger than the yields from the single k鄄
mer assemblers (ABySS, Velvet, SOAPdenovo, Oa鄄
ses, Trinity). This suggested that the assembly in鄄
tegrity produced by the multiple k鄄mer assemblers
was higher, i. e., these assemblers avoided more
TCTS losses than the single k鄄mer assemblers.
Mapping of reads
Assembler quality depends on the maximal usage
of reads to unambiguously map the TCTSs. There鄄
fore, to compare all assemblers, each of the paired鄄
end reads from the Drosophila MS and cultivated rice
datasets was mapped separately back to each assem鄄
bled transcriptome using Bowtie software ( http: / /
bowtie鄄bio. sourceforge. net / index. shtml ) ( Lang鄄
mead et al., 2009). The Bowtie mapping algorithm
was used to map short reads with a seed length of
28, a maximum of two allowed mismatches in the
seed, and the maximum permitted total quality val鄄
ues of 70 at all mismatched read positions throughout
the entire alignment ( cf. Bowtie instruction manu鄄
al) . The assembly with the highest percentage of
back鄄mapped reads was produced by Trinity for both
datasets (78. 6% and 75. 4% in Drosophila MS and
cultivated rice, respectively ), followed by Oases
(56. 6% and 59. 8% ) and Trans鄄ABySS (64. 3%
and 54. 4% ) (Tables 3 and 4). The fewest reads
were mapped to the Velvet assembled transcriptome,
while ABySS, SOAPdenovo, Multiple鄄k and T鄄IDBA
fared comparably well in terms of the ratio of mapped
reads. The assemblies of Trinity, ABySS, Velvet,
Multiple鄄k and T鄄IDBA showed no reads with multi鄄
ple matches (逸5), suggesting that they generated
the least redundancy of all assemblers. Conversely,
Trans鄄ABySS, Oases and SOAPdenovo seemed to
produce high assembly redundancy.
Accuracy, integrity and coverage
Transcriptome assemblies have been described
in terms of their accuracy, integrity, and coverage
with respect to the actual transcriptome being ana鄄
lyzed. In this study, accuracy was measured by
comparing each TCTS with all known coding se鄄
quences from D. melanogaster or O. sativa genome
closely related Drosophila MS and cultivated rice
species, respectively, and then determining the per鄄
centage of all TCTSs which overlapped with known
coding sequences. Integrity was defined as the total
percentage of all known proteins for which TCTSs
were present, i. e., the degree of proteome represen鄄
tation by the complete TCTS population. The optimal
assembler in terms of transcriptome assembly quality
would be the program that produces the highest com鄄
bined accuracy and integrity of the transcriptome.
After de novo assembling the transcriptomes from the
Drosophila MS and cultivated rice reads, we used
BLASTX (ftp: / / ftp. ncbi. nlm. nih. gov / blast / execut鄄
ables / release / 2. 2. 25 / blast鄄2. 2. 25鄄x64鄄linux. tar. gz)
to compare them to 23 711 D. melanogaster protein se鄄
quences (ftp: / / ftp. flybase.net / genomes / Drosophila_
melanogaster / dmel_ r5. 38 _FB2011 _06 / fasta / dmel鄄
all鄄translation鄄r5. 38. fasta. gz) and 67 393 O. sativa
protein sequences (ftp: / / ftp. plantbiology.msu. edu /
pub / data / Eukaryotic_Projects / o_sativa / annotation_
dbs / pseudomolecules / version_6. 1 / all. dir / all. pep)
to calculate assembler accuracy. To compute integri鄄
ty, we used TBLASTN to compare the same D. mela鄄
nogaster and O. sativa protein sequences to the de no鄄
vo transcriptomes assembled by each assembler (Alts鄄
chul et al., 1990; Mount, 2007). For the alignment
of each query sequence, only the smallest E value hit
(cutoff =1e鄄5) was considered. Trinity performed best
in terms of assembly accuracy (92. 5% and 82. 1%
for Drosophila MS and cultivated rice, respectively),
and Trans鄄ABySS (81. 6% and 68. 5%) and Oases
(76.5% and 62.9%) also fared well (Tables 5 and 6).
Trinity also obtained the highest integrity (87. 6%)
on the assembly of sequence reads from cultivated
5945 期摇 摇 摇 摇 ZHAO Lei et al. : Comparing De Novo Transcriptome Assemblers Using Illumina RNA鄄Seq Reads摇 摇 摇 摇 摇
rice, followed by SOAPdenovo ( 84. 5%), Oases
(84. 0% ) and Trans鄄ABySS (83. 5% ) (Table 6).
Assembly of the Drosophila MS reads led to different
results for integrity, with Oases first at 87. 1% and
Multiple鄄k second at 86. 4% , while Trinity was to鄄
ward the bottom (76. 9% ) (Table 5).
For further analysis of the quality of de novo
transcriptome assemblies produced by the eight test鄄
ed programs, the TCTS coverage and gene coverage
were calculated (Tables 7 and 8). The UCSC Blat
(Kent, 2002) software ( http: / / genome鄄test. cse.
ucsc. edu / ~ kent / exe / linux / blatSuite. 34. zip ) was
used to align all TCTSs to 15 222 D. melanogaster
gene sequences ( ftp: / / ftp. flybase. net / genomes /
Drosophila_melanogaster / dmel_r5. 38_FB2011_06 /
fasta / dmel鄄all鄄gene鄄r5. 38. fasta. gz) and 56 797 O.
sativa gene sequences ( ftp: / / ftp. plantbiology. msu.
edu / pub / data / Eukaryotic_Projects / o_sativa / annota鄄
tion_dbs / pseudomolecules / version_6. 1 / all. dir / all.
seq) to determine TCTS coverage. Meanwhile, the
same D. melanogaster and O. sativa gene sequences
were aligned to the actual transcriptomes to deter鄄
mine gene coverage. For the alignment of each query
sequence, only the highest鄄scoring hit ( see Materi鄄
als and Methods ) was used. Trinity showed the
highest percentage of all TCTSs that covered at least
95% of each TCTS (72. 9% and 56. 2% for Dro鄄
sophila MS and cultivated rice, respectively), fol鄄
lowed by Trans鄄ABySS (55. 8% and 36. 7% ) and
Oases (40. 0% for both) (Tables 7 and 8). These
results resembled those obtained by using BLASTX,
supporting the conclusion that the de novo assembly
of transcriptome obtained from Trinity possessed
the highest accuracy. The program yielding TCTSs
covering at least 80% of the most genes was Oases
(1 097 and 1 172), followed by Trinity (549 and
906), Trans鄄ABySS (146 and 627) and SOAPde鄄
novo ( 127 and 625 ), indicating these programs
have better capacity to reconstruct more full鄄length
transcripts (Tables 7 and 8).
Mismatches ratios and alternative splicing
Mismatches ratios were compared as another
way to assess the performance of these de novo as鄄
semblers. These data were obtained by aligning each
de novo assembly transcriptome to the protein and
gene databases and using BLASTX (Tables 5 and 6)
and Blat (Tables 7 and 8) (cf. Materials and Meth鄄
ods). Oases showed the highest mismatches rate
(0. 53% ), and ABySS displayed the lowest mismat鄄
ches rate (0. 30% ) for the Drosophila MS transcrip鄄
tome ( Table 5), while Trinity showed the highest
mismatches rate (0. 78% ), and SOAPdenovo dis鄄
played the lowest mismatches rate (0. 51% ) for cul鄄
tivated rice (Table 6). The mismatches rates pro鄄
duced by the other assemblers were comparably low
to the least effective for each organism. Additional鄄
ly, when Blat aligner was used, Multiple鄄k showed
the highest mismatches rate (0. 86% ), and SOAPd鄄
enovo displayed the lowest mismatches rate (0. 21%)
for Drosophila MS (Table 7), while T鄄IDBA showed
the highest (0. 31% ), and SOAPdenovo displayed
the lowest mismatches rate (0. 10% ) for cultivated
rice (Table 8).
Alternative splicing is the one of the biggest
challenges in RNA鄄Seq de novo assembly ( Martin
and Wang, 2011; Wang et al., 2010a). As such,
we compared the efficiency of alternative鄄splicing
resolution among these de novo assemblers. Each as鄄
sembled transcriptome was mapped separately to D.
melanogaster ( ftp: / / ftp. flybase. net / genomes / Dro鄄
sophila_melanogaster / dmel_r5. 38_FB2011_06 / fas鄄
ta / dmel鄄all鄄chromosome鄄r5. 38. fasta. gz) and O. sa鄄
tiva ( ftp: / / ftp. ensemblgenomes. org / pub / plants /
release鄄11 / fasta / oryza _ sativa / dna / Oryza _ sativa.
MSU6. dna. toplevel. fa. gz) genome sequences to
detect novel and known splices using GMAP (Wu
and Watanabe, 2005) software ( cf. Materials and
Methods). We compared the numbers of novel and
annotated splice sites captured by each assembler
(Fig. 5). In Drosophila MS, Trans鄄ABySS identi鄄
fied the most known splices (73 320), followed by
Oases (44 810) and Multiple鄄k (34 430 ); then
Trans鄄ABySS identified the largest number of novel
splices (32 876), followed by Oases (24 162) and
694摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 植 物 分 类 与 资 源 学 报摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 第 34 卷
Multiple鄄k (18 645). In cultivated rice, Trans鄄ABySS
still identified the largest number of known splices
(265 228), followed by Oases (136 884), Trinity
(132 213 ); then Trans鄄ABySS still the largest
number of novel splices (122 613), followed Oases
(77 050), and Trinity (75 170). In terms of process鄄
ing alternatively spliced transcripts, Trans鄄ABySS,
Oases and Trinity showed the best performance; oth鄄
ers were closely comparable (Fig. 5).
Combining assembler programs into a novel workflow
Despite the burgeoning number of de novo
RNA鄄Seq assemblers, transcriptome assembly re鄄
mains challenging (Wang et al., 2010a; Nagarajan
and Pop, 2009; Ozsolak and Milos, 2011; Martin
and Wang, 2011). Many complicating factors are
regularly encountered, such as short read lengths,
poor read quality, low sequencing depth, and alter鄄
natively spliced transcripts ( Paszkiewicz and Stud鄄
holme, 2010; Wall et al., 2009). As such, the us鄄
er鄄defined single k鄄mer programs are unlikely to
yield optimal performance for de novo transcriptome
assembly, but even the multiple k鄄mer programs can
fall short of satisfactory results (Garber et al., 2011;
Zerbino, 2010). Therefore, a strategy employing a
combination of multiple k鄄mer programs should be
produced to increase assembly success (Surget鄄Gro鄄
ba and Montoya鄄Burgos, 2010; Robertson et al.,
2010; Martin et al., 2010). To take full advantage
of the assembly capabilities of different multiple k鄄
mer assemblers, we have designed the “ ETM 冶
method. It consists of three components, beginning
with the correction of read errors, followed by the
production of a series of k鄄mer Trinity assemblies,
and concluding with a merging step and the removal
of redundancy (Fig. 6).
To maximize the quality of the short sequence
reads, the ETM process begins by applying the ECH鄄
O algorithm to correct read errors. ECHO uses a max鄄
imum a posteriori model on read overlaps, assigning
quality scores to the bases to avoid false overlaps
(Kao et al., 2011). ETM then employs the Trinity
assembler on the corrected reads, combining differ鄄
ent series of k鄄mers. Finally, these incremental as鄄
semblies are merged into the final assembly using
Multiple鄄k爷 s CD鄄HIT鄄EST function ( http: / / www.
bioinformatics. org / cd鄄hit / ) (Li and Godzik, 2006)
Fig. 5摇 The results of the alternative
splicing analysis
Fig. 6摇 Summary of the ETM摇
assembly workflow摇
7945 期摇 摇 摇 摇 ZHAO Lei et al. : Comparing De Novo Transcriptome Assemblers Using Illumina RNA鄄Seq Reads摇 摇 摇 摇 摇
to remove redundancy and retain the best possible
assembly (see Additional file 5, http: / / journal.kib.
ac.cn / UserFiles / File / 2012鄄084_editing(1). rar). In
this study, four ETM assemblies from each of the
two datasets were produced using the incremental
Trinity assemblies as input (Table S1, http: / / jour鄄
nal.kib. ac. cn / UserFiles / File / 2012鄄084_editing(1).
rar), and then the best of the four assemblies was
chosen. We found that our ETM assemblies showed
better TCTS summary statistics than those resulting
from any of the eight assembler programs, including
higher values for total bases, N50, N90, and mean
and median lengths, indicating a substantial im鄄
provement in contiguity. Surprisingly, the multi鄄hit
reads (逸5) remained absent from the ETM results,
indicating that little redundancy was produced by
this method. Accuracy, integrity, coverage and mis鄄
matches rates were also superior in the ETM products
(Table 9). These results suggested that our strategy
of combining multiple k鄄mer programs in a serial
workflow was advantageous for de novo transcriptome
assembly on non鄄model organisms.
Discussion and conclusions
We compared the de novo assembly performance
of eight assembly programs, using two datasets pro鄄
duced by the Illumina sequencing platform. For
each de novo assembly program, we varied the k鄄mer
parameters ( see Additional File 1, Table S1, ht鄄
tp: / / journal. kib. ac. cn / UserFiles / File / 2012鄄084 _
editing(1). rar) and chose the best assembly results
for the comparisons. Trinity showed the best per鄄
formance of the single k鄄mer tools, followed by Oa鄄
ses and SOAPdenovo. Trans鄄ABySS performed best
of the multiple k鄄mer programs, although the others
were comparably effective. Additionally, we designed
the ETM method, and found it to be well鄄suited for
de novo transcriptome assembly.
The quality of de novo assembly transcriptome
can be estimated by using the information of closely
related species (Kumar and Blaxter, 2010; Wang et
al., 2010c; Ness et al., 2011; Feldmeyer et al.,
Table 9摇 Analysis of ETM assemblies for Drosophila MS
and cultivated rice transcriptomes
Drosophila MS Cultivated rice
TCTS逸100 bp 45 424 138 637
TCTS臆200 bp 0 0
TCTS逸1 kb 10 723 46 276
Max 8 519 13 038
Mean 796 969
Median 588 670
N50 972 1 344
N90 393 429
% of reads mapped 80. 2 77. 0
Multi鄄hit reads (逸5) 0 0
% of TCTSs hit 93. 5 84. 5
% of TCTSs sharing 80%
identity to a known protein 92. 9 74. 7
% of proteins hit 78. 2 89. 2
% of proteins showing 80%
identity with a TCTS 58. 9 57. 5
Mismatches ratio (% )
using BLASTX 0. 34 0. 77
% of TCTSs with full TCTS
coverage 73. 9 58. 7
Number of genes with full
gene coverage 571 1 173
Mismatches ratio (% )
Using Blat 0. 84 0. 18
Novel Splices 21 427 39 893
Know Splices 123 888 216 506
Notes: E鄄value cutoff =1e鄄5. TCTS coverage represents the percentage
by which TCTS length exactly matched the gene coding sequence,
computed by match length / TCTS length; gene coverage is at least
80% of each gene length matched the TCTS, and is computed by
match length / gene length. Mismatches ratio is computed by total mis鄄
matched bases / total de novo assembled bases
2011). To evaluate and compare the performance of
the programs as well as our ETM workflow, we took
advantage of D. melanogaster or O. sativa for refer鄄
ence. Due to no true genomes for Drosophila MS and
cultivated rice species, and exist the divergence of
sequence and structure to D. melanogaster and O. sa鄄
tiva from the aspect of evolution, so, genome鄄guided
assemblers (Cufflinks and Scripture) were not used
for getting the TCTSs of Drosophila MS and cultivated
rice; otherwise, if done, wrong transcripts may be
produced, then the assessment may also not be accu鄄
rate. Therefore, genome鄄guided assemblers ( Cuf鄄
flinks and Scripture) were not compared with de novo
894摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 植 物 分 类 与 资 源 学 报摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 第 34 卷
assemblers. Since this study was designed to compare
de novo transcriptome assembly programs using mR鄄
NA sequencing reads, D. melanogaster and O. sativa
protein and gene datasets can provide excellent refer鄄
ence sequences. To support this conclusion, each
transcriptome de novo assembled was also aligned to
D. melanogaster and O. sativa genome sequences, and
the results were shown (see Additional File 1, Tables
S2鄄S5, http: / / journal. kib. ac. cn / UserFiles / File /
2012鄄084_editing(1). rar). As we known, sequence
depth has various effects on the accuracy of de novo
assembly (Garber et al., 2011; Martin and Wang,
2011). The main purposes of this study were to com鄄
pare the efficiency of each de novo assembler based
the same trial datasets, so sequence depth was not
deeply explored. Generally speaking, sequence depth
should achieve at least 30 伊 to obtain a high quality
de novo transcriptome assembly (Martin and Wang,
2011). When the reads were mapped back to the as鄄
sembly, the threshold value for multiple hit reads
(逸5) allowed estimation of the possible redundancy
of the assemblies, since many sequence reads likely
derived from transcripts of gene family members.
Here, we performed a systematic comparison
study to evaluate de novo assembly programs accord鄄
ing to a series of summary statistics, including conti鄄
guity, accuracy, integrity, and coverage and mis鄄
matches rate. Using our results, researchers may
choose the appropriate assembler based on their spe鄄
cific experimental data and the computational re鄄
sources available. Many factors may affect de novo
transcriptome assembly, however, such as the pres鄄
ence of gene duplications and alternatively spliced
transcripts. To accommodate the variability in tran鄄
script abundance intrinsic to gene expression, and to
improve on assembly quality, we recommend a serial
multiple k鄄mer strategy such as our ETM method to
obtain maximal transcriptome coverage. Such a strat鄄
egy should also be considered when designing future
assembly programs, since our ETM method performed
better transcriptome assembly than any single program
currently available. In addition, each NGS platform
has certain advantages and disadvantages, and assem鄄
bly performance could be improved by combining the
best complementary advantages of the different NGS
technologies and de novo assembly tools (Adamidi et
al., 2011; Iorizzo et al., 2011; Su et al., 2011;
Angeloni et al., 2011; Sandmann et al., 2011).
Acknowledgements: We thank the assistance of Liu Qi, In鄄
stitute of Genomic Medicine / Zhejiang Provincial Key Labora鄄
tory of Medical Genetics, Zhuang Hui鄄Fu, Kunming Institute
of Botany, Chinese Academy of Sciences.
References:
Adamidi C, Wang Y, Gruen D et al., 2011. De novo assembly and
validation of planaria transcriptome by massive parallel sequen鄄
cing and shotgun proteomics [ J] . Genome Research, 21 (7):
1193—1200
Altschul SF, Gish W, Miller W et al., 1990. Basic local alignment
search tool [J] . Journal of Molecular Biology, 215 (3): 403—
410
Angeloni F, Wagemaker CA, Jetten MS et al., 2011. De novo tran鄄
scriptome characterization and development of genomic tools for
Scabiosa columbaria L. using next鄄generation sequencing tech鄄
niques [J] . Molecular Ecology Resources, 11 (4): 662—674
Ansorge WJ, 2009. Next鄄generation DNA sequencing techniques [J].
New Biotechnology, 25 (4): 195—203
Atambayeva SA, Khailenko VA, Ivashchenko AT, 2008. Intron and
exon length variation in Arabidopsis, rice, nematode, and human
[J] . Mathematical and Systems Biology, 42 (2): 312—320
Birol I, Jackman SD, Nielsen CB et al., 2009. De novo transcriptome
assembly with ABySS [J]. Bioinformatics, 25 (21): 2872—2877
Birzele F, Schaub J, Rust W et al., 2010. Into the unknown: expres鄄
sion profiling without genome sequence information in CHO by
next generation sequencing [ J ] . Nucleic Acids Research, 38
(12): 3999—4010
Chen S, Yang P, Jiang F et al., 2010. De novo analysis of transcrip鄄
tome dynamics in the migratory locust during the development of
phase traits [J] . PloS One, 5 (12): e15633
Collins LJ, Biggs PJ, Voelckel C et al., 2008. An approach to tran鄄
scriptome analysis of non鄄model organisms using short鄄read se鄄
quences [J] . Genome Inform, 21: 3—14
Crawford JE, Guelbeogo WM, Sanou A et al., 2010. De novo tran鄄
scriptome sequencing in Anopheles funestus using Illumina RNA鄄
seq technology [J] . PloS One, 5 (12): e14202
Deutsch M, Long M, 1999. Intron鄄exon structures of eukaryotic model
organisms [J]. Nucleic Acids Research, 27 (15): 3219—3228
Feldmeyer B, Wheat CW, Krezdorn N et al., 2011. Short read Illu鄄
mina data for the de novo assembly of a non鄄model snail species
9945 期摇 摇 摇 摇 ZHAO Lei et al. : Comparing De Novo Transcriptome Assemblers Using Illumina RNA鄄Seq Reads摇 摇 摇 摇 摇
transcriptome (Radix balthica, Basommatophora, Pulmonata),
and a comparison of assembler performance [ J] . BMC Genom鄄
ics, 12: 317. 10. 1186 / 1471鄄2164鄄12鄄317
Fullwood MJ, Wei CL, Liu ET et al., 2009. Next鄄generation DNA
sequencing of paired鄄end tags ( PET) for transcriptome and ge鄄
nome analyses [J] . Genome Research, 19 (4): 521—532
Garber M, Grabherr MG, Guttman M et al., 2011. Computational
methods for transcriptome annotation and quantification using
RNA鄄seq [J] . Nature Methods, 8 (6): 469—477
Garcia TI, Shen Y, Catchen J et al., 2012. Effects of short read
quality and quantity on a de novo vertebrate transcriptome assem鄄
bly [ J ] . Comparative Biochemistry and Physiology鄄Part C:
Taxicology & Pharmacology, 155 (1): 95—101
Garg R, Patel RK, Tyagi AK et al., 2011. De novo assembly of
chickpea transcriptome using short reads for gene discovery and
marker identification [J] . DNA Research, 18 (1): 53—63
Ghiselli F, Milani L, Chang PL et al., 2012. De novo assembly of
the manila clam Ruditapes philippinarum transcriptome provides
new insights into expression bias, mitochondrial doubly uniparen鄄
tal inheritance and sex determination [ J] . Molecular Biology
and Evolution, 29 (2): 771—786
Gibbons JG, Janson EM, Hittinger CT et al., 2009. Benchmarking
next鄄generation transcriptome sequencing for functional and evo鄄
lutionary genomics [ J] . Molecular Biology and Evolution, 26
(12): 2731—2744
Grabherr MG, Haas BJ, Yassour M et al., 2011. Full鄄length tran鄄
scriptome assembly from RNA鄄Seq data without a reference ge鄄
nome [J] . Nature Biotechnology, 29 (7): 644—652
Guttman M, Garber M, Levin JZ et al., 2010. Ab initio reconstruc鄄
tion of cell type鄄specific transcriptomes in mouse reveals the con鄄
served multi鄄exonic structure of lincRNAs [ J] . Nature Biotech鄄
nology, 28 (5): 503—510
Haas BJ, Zody MC, 2010. Advancing RNA鄄Seq analysis [ J] . Na鄄
ture Biotechnology, 28 (5): 421—423
Hamilton JP, Hansey CN, Whitty BR et al., 2011. Single nucleotide
polymorphism discovery in elite north american potato germplasm
[J] . BMC Genomics, 12: 302
Hao da C, Ge G, Xiao P et al., 2011. The first Insight into the tissue
specific Taxus transcriptome via illumina second generation se鄄
quencing [J] . PloS One, 6 (6): e21220
Hu H, Bandyopadhyay PK, Olivera BM et al., 2011. Characteriza鄄
tion of the Conus bullatus genome and its venom鄄duct transcrip鄄
tome [J] . BMC Genomics, 12: 60
Iorizzo M, Senalik DA, Grzebelus D et al., 2011. De novo assembly
and characterization of the carrot transcriptome reveals novel
genes, new markers, and genetic diversity [ J] . BMC Genomics,
12 (1): 389
Jager M, Ott CE, Grunhagen J et al., 2011. Composite transcriptome
assembly of RNA鄄seq data in a sheep model for delayed bone
healing [J] . BMC Genomics, 12: 158
Kao WC, Chan AH, Song YS, 2011. ECHO: A reference鄄free short鄄
read error correction algorithm [J] . Genome Research, 21 (7):
1181—1192
Kawahara鄄Miki R, Wada K, Azuma N et al., 2011. Expression profi鄄
ling without genome sequence information in a non鄄model spe鄄
cies, pandalid shrimp (Pandalus latirostris), by next鄄generation
sequencing [J] . PloS One, 6 (10): e26043
Kent WJ, 2002. BLAT—the BLAST鄄like alignment tool [ J] . Ge鄄
nome Research, 12 (4): 656—664
Kumar S, Blaxter ML, 2010. Comparing de novo assemblers for 454
transcriptome data [J] . BMC Genomics, 11: 571
Langmead B, Trapnell C, Pop M et al., 2009. Ultrafast and memory鄄
efficient alignment of short DNA sequences to the human genome
[J] . Genome Biology, 10 (3): R25
Li R, Yu C, Li Y et al., 2009. SOAP2: an improved ultrafast tool for
short read alignment [J]. Bioinformatics, 25 (15): 1966—1967
Li W, Godzik A, 2006. Cd鄄hit: a fast program for clustering and
comparing large sets of protein or nucleotide sequences [ J] .
Bioinformatics, 22 (13): 1658—1659
Martin J, Bruno VM, Fang Z et al., 2010. Rnnotator: an automated
de novo transcriptome assembly pipeline from stranded RNA鄄Seq
reads [J] . BMC Genomics, 11: 663
Martin JA, Wang Z, 2011. Next鄄generation transcriptome assembly
[J] . Nature Reviews Genetics, 12 (10): 671—682
McManus CJ, Coolon JD, Duff MO et al., 2010. Regulatory diver鄄
gence in Drosophila revealed by mRNA鄄seq [ J] . Genome Re鄄
search, 20 (6): 816—825
Metzker ML, 2010. Sequencing technologies 鄄 the next generation
[J] . Nature Reviews Genetics, 11 (1): 31—46
Miller JR, Koren S, Sutton G, 2010. Assembly algorithms for next鄄
generation sequencing data [J] . Genomics, 95 (6): 315—327
Mizrachi E, Hefer CA, Ranik M et al., 2010. De novo assembled ex鄄
pressed gene catalog of a fast鄄growing Eucalyptus tree produced by
Illumina mRNA鄄Seq [J] . BMC Genomics, 11: 681
Mortazavi A, Williams BA, McCue K et al., 2008. Mapping and
quantifying mammalian transcriptomes by RNA鄄Seq [ J] . Nature
Methods, 5 (7): 621—628
Mount DW, 2007. Using the Basic Local Alignment Search Tool
(BLAST) [Z]. CSH Protoc 2007: pdb top17
Nagalakshmi U, Wang Z, Waern K et al., 2008. The transcriptional
landscape of the yeast genome defined by RNA sequencing [J] .
Science, 320 (5881): 1344—1349
Nagarajan N, Pop M, 2009. Parametric complexity of sequence as鄄
sembly: theory and applications to next generation sequencing
[J] . Journal of Computational Biology, 16 (7): 897—908
Ness RW, Siol M, Barrett S, 2011. De novo sequence assembly and
characterization of the floral transcriptome in cross鄄 and self鄄fer鄄
tilizing plants [J] . BMC Genomics, 12: 298
Ozsolak F, Milos PM, 2011. RNA sequencing: advances, challenges
and opportunities [J]. Nature Reviews Genetics, 12 (2): 87—98
005摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 植 物 分 类 与 资 源 学 报摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 第 34 卷
Paszkiewicz K, Studholme DJ, 2010. De novo assembly of short se鄄
quence reads [J] . Brief Bioinform, 11 (5): 457—472
Peng Y, Leung HCM, Yiu SM et al., 2011. T鄄IDBA: a de novo iter鄄
ative de bruijn graph assembler for transcriptome [J] . Research
in Computational Molecular Biology, 6577: 337—338
Qin YF, Fang HM, Tian QN et al., 2011. Transcriptome profiling
and digital gene expression by deep鄄sequencing in normal / regen鄄
erative tissues of planarian Dugesia japonica [ J] . Genomics, 97
(6): 364—371
Qiu Q, Ma T, Hu Q et al., 2011. Genome鄄scale transcriptome analy鄄
sis of the desert poplar, Populus euphratica [ J] . Tree Physiolo鄄
gy, 31 (4): 452—461
Robertson G, Schein J, Chiu R et al., 2010. De novo assembly and a鄄
nalysis of RNA鄄seq data [J]. Nat Methods, 7 (11): 909—912
Sandmann T, Vogg MC, Owlarn S et al., 2011. The head鄄regenera鄄
tion transcriptome of the planarian Schmidtea mediterranea [ J] .
Genome Biology, 12: R76
Santure AW, Gratten J, Mossman JA et al., 2011. Characterisation of
the transcriptome of a wild great tit Parus major population by next
generation sequencing [J] . BMC Genomics, 12: 283
Shendure J, Ji H, 2008. Next鄄generation DNA sequencing [J] . Nat
Biotechnol, 26 (10): 1135—1145
Shi CY, Yang H, Wei CL et al., 2011. Deep sequencing of the Ca鄄
mellia sinensis transcriptome revealed candidate genes for major
metabolic pathways of tea鄄specific compounds [J] . BMC Genom鄄
ics, 12: 131
Simon SA, Zhai J, Nandety RS et al., 2009. Short鄄read sequencing
technologies for transcriptional analyses [ J] . Annual Review of
Plant Biology, 60: 305—333
Simpson JT, Wong K, Jackman SD et al., 2009. ABySS: a parallel
assembler for short read sequence data [ J] . Genome Research,
19 (6): 1117—1123
Su CL, Chao YT, Alex Chang YC et al., 2011. De novo assembly of
expressed transcripts and global analysis of Phalaenopsis aphrodi鄄
te transcriptome [ J] . Plant and Cell Physiology, 52 ( 9 ):
1501—1514
Surget鄄Groba Y, Montoya鄄Burgos JI, 2010. Optimization of de novo
transcriptome assembly from next鄄generation sequencing data
[J] . Genome Research, 20 (10): 1432—1440
Tang Q, Ma XJ, Mo CM et al., 2011. An efficient approach to find鄄
ing Siraitia grosvenorii triterpene biosynthetic genes by RNA鄄seq
and digital gene expression analysis [ J] . BMC Genomics, 12
(1): 343
Trapnell C, Williams BA, Pertea G et al., 2010. Transcript assembly
and quantification by RNA鄄Seq reveals unannotated transcripts
and isoform switching during cell differentiation [J] . Nature Bi鄄
otechnology, 28 (5): 511—515
Wall CE, Cozza S, Riquelme CA et al., 2011. Whole transcriptome
analysis of the fasting and fed Burmese python heart: insights in鄄
to extreme physiological cardiac adaptation [ J] . Physiological
Genomics, 43 (2): 69—76
Wall PK, Leebens鄄Mack J, Chanderbali AS et al., 2009. Comparison
of next generation sequencing technologies for transcriptome char鄄
acterization [J] . BMC Genomics, 10: 347
Wang L, Li P, Brutnell TP, 2010a. Exploring plant transcriptomes u鄄
sing ultra high鄄throughput sequencing [J] . Brief Funct Genom鄄
ics, 9 (2): 118—128
Wang XW, Luan JB, Li JM et al., 2010b. De novo characterization
of a whitefly transcriptome and analysis of its gene expression
during development [J] . BMC Genomics, 11: 400
Wang Z, Fang B, Chen J et al., 2010c. De novo assembly and char鄄
acterization of root transcriptome using Illumina paired鄄end se鄄
quencing and development of cSSR markers in sweet potato (Ip鄄
omoea batatas) [J] . BMC Genomics, 11: 726
Wang Z, Gerstein M, Snyder M, 2009. RNA鄄Seq: a revolutionary
tool for transcriptomics [J] . Nature Reviews Genetics, 10 (1):
57—63
Wenping H, Yuan Z, Jie S et al., 2011. De novo transcriptome se鄄
quencing in Salvia miltiorrhiza to identify genes involved in the
biosynthesis of active ingredients [J]. Genomics, 98: 272—279
Wickramasinghe S, Hua S, Rincon G et al., 2011. Transcriptome
profiling of bovine milk oligosaccharide metabolism genes using
RNA鄄sequencing [J] . PLoS One, 6 (4): e18895
Wong MM, Cannon CH, Wickneswari R, 2011. Identification of lignin
genes and regulatory sequences involved in secondary cell wall
formation in Acacia auriculiformis and Acacia mangium via de no鄄
vo transcriptome sequencing [J]. BMC Genomics, 12 (1): 342
Wu TD, Watanabe CK, 2005. GMAP: a genomic mapping and align鄄
ment program for mRNA and EST sequences [J] . Bioinformat鄄
ics, 21 (9): 1859—1875
Xiang LX, He D, Dong WR et al., 2010. Deep sequencing鄄based
transcriptome profiling analysis of bacteria鄄challenged Lateolabrax
japonicus reveals insight into the immune鄄relevant genes in marine
fish [J] . BMC Genomics, 11: 472
Xue J, Bao YY, Li BL et al., 2010. Transcriptome analysis of the
brown planthopper Nilaparvata lugens [J] . PloS One, 5 (12):
e14233
Yang SS, Tu ZJ, Cheung F et al., 2011. Using RNA鄄Seq for gene i鄄
dentification, polymorphism detection and transcript profiling in
two alfalfa genotypes with divergent cell wall composition in
stems [J] . BMC Genomics, 12: 199
Zerbino DR, 2010. Using the Velvet de novo assembler for short鄄read
sequencing technologies [ J] . Current Protocols in Bioinformat鄄
ics, 31: 11. 15. 11—11. 15. 12
Zerbino DR, Birney E, 2008. Velvet: algorithms for de novo short
read assembly using de Bruijn graphs [J] . Genome Research, 18
(5): 821—829
Zhang G, Guo G, Hu X et al., 2010. Deep RNA sequencing at sin鄄
gle base鄄pair resolution reveals high complexity of the rice tran鄄
scriptome [J] . Genome Research, 20 (5): 646—654
1055 期摇 摇 摇 摇 ZHAO Lei et al. : Comparing De Novo Transcriptome Assemblers Using Illumina RNA鄄Seq Reads摇 摇 摇 摇 摇