Optimizing and benchmarking de novo transcriptome sequencing: from library preparation to assembly evaluation

Large-size library inserts improve connectivity of transcriptome assembly

We introduced two major changes to the standard protocol of Illumina TruSeq RNA Sample Prep Kit in regards to RNA fragmentation and size selection of inserts. First, we shortened RNA fragmentation time to increase overall insert lengths. Second, in DNA purification, we decreased the volume of Agencourt AMPure XP, aiming to retain libraries with large-size fragments, e.g., those longer than 300 bp. Using total RNA extracted from a gecko embryo four days before the estimated date of oviposition [−4 days post oviposition (dpo)], we prepared libraries without and with the above modifications to the standard protocol (Library A and Library B, respectively; Table 1). For Library A and B, we loaded the same number of DNA molecules, based on prior quantification, in the same numbers of lanes (Table 1), and confirmed that it resulted in comparable numbers of reads (Additional file 1). Mapping of the reads to the de novo assemblies (see Methods) indicated that Library B had larger overall fragment sizes (665 bp on average) and a broad size distribution compared to Library A (349 bp on average; Fig. 2). Additionally, the fraction corresponding to reads shorter than 300 bp was largely reduced in the size distribution of Library B.

Library	Library preparation				Sequencing
	RNA source	Duration of RNA fragmentation (min)	× AMPure volume (targeted fraction to retain)	PCR cycles
A	-4 dpo whole embryoa	8	× 1.6 (>100 bp)	6	HiSeq 1/4 lanes, 171 cycles paired-end
B		2	× 0.7 (>300 bp)	6
C				6	MiSeq 1/4 lanes, 250 cycles paired-end
D	9 dpo whole embryo	8	× 1.6 (>100 bp)	2	HiSeq 1/4 lanes, 171 cycles paired-end
E		2	× 0.7 (>300 bp)	4
F	30 dpo head	4	× 1.0 (>150 bp)	6	HiSeq 2/3 lanes, 151 cycles paired-end
G	30 dpo liver			6
H	30 dpo tail			6

^aEmbryo of 4 days before the estimated day of oviposition

In order to evaluate the performance of the protocol modifications, we sequenced the libraries on the Illumina HiSeq with 171 cycles [23]. It is expected that longer reads should increase the proportion of large-size fragments, for which paired reads will cover the full stretch. We inferred the size of the sequenced inserts by mapping the paired reads to transcript contigs that were assembled from the obtained reads (see Methods). The result showed that the sizes of the sequenced fragments were smaller than the sizes of the majority of the fragments in the prepared library (Fig. 2). Nevertheless, the average size of sequenced fragments of Library B (360 bp) was larger than that of Library A (180 bp) (Fig. 2), which was consistent with the observations of the average sizes of the fragments in the prepared library.

Derivation of new reference gene set for vertebrates

Coverage of the protein-coding landscape, which we call ‘completeness’ in this article, is one of the typical quality measures of de novo transcriptome assembly [13, 24]. The pipeline CEGMA has been used for this purpose in eukaryotes [16, 17]. A subset of 458 CEGs possessing no or minimal paralogs (248 CEGs) is used as default in completeness assessment by CEGMA [17]. Each of the CEGs consists of only one gene per species, even if the species has its paralog due to a lineage-specific duplication. To examine potential effects of such additional paralogs, we analyzed molecular phylogenies of all the 248 CEGs based on the gene trees provided in Ensembl [25]. Overall, 71 duplication events in 64 CEGs, out of the 248 CEGs, were revealed to have occurred in the lineage leading to vertebrates, and were dated at different evolutionary periods by referring to Ensembl Compara [25] (Additional file 4: Figure S3A). These paralogs could potentially cause substantial overestimation of completeness in vertebrate transcriptome or genome sequences. One of the CEGs prone to such effect, the glucose-6-phosphate dehydrogenase (G6PD) gene, is shown as an example in Additional file 4: Figure S3B. Although the latest chicken genome assembly harbors no G6PD ortholog, CEGMA misidentified its ancient duplicate, hexose-6-phosphate dehydrogenase (H6PD), as a G6PD ortholog in the chicken genome. These caveats demonstrate the need of a carefully validated gene set customized for vertebrates to improve the accuracy of completeness assessment.

