Revising transcriptome assemblies with phylogenetic information

Abstract

A common transcriptome assembly error is to mistake different transcripts of the same gene as transcripts from multiple closely related genes. This error is difficult to identify during assembly, but in a phylogenetic analysis such errors can be diagnosed from gene phylogenies where they appear as clades of tips from the same species with improbably short branch lengths. treeinform is a method that uses phylogenetic information across species to refine transcriptome assemblies within species. It identifies transcripts of the same gene that were incorrectly assigned to multiple genes and reassign them as transcripts of the same gene. The treeinform method is implemented in Agalma, available at https://bitbucket.org/caseywdunn/agalma, and the general approach is relevant in a variety of other contexts.

Fig 1. An example gene phylogeny from the test dataset before running treeinform.

Each tip is an exemplar transcript that was initially assigned to a different gene. In front, corresponding multiple sequence alignment, with sites ordered from highest to lowest identity to the inferred ancestral site for clarity on sequence diversity. Black indicates a difference from the ancestral sequence. The four Hydra transcripts in color were assigned to different genes by Trinity [2] despite two of the transcripts sharing the exact same sequences, and the two other transcripts differing by a small gap. After treeinform, all transcripts from these four genes are reassigned to a single gene.

Fig 2. Histogram of subtree lengths for internal nodes in each gene phylogeny from the test dataset containing tip descendants from the same species.

Subtree lengths greater than 1 were filtered out for clarity.

Fig 3. Cluster size counts for Trinity assembly and Corset clustering algorithm on Trinity contigs.

There are 3 Trinity clusters with size greater than 30, while there are 20 Corset clusters with size greater than 30.

Fig 4. Histogram of subtree lengths for internal nodes in each Siphonophora subset gene tree from Agalma with Corset clusterings containing tip descendants from the same species.

Subtree lengths greater than 1 were filtered out for clarity.

Fig 5. Histogram of the inferred duplication times with an overlaid mixture model.

Component 1 of the mixture model (red) captures the technical issues we address here, where transcripts from the same gene are assigned to different genes, and component 2 (blue) captures the true biological pattern, where transcripts from different genes are correctly assigned to different genes. We first ran phyldog [16] on the test dataset using the multiple sequence alignments and a given species phylogeny [20]. This provided gene phylogenies with internal nodes annotated as duplication or speciation events. We then used the annotations to time-calibrate the gene phylogenies for the mixture model.

Fig 6. Percentage of reassigned transcripts (log scale).

47,688 genes were included in the gene phylogenies, of which 23,396 (49.06%) were in gene families of 2 or more, and thus candidates for reassignment. The default threshold for treeinform is marked by the grey vertical dashed line.

Fig 7. Density from theoretical and the empirical density under 3 different thresholds before treeinform was run.

The distribution before treeinform has a large peak on the left that is removed by treeinform with all examined thresholds. Black line represents theoretical density.

Fig 8. Histogram of subtree lengths for internal nodes in each Drosophila and Echinoidea subset gene tree from Agalma containing tip descendants from the same species.

Top is Drosophila and bottom is Echinoidea. Subtree lengths greater than 1 were filtered out for clarity.

Fig 9. Precision vs. recall for pairs of transcripts with regards to known CDS as treeinform threshold increases.

Top: Plot for Drosophila with CDS. Precision does not increase with any threshold, only recall. Bottom: Plot for Strongylocentrotus purpuratus. The biggest improvement is made at a threshold value 5e-06, with precision and recall both increasing up to a threshold of 0.05.