Instead of using the 248 CEGs, we adopted 233 core vertebrate genes (CVGs; Fig. 3a; Additional file 5). Each of the CVGs is composed of one-to-one orthologs based on eggNOG of 29 species including jawed vertebrates as well as cartilaginous fish and cyclostomes [26] (Fig. 3a; see Methods). All of the 233 CVGs were required to have tunicate orthologs, in addition to having the one-to-one orthologies validated by another ortholog database, Ensembl Compara. These conditions ensured that the one-to-one orthology was retained throughout vertebrate evolution (Fig. 3a), despite large-scale gene (or genome) duplication events in early vertebrates. To consolidate this gene set for CEGMA, eight vertebrate species were selected from the CVG set (Fig. 3b; also Additional file 5): human (Homo sapiens), platypus (Ornithorhynchus anatinus), chicken (Gallus gallus), Xenopus tropicalis, zebrafish (Danio rerio), stickleback (Gasterosteus aculeatus), elephant shark (Callorhinchus milii), and sea lamprey (Petromyzon marinus). Out of the 233 CVGs, only 17 are also included in the 248 CEGs (Additional file 5), and thus there is a substantial difference between the compositions of CVG and CEG. On the other hand, the proportion of the genes identified as ‘complete’ by CEGMA (completeness score) was positively correlated between the CEG and CVG sets (Fig. 3c; R = 0.95; p = 1.0 × 10⁻⁶). Generally, CEGMA referring to the CVG yielded lower completeness scores than with the CEG, and the completeness scores with the CVG showed a higher variation than with the CEG (Fig. 3c). These results indicate that the CVG set has enabled completeness assessment with fewer false positives and higher resolution.

Assessment of assemblies reconstructed from multiple samples

We sequenced eight RNA-seq libraries using total RNA extracted from three embryos at different stages (Fig. 1a; Tables 1 and 2). Because their expression profiles and genetic backgrounds were thought to vary across the samples, we demonstrated several approaches to integrating individual assemblies derived from all of the samples (Fig. 4). All-in-one approach assembles reads from multiple samples at a time (Fig. 4a). In silico normalization removes redundant reads prior to the all-in-one assembly (Fig. 4b), which should lead to reduced computational costs [13, 27]. Another approach is assembly and clustering: assemblies from individual samples are clustered based on sequence similarity (Fig. 4c). We examined five approaches and evaluated their products (See Methods; Assembly 9-13 in Table 2). According to the completeness assessment by CEGMA, the completeness scores of the integrated assemblies varied substantially even though these assemblies were produced from the same set of sequence reads (Fig. 3c and Additional file 6). An integrated assembly is expected to have a completeness score equal to or larger than those of any individual assemblies. This expectation was satisfied only with Assembly 13, which was built with the assembly and clustering approach (Fig. 3c). Assembly 10, based on the all-in-one approach by SOAPdenovo-trans employing multiple k-mer lengths, showed the second largest score. For this assembly however, the completeness score with the CEG did not exceed that of the individual assemblies, whereas it did exceed them with the CVG (Fig. 3c and Additional file 6).

Assembly following in silico normalization can be an alternative approach when sequence data is too large to handle with available computational resources (e.g., more than 200 million reads) [13, 27]. In silico normalization reduced redundancy from the input sequence reads and dramatically decreased the data size down to approximately 10 % (Table 2). However, the completeness scores of Assembly 11 and Assembly 12 were much poorer than even that of Assembly 9, the all-in-one assembly by Trinity.

In order to evaluate the proportion of short read pairs that properly participated in an assembly, we mapped the reads to the assembled contigs and calculated the proportions of the pairs with the ‘properly-paired’ flag in the mapping results (Additional file 7). The result showed that Assembly 10 had the highest mapping rate among the integrated assemblies. In addition, Assembly 10 also exhibited the largest N50 length among all the integrated assemblies, supporting its high overall connectivity (Table 2). Assembly 13 showed the second largest mapping rate and N50 length (Table 2 and Additional file 7). While the qualities of Assembly 10 and Assembly 13 were comparable to each other, we finally chose Assemble 13 as the transcript sequence set representing this study because of its higher completeness score than any other individual assembly. Assembly 13 was composed of 996,336 contigs and 457,323 subcomponents, of which 444,832 were homologous to annotated protein-coding genes of vertebrate genomes that were selected for this purpose (E-value < 1E-10 based on BLASTX; see Methods).

Adapting CVG to the newly introduced completeness assessment tool, BUSCO

The CVG is applicable as a core gene set to not only CEGMA, but widely to other complete assessment methods. We applied the CVG to BUSCO, a recently introduced tool [28] (See Methods). Using the 13 gecko transcriptome assemblies built in this study, we demonstrated that the completeness scores with BUSCO referring to the 233 CVGs were comparable to those with CEGMA, and Assembly 13 again showed the highest completeness score with BUSCO (Additional files 6 and 8).

Originally, BUSCO offers 3,023 ortholog groups of 41 vertebrates that retain one-to-one orthology in almost all of these species based on OrthoDB [29]. However, the BUSCO’s vertebrate gene set is composed of bony vertebrates only, which potentially could fail to identify orthologs of cartilaginous fish and cyclostomes. For example, using the Japanese lamprey (Lethenteron japonicum) genome assembly, BUSCO referring to its original vertebrate gene set showed a completeness score of only 21 %. On the other hand, the completeness scores of BUSCO referring to the metazoan gene set and the CVG (70 % and 73 %, respectively) were comparable to those based on CEGMA (83 % and 77 % referring to the CVG and CEG, respectively). Remarkably, execution of BUSCO referring to the CVG dramatically reduces computational time. Using the Assembly 13, a BUSCO run referring to the CVG is completed in approximately 15 min, while a BUSCO run referring to the vertebrate gene set and a CEGMA run referring to the CVG takes much longer (116 min and 438 min, respectively). We propose the CVG as a handy gene set to be used with BUSCO for completeness assessment of vertebrate genome and transcriptome assemblies.

Coordination between library insert size and read length

In this study, we modified the library preparation protocol for RNA-seq on the Illumina platform to keep up with further developing sequencers. Although the protocol modifications increased the connectivity of the transcriptome assemblies, there still remains room for further improvement. The fragment size distribution of the library prepared with modifications shifted towards the fraction of small-sized fragments (Fig. 2). This shift was notable for the size distribution of Library B: the average size of the actually sequenced fragments was more than 300 bp shorter than that of the fragments in the prepared libraries. This shift of the size distribution was observed regardless of the choice of sequencer models (Fig. 2; Additional file 2). One possible explanation for this is insufficient size selection with AMPure beads. Although its volume was reduced to retain only fragments longer than 300 bp, those shorter than that may not have been sufficiently removed (Fig. 2). The other possible explanation is that ‘clusters’ in Illumina chemistry might not be generated evenly on a flow cell because of variable fragment lengths, although the total size of the output does not differ among sequencer runs for libraries with different fragment length distributions. The Rapid SBS kit v2, which has been released recently, improves robustness against high cluster density through the use of an upgraded HiSeq Control Software (HCS).

Improved completeness assessment based on the 233 CVGs

CEGMA has been a standard for assessing genome assembly and is now widely used for assessment of transcriptome assemblies [30]. However, the presence of the paralogs in the 248 CEGs potentially results in overestimation of the completeness assessment. To improve the accuracy and resolution of the completeness assessment of the gecko transcriptome assemblies, we derived a new gene set composed of one-to-one orthologs of vertebrates (CVG). More conservative completeness scores were computed with CEGMA referring to the 233 CVGs than to the 248 CEGs, suggesting more accurate assessment based on the former (Fig. 3c; Additional file 6). Indeed, none of the CEGMA executions on the 13 assemblies produced false-positive orthologs to the CVGs, as shown in the example of G6PD tree. In contrast, each assessment based on the CEG had approximately 2–5 % of such false positives. False positives in completeness assessment are considered to be more problematic in evaluating transcriptome assemblies than genome assemblies. This is because an absence of a particular transcript in transcriptome assemblies can be caused by no or extremely low gene expression in addition to insufficient sequence read numbers. Establishing a set of one-to-one orthologs for a particular taxonomy group, as demonstrated in the present study, would be applicable to other taxonomic groups.

Members of a one-to-one ortholog group of a specific taxon consequently share high sequence similarity, leading to fewer false positives caused by ancient paralogs being misidentified as orthologs. CEGMA detects genes homologous to a given reference gene sequence and recognizes those satisfying the HMMER score cutoff for the given gene as orthologs. The 233 CVGs have higher HMMER score cutoffs (median, 151.8) for ortholog detection than the 248 CEGs (median, −68.76) since their members are widespread across eukaryotes (Fig. 3b). In addition, orthologs of longer sequences result in a complete assessment with higher resolution because contigs with more than 70 % coverage in length to the HMMER profiles are recognized as ‘complete’ by CEGMA [17]. The larger the lengths of the HMMER profiles are, the less likely the contigs satisfying the ‘completeness’ are reconstructed. Indeed, the HMMER profiles of the CVGs had larger length than those of the 248 CEGs (medians, 557 aa and 379 aa, respectively), leading to the higher resolution of the assessment based on CVG. This high resolution can also be achieved by high sequence similarity among members in a CVG because truncated sequences result in lower HMM scores than the given thresholds, leading to exclusion from ortholog candidates. It is noted that the completeness score of Assembly 13 based on the CVG almost reached 100 %, demonstrating that low completeness scores of the CVG (Fig. 3c) is not caused by the absence of the expression of orthologs in transcriptomes. This suggests that the set of the 233 CVGs better assesses the contents of embryonic transcriptomes of vertebrates, at least.

The new tool, BUSCO, will provide versatile solutions for completeness assessment at various levels. In our analysis, completeness assessment using BUSCO referring to the CVG, which also showed the highest completeness score for Assembly 13, performed similarly to the assessment using CEGMA (Additional files 6 and 8). BUSCO originally offers a vertebrate ortholog set (3,032 groups). This gene set generally performs with similar resolution in completeness assessment to CVG (Additional file 6). In addition, CVG and the BUSCO’s vertebrate ortholog set showed similar performances to identify orthologs in comprehensive transcript sequence sets of diverse bony vertebrates (Additional file 9). Running BUSCO using CVG as a reference has two advantages. The first is that the CVG set includes a cartilaginous fish and a cyclostome, providing a wide taxonomic range for completeness assessment of high sensitivity. Secondly, the fewer components of the CVG allow a much more rapid computation by BUSCO. For public uses, we included the BUSCO-compatible CVG data set in the CVG suite.

Transcriptome assembly of Madagascar ground gecko

Among the five integrated assemblies, Assembly 13 and 10 showed markedly high completeness, connectivity, and accuracy (Fig. 3c; Table 2; Additional files 6 and 7). On the other hand, it was demonstrated that the other assemblies failed to reconstruct considerable numbers of transcripts, which were retained in the individual assemblies (Fig. 3c). Assembly 13 was a merger of the individual assemblies, while Assembly 10 was one of the all-in-one assemblies with different k-mer lengths. Employing multiple k-mer lengths is advantageous in transcriptome assembly because different k-mer lengths show different effective ranges of sequence depths for transcriptome assembly and cover different parts of a transcriptome [5, 9]. The superiority of clustering of individual assemblies to all-in-one assemblies can be explained by variation of expression profiles among samples. This variation of sequence depths for a transcript should increase the possibility of reconstructing a full-length transcript sequence, even based on a single k-mer length.

In general, when typical de novo transcriptome assembly programs are used, the number of resultant contigs constantly increases along with the number of reads fed into the program [30] (Table 2). One possible reason for this is the abundance of so-called ‘leaky’ transcripts from intergenic regions [31] and unprocessed mRNA [32] in addition to contaminated genomic DNA [13]. Such molecules can be contained in the starting materials of library preparation, which will eventually result in low counts of reads after sequencing. To combat this, we excluded contigs with mapping counts of fewer than five in the assembly procedure (see Methods), which reduced the number of contigs of Assembly 13 down to 64 % (Additional file 1). In addition, N50 length of the filtered assembly increased by 336 bp compared with that of Assembly 13. Importantly, the completeness scores of this filtered assembly remained unchanged, suggesting that contigs with low mapping counts rarely include sequences derived from protein-coding genes conserved across vertebrate species. Indeed, in Assembly 13, only 5.6 % of the contigs homologous to the vertebrate genes were included in the discarded contigs. The representative assembly, Assembly 13, is available in the online gateway Reptiliomix linked from our laboratory website (http://www2.clst.riken.jp/phylo/reptiliomix.html), allowing data downloading and similarity searches.

Library preparation and sequencing

All animal experiments and housing were conducted in accordance with guidelines approved by the RIKEN Animal Experiments Committee (Approval IDs AH25-05-1 and AH14-05-47). Embryos of the Madagascar ground gecko, Paroedura picta, were provided by Animal Resource Development Unit, RIKEN CLST. Eggs of 9 and 30 dpo were collected from natural ovipositions, respectively, and an egg of four days before the estimated date of oviposition [−4 days post oviposition (dpo)] was extracted from a egg-bearing female. Total RNA was extracted using TRIzol reagent (Life Technologies), and RNA-seq libraries for non-stranded paired ends were prepared by Illumina TruSeq Total RNA Sample Prep Kit according to its standard protocol unless specifically described below. After confirming the reproducibility in a different species (data not shown), we shortened RNA fragmentation time from eight to two (Library B, C, and E) or four (Library F—H) minutes. In DNA purification, we applied x0.7 (Library B, C, and E) or x1.0 (Library B, C, and E) volume of Agencourt AMPure XP, instead of x1.6. These conditions were summarized in Table 1. The libraries were sequenced with HiSeq 1500 (Illumina inc.) operated by HiSeq Control Software v2.0.12.0 using Rapid SBS kit v1 as well as MiSeq using MiSeq Reagent kit v2, in the read lengths designated in Table 1.

Quality control of sequenced reads

Raw nucleotide bases were called with RTA 1.17.21.3 and converted to the Fastq-format files with bcl2fastq 1.8.3 (Illumina Inc.). The short reads were deposited in the DDBJ Short Read Archive (DRA) database under the accession numbers PRJDB4004. Initially, qualities of sequence reads were checked with FastQC v0.10.1 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/), and no marked abnormalities were observed for all of the samples. Adapter sequences and low quality bases (<Q30) were trimmed from the 3′-ends by trim_galore (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/), in which cutadapt is implemented [33], discarding the reads of shorter than 50 bp after adapter and quality trimming. Low quality reads in which proportion of the bases ≥ Q30 was less than 80 % were discarded by the program fastq_quality_filter in FASTX Toolkit 0.0.13 (http://hannonlab.cshl.edu/fastx_toolkit/index.html). Fragments retaining paired reads were used for de novo assembly.

De novo assembly

The assemblies for individual libraries (Assembly 1–8 in Table 2) were built by Trinity r20131110 [6] which employs a single k-mer length (k = 25). The two all-in-one assemblies were built using the reads of all the libraries by Trinity (Assembly 9) and SOAPdenovo-trans v1.03 [8] based on multiple k-mer lengths (k = 21, 25, 31, 41, 51, 61, 71, 81, and 91) (Assembly 10). For Trinity, we employed default parameters except the ‘group_pairs_distance’ option setting at 1000 for assembling the fragments with long inserts (Assembly 2, 3, 5, 6, 7, 8, 11, and 12). For SOAPdenovo-trans, we set parameters as follows: max_rd_len = 250, rd_len_cutof = 250, avg_ins = 300, reverse_seq = 0, asm_flags = 3, and map_len = 32. The contigs shorter than 200 bp were discarded from Assembly 10, as the default setting of Trinity. In order to remove redundancy of multiple sequences derived presumably from identical transcripts, contigs assembled by SOAPdenovo-trans with multiple k-mer lengths were merged by cd-hit-est v4.6.1 [34] with the similarity threshold of 99 % and word size at eight nucleotides, followed by clustering using gicl v0.0.1 (http://sourceforge.net/projects/gicl/) [35] with the similarity threshold of 95 % and the overlap length threshold of 50 bp. In silico normalizations were performed using the reads of all the libraries by two different programs in order to produce sequences showing unimodal k-mer coverage distributions with averages around 25: khmer 0.2 with options of “-K 20 –C 20” [27] and normalize_by_kmer_coverage.pl implemented in Trinity with options of “--KMER_SIZE 25 --max_cov 50” [13]. The resultant reads from the in silico normalizations were assembled with Trinity into Assembly 11 and Assembly 12. The assemblies for the individual libraries were merged into Assembly 13, as performed for Assembly 10. For Assembly 10 and Assembly 13, which underwent a post-assembly merge, we reassigned gene-transcript relationships by performing single linkage clustering of the contigs derived from the same locus. In order to remove the contigs with minimal read depths that possibly resulted from so-called ‘leaky’ transcription [31], all assemblies were subjected to further modification as follows. The reads were mapped to the contigs in each assembly with Bowtie2 version 2.2.2 end-to-end mode [36]. Using the mapping results, we excluded from the assemblies the contigs on which fewer than five reads were mapped, based on read counts with eXpress v1.5.1 [37]. The mapping and read count were processed using the wrapper align_and_estimate_abundance.pl in Trinity with default parameters except the ‘max_ins_size’ option set at 1500 for Assembly 2-3, and 5-13: the option corresponds to the ‘maxin’ option of Bowtie2 and is set at 800 as default in the wrapper.

In order to confirm that the discarded contigs contain little substantial information for transcriptome analyses, those contig sequences were subjected to searches for possible protein-coding regions with homologs in annotated vertebrate protein databases. For this purpose, similarity searches were carried out using BLASTX [38] in protein sequences of 13 vertebrates with an E-value cutoff of 1E-10: human (Homo sapiens), dog (Canis familiaris), opossum (Monodelphis domestica), chicken (Gallus gallus), zebra finch (Taeniopygia guttata), Chinese alligator (Alligator sinensis), Chinese softshell turtle (Pelodiscus sinensis), green sea turtle (Chelonia mydas), green anole (Anolis carolinensis), Burmese python (Python molurus bivittatus), Xenopus tropicalis, stickleback (Gasterosteus aculeatus), and elephant shark (Callorhinchus milii). Protein sequences of Chinese alligator and green sea turtle were obtained from GigaDB, and those of Burmese python and elephant shark were obtained from NCBI Genbank and Elephant Shark Genome Project web page (http://esharkgenome.imcb.a-star.edu.sg/), respectively. Protein sequences of the other species were downloaded from Ensembl release 75.

Selection of components of CVG

The core vertebrate genes (CVG) were defined as one-to-one orthologs selected based on eggNOG [26] as follows. Initially, we extracted 463 chordate ortholog groups of eggNOG v4.0 (ChorNOGs) [26] that were composed of one-to-one genes of the 26 ‘core’ vertebrates defined by eggNOG, zebrafish (Danio rerio), and sea lamprey (Petromyzon marinus). From them, 292 groups possessing at least one ortholog of either Ciona intestinalis or C. savignyi were selected. Orthologs of elephant shark, whose genome assembly was released later, were added to the gene sets based on the BLASTP [38] reciprocal best-hit approach, and one-to-one elephant shark orthologs were identified in 270 groups out of the 292. The one-to-one orthology of 233 gene sets were systematically validated by gene trees produced by Ensembl release 70 [25], and examined with manual curation when necessary (also see Fig. 3a). Finally we extracted the one-to-one orthologs of eight species, human, platypus (Ornithorhynchus anatinus), chicken, Xenopus tropicalis, zebrafish (Danio rerio), stickleback, elephant shark and sea lamprey from each of the gene sets and defined the selected gene set as 233 CVGs for CEGMA (Additional file 4).

Assessment of assembly

N50 lengths were computed by TrinityStats.pl embedded in Trinity [13]. Short reads were mapped to contigs by Bowtie2 [36], and mapping rates were obtained from its summary output. For completeness assessment, we applied CEGMA version 2.4 [17] based on the 248 CEGs and 233 CVGs separately. Using the CVG gene set for CEGMA, containing the eight vertebrates, HMMER profiles [39] were generated by HMMER 3.0 based on the multiple amino acid alignments processed by MAFFT v7.158b [40] following format conversion into HMMER 2.X. The HMMER bit score cutoffs of the CVGs were computed according to the criterion proposed by Parra et al. [16, 17]: the cutoff values for the standard CEGMA and the completeness analysis, corresponding to the ‘profiles_cutoff.tbl’ and ‘completeness_cutoff.tbl’ files in the original CEGMA package, respectively. We computed the cutoff values of the completeness analysis as maximum hmmsearch bit scores between the HMMER profiles and proteins from any transcripts of the eight species instead of the proteins from the representative transcripts of the genes. In order to conduct the complete assessment using custom gene sets, we modified the scripts of ‘cegma’ and ‘completeness’ implemented in the CEGMA.

We identified false positives in the CEGMA results for the gecko transcriptome assemblies based on a BLASTP search [38]. Using a gecko protein that was predicted as an ortholog of a CVG/CEG by CEGMA, we searched for its best-hit homolog in the human proteins. If this human best-hit was a paralog to the human protein of the CVG/CEG and if these paralogs were duplicated before the split of mammals and sauropsids, the protein predicted by CEGMA was recognized as a misidentified ortholog to the gene group. In this analysis, the human protein sequences and inferred timings of gene duplications were obtained from Ensembl release 70.

In addition to CEGMA, we conducted completeness assessment using BUSCO v1.1 [28] referring to the CVG, and the CVG dataset for BUSCO was prepared as follows. Using the whole sequence set (29 species) of the CVG, HMMER profiles [39] were generated by HMMER 3.1b2 based on the multiple amino acid alignments processed by MAFFT v7.158b [40]. Protein profiles of the CVG for Augustus were generated with msa2prfl.pl in Augustus 3.1 [41] based on the multiple alignments. The consensus sequence of each CVG was inferred by hmmemit in the HMMER. The cutoff values of sequence lengths were computed according to the criterion described previously [28]. As for cutoff values of HMMER bit scores, we used the values for CEGMA described above, instead of those according to the original criterion by BUSCO.

Mapping rates of the assemblies were computed by SAMtools version 0.1.19 [42] using the mapping files that were made for counting mapped reads to the contigs in the previous subsection. Insert lengths of the fragments were estimated with CollectInsertSizeMetrics in Picard Tools version 1.90 (http://broadinstitute.github.io/picard/). For this purpose, in mapping, we used first 50 nucleotides of each paired read so that paired-reads were not mapped overlapping each other. These reads were mapped to the assemblies using Bowtie2 with the same parameters described in the subsection “De novo assembly” in Methods.

Molecular phylogenetic analysis

Peptide sequences of G6PD and H6PD were collected from the gene set of KOG0563 in CEGMA as well as NCBI Genbank and Ensembl release 70 with the assistance of aLeaves [43]. The homologous peptides were aligned with six different approaches: forward and reverse directions by MAFFT v7.158b [40], Clustal Omega 1.2.0 [44], and T-Coffee 10.00.r1613A [45]. The consensus multiple alignment from the six procedures was made by M-Coffee [46] implemented in the T-Coffee package. Unambiguous alignment sites were selected based on trimAl version 1.4 with the automated1 option following removal of gapped sites. Molecular phylogenetic trees were reconstructed based on RAxML version 7.5.7 [47] assuming the PROTCATWAG model with 1,000 bootstrap replicates ("−f a" option) and PhyloBayes 3.3f assuming the CAT-GTR model [48].