Comprehensive evaluation of RNA-seq analysis pipelines in diploid and polyploid species

Machine learning on alignment features for parent-of-origin classification of simulated hybrid RNA-seq

BACKGROUND

Parent-of-origin allele-specific gene expression (ASE) can be detected in interspecies hybrids by virtue of RNA sequence variants between the parental haplotypes. ASE is detectable by differential expression analysis (DEA) applied to the counts of RNA-seq read pairs aligned to parental references, but aligners do not always choose the correct parental reference.

RESULTS

We used public data for species that are known to hybridize. We measured our ability to assign RNA-seq read pairs to their proper transcriptome or genome references. We tested software packages that assign each read pair to a reference position and found that they often favored the incorrect species reference. To address this problem, we introduce a post process that extracts alignment features and trains a random forest classifier to choose the better alignment. On each simulated hybrid dataset tested, our machine-learning post-processor achieved higher accuracy than the aligner by itself at choosing the correct parent-of-origin per RNA-seq read pair.

CONCLUSIONS

For the parent-of-origin classification of RNA-seq, machine learning can improve the accuracy of alignment-based methods. This approach could be useful for enhancing ASE detection in interspecies hybrids, though RNA-seq from real hybrids may present challenges not captured by our simulations. We believe this is the first application of machine learning to this problem domain.

Transcriptomics of mussel transmissible cancer MtrBTN2 suggests accumulation of multiple cancer traits and oncogenic pathways shared among bilaterians

Transmissible cancer cell lines are rare biological entities giving rise to diseases at the crossroads of cancer and parasitic diseases. These malignant cells have acquired the amazing capacity to spread from host to host. They have been described only in dogs, Tasmanian devils and marine bivalves. The Mytilus trossulus bivalve transmissible neoplasia 2 (MtrBTN2) lineage has even acquired the capacity to spread inter-specifically between marine mussels of the Mytilus edulis complex worldwide. To identify the oncogenic processes underpinning the biology of these atypical cancers we performed transcriptomics of MtrBTN2 cells. Differential expression, enrichment, protein-protein interaction network, and targeted analyses were used. Overall, our results suggest the accumulation of multiple cancerous traits that may be linked to the long-term evolution of MtrBTN2. We also highlight that vertebrate and lophotrochozoan cancers could share a large panel of common drivers, which supports the hypothesis of an ancient origin of oncogenic processes in bilaterians.

How genomics can help biodiversity conservation

The availability of public genomic resources can greatly assist biodiversity assessment, conservation, and restoration efforts by providing evidence for scientifically informed management decisions. Here we survey the main approaches and applications in biodiversity and conservation genomics, considering practical factors, such as cost, time, prerequisite skills, and current shortcomings of applications. Most approaches perform best in combination with reference genomes from the target species or closely related species. We review case studies to illustrate how reference genomes can facilitate biodiversity research and conservation across the tree of life. We conclude that the time is ripe to view reference genomes as fundamental resources and to integrate their use as a best practice in conservation genomics.

Optimization and application of non-native Phragmites australis transcriptome assemblies

Phragmites australis (common reed) has a cosmopolitan distribution and has been suggested as a model organism for the study of invasive plant species. In North America, the non-native subspecies (ssp. australis) is widely distributed across the contiguous 48 states in the United States and large parts of Canada. Even though millions of dollars are spent annually on Phragmites management, insufficient knowledge of P. australis impeded the efficiency of management. To solve this problem, transcriptomic information generated from multiple types of tissue could be a valuable resource for future studies. Here, we constructed forty-nine P. australis transcriptomes assemblies via different assembly tools and multiple parameter settings. The optimal transcriptome assembly for functional annotation and downstream analyses was selected among these transcriptome assemblies by comprehensive assessments. For a total of 422,589 transcripts assembled in this transcriptome assembly, 319,046 transcripts (75.5%) have at least one functional annotation. Within the transcriptome assembly, we further identified 1,495 transcripts showing tissue-specific expression pattern, 10,828 putative transcription factors, and 72,165 candidates for simple sequence repeats markers. The identification and analyses of predicted transcripts related to herbicide- and salinity-resistant genes were shown as two applications of the transcriptomic information to facilitate further research on P. australis. Transcriptome assembly and selection would be important for the transcriptome annotation. With this optimal transcriptome assembly and all relative information from downstream analyses, we have helped to establish foundations for future studies on the mechanisms underlying the invasiveness of non-native P. australis subspecies.

Unleashing the power within short-read RNA-seq for plant research: Beyond differential expression analysis and toward regulomics

RNA-seq has become a state-of-the-art technique for transcriptomic studies. Advances in both RNA-seq techniques and the corresponding analysis tools and pipelines have unprecedently shaped our understanding in almost every aspects of plant sciences. Notably, the integration of huge amount of RNA-seq with other omic data sets in the model plants and major crop species have facilitated plant regulomics, while the RNA-seq analysis has still been primarily used for differential expression analysis in many less-studied plant species. To unleash the analytical power of RNA-seq in plant species, especially less-studied species and biomass crops, we summarize recent achievements of RNA-seq analysis in the major plant species and representative tools in the four types of application: (1) transcriptome assembly, (2) construction of expression atlas, (3) network analysis, and (4) structural alteration. We emphasize the importance of expression atlas, coexpression networks and predictions of gene regulatory relationships in moving plant transcriptomes toward regulomics, an omic view of genome-wide transcription regulation. We highlight what can be achieved in plant research with RNA-seq by introducing a list of representative RNA-seq analysis tools and resources that are developed for certain minor species or suitable for the analysis without species limitation. In summary, we provide an updated digest on RNA-seq tools, resources and the diverse applications for plant research, and our perspective on the power and challenges of short-read RNA-seq analysis from a regulomic point view. A full utilization of these fruitful RNA-seq resources will promote plant omic research to a higher level, especially in those less studied species.

Genome-Wide Identification and Characterization of RNA/DNA Differences Associated with Fusarium graminearum Infection in Wheat

RNA/DNA difference (RDD) is a post-transcriptional modification playing a crucial role in regulating diverse biological processes in eukaryotes. Although it has been extensively studied in plant chloroplast and mitochondria genomes, RDDs in plant nuclear genomes are not well studied at present. Here, we investigated the RDDs associated with fusarium head blight (FHB) through a novel method by comparing the RNA-seq data between Fusarium-infected and control samples of four wheat genotypes. A total of 187 high-confidence unique RDDs in 36 genes were identified, representing the first landscape of the FHB-responsive RDD in wheat. The majority (26) of these 36 RDD genes were correlated either positively or negatively with FHB levels. Effects of these RDDs on RNA and protein sequences have been identified, their editing frequency and the expression level of the corresponding genes provided, and the prediction of the effect on the minimum folding free energy of mRNA, miRNA binding, and colocation of RDDs with conserved domains presented. RDDs were predicted to induce modifications in the mRNA and protein structures of the corresponding genes. In two genes, TraesCS1B02G294300 and TraesCS3A02G263900, editing was predicted to enhance their affinity with tae-miR9661-5p and tae-miR9664-3p, respectively. To our knowledge, this study is the first report of the association between RDD and FHB in wheat; this will contribute to a better understanding of the molecular basis underlying FHB resistance, and potentially lead to novel strategies to improve wheat FHB resistance through epigenetic methods.

A method for simultaneous detection of small and long RNA biotypes by ribodepleted RNA-Seq

RNA sequencing offers unprecedented access to the transcriptome. Key to this is the identification and quantification of many different species of RNA from the same sample at the same time. In this study we describe a novel protocol for simultaneous detection of coding and non-coding transcripts using modifications to the Ion Total RNA-Seq kit v2 protocol, with integration of QIASeq FastSelect rRNA removal kit. We report highly consistent sequencing libraries can be produced from both frozen high integrity mouse hippocampal tissue and the more challenging post-mortem human tissue. Removal of rRNA using FastSelect was extremely efficient, resulting in less than 1.5% rRNA content in the final library. We identified > 30,000 unique transcripts from all samples, including protein-coding genes and many species of non-coding RNA, in biologically-relevant proportions. Furthermore, the normalized sequencing read count for select genes significantly negatively correlated with Ct values from qRT-PCR analysis from the same samples. These results indicate that this protocol accurately and consistently identifies and quantifies a wide variety of transcripts simultaneously. The highly efficient rRNA depletion, coupled with minimized sample handling and without complicated and high-loss size selection protocols, makes this protocol useful to researchers wishing to investigate whole transcriptomes.

Trimming and Validation of Illumina Short Reads Using Trimmomatic, Trinity Assembly, and Assessment of RNA-Seq Data

Next-generation sequencing (NGS) technologies can generate billions of reads in a single sequencing run. However, with such high-throughput comes quality issues which have to be addressed before undertaking downstream analysis. Quality control on short reads is usually performed at default settings due to a lack of in-depth understanding of a particular software's parameters and their effect if changed on the output. Here we demonstrate how to optimize read trimming using Trimmomatic. We highlight the benefits of trimming by comparing the quality of transcripts assembled using trimmed and untrimmed reads.

Improving the quality of barley transcriptome de novo assembling by using a hybrid approach for lines with varying spike and stem coloration

De novo transcriptome assembly is an important stage of RNA-seq data computational analysis. It allows the researchers to obtain the sequences of transcripts presented in the biological sample of interest. The availability of accurate and complete transcriptome sequence of the organism of interest is, in turn, an indispensable condition for further analysis of RNA-seq data. Through years of transcriptomic research, the bioinformatics community has developed a number of assembler programs for transcriptome reconstruction from short reads of RNA-seq libraries. Different assemblers makes it possible to conduct a de novo transcriptome reconstruction and a genome-guided reconstruction. The majority of the assemblers working with RNA-seq data are based on the De Bruijn graph method of sequence reconstruction. However, specif ics of their procedures can vary drastically, as do their results. A number of authors recommend a hybrid approach to transcriptome reconstruction based on combining the results of several assemblers in order to achieve a better transcriptome assembly. The advantage of this approach has been demonstrated in a number of studies, with RNA-seq experiments conducted on the Illumina platform. In this paper, we propose a hybrid approach for creating a transcriptome assembly of the barley Hordeum vulgare isogenic line Bowman and two nearly isogenic lines contrasting in spike pigmentation, based on the results of sequencing on the IonTorrent platform. This approach implements several de novo assemblers: Trinity, Trans-ABySS and rnaSPAdes. Several assembly metrics were examined: the percentage of reference transcripts observed in the assemblies, the percentage of RNA-seq reads involved, and BUSCO scores. It was shown that, based on the summation of these metrics, transcriptome meta-assembly surpasses individual transcriptome assemblies it consists of.

Comparing de novo transcriptome assembly tools in di- and autotetraploid non-model plant species

Background

Polyploidy is very common in plants and can be seen as one of the key drivers in the domestication of crops and the establishment of important agronomic traits. It can be the main source of genomic repatterning and introduces gene duplications, affecting gene expression and alternative splicing. Since fully sequenced genomes are not yet available for many plant species including crops, de novo transcriptome assembly is the basis to understand molecular and functional mechanisms. However, in complex polyploid plants, de novo transcriptome assembly is challenging, leading to increased rates of fused or redundant transcripts. Since assemblers were developed mainly for diploid organisms, they may not well suited for polyploids. Also, comparative evaluations of these tools on higher polyploid plants are extremely rare. Thus, our aim was to fill this gap and to provide a basic guideline for choosing the optimal de novo assembly strategy focusing on autotetraploids, as the scientific interest in this type of polyploidy is steadily increasing.

Results

We present a comparison of two common (SOAPdenovo-Trans, Trinity) and one recently published transcriptome assembler (TransLiG) on diploid and autotetraploid species of the genera Acer and Vaccinium using Arabidopsis thaliana as a reference. The number of assembled transcripts was up to 11 and 14 times higher with an increased number of short transcripts for Acer and Vaccinium, respectively, compared to A. thaliana. In diploid samples, Trinity and TransLiG performed similarly good while in autotetraploids, TransLiG assembled most complete transcriptomes with an average of 1916 assembled BUSCOs vs. 1705 BUSCOs for Trinity. Of all three assemblers, SOAPdenovo-Trans performed worst (1133 complete BUSCOs).

Conclusion

All three assembly tools produced complete assemblies when dealing with the model organism A. thaliana, independently of its ploidy level, but their performances differed extremely when it comes to non-model autotetraploids, where specifically TransLiG and Trinity produced a high number of redundant transcripts. The recently published assembler TransLiG has not been tested yet on any plant organism but showed highest completeness and full-length transcriptomes, especially in autotetraploids. Including such species during the development and testing of new assembly tools is highly appreciated and recommended as many important crops are polyploid.

Supplementary Information

The online version contains supplementary material available at 10.1186/s12859-021-04078-8.

Validation of reference genes for quantitative gene expression in the Lippia alba polyploid complex (Verbenaceae)

Lippia alba (Verbenaceae) is one of the most studied species of the genus Lippia, mainly due to its medicinal properties. The species was described as a polyploid complex with five cytotypes. The comparison of gene expression in species with several ploidal levels needs to be conducted carefully due to possible changes in gene regulation. Quantitative reverse transcription PCR (qRT-PCR) is a widely used method for transcript abundance analyses in plants. Besides being an extremely powerful technique, relative quantification by Real-Time quantitative PCR (RT-qPCR) needs the normalization with a stable reference gene. We evaluated the stability of nine candidate reference genes in Lippia alba with different ploidal levels using NormFinder, geNorm, and RefFinder software. The product of each primer showed a single peak in the melting curve. The R² value ranged from 0.998 to 1000 and primers efficiency ranged from 98.95% to 129%. The CIT gene came up as a stable housekeeping gene, being appropriate for studies in polyploid accessions of Lippia alba. Considering that polyploidy is widely documented in Angiosperms, the results can be used not only for further gene expression studies in L. alba but also as a possible reference gene for other polyploid complexes. Differential stability among different genes highlights the importance of the validation of reference genes used for RT-qPCR approach in polyploid studies.

RNA-Seq in Nonmodel Organisms

RNA-Seq is nowadays an indispensable approach for comparative transcriptome profiling in model and nonmodel organisms. Analyzing RNA-Seq data from nonmodel organisms poses unique challenges, due to unavailability of a high-quality genome reference and to relative sparsity of tools for downstream functional analyses. In this chapter, we provide an overview of the analysis steps in RNA-Seq projects of nonmodel organisms, while elaborating on aspects that are unique to this analysis. These will include (1) strategic decisions that have to be made in advance, regarding sequencing technology and reference to use; (2) how to search for available draft genomes, and, if necessary, how to improve their gene prediction and annotation; (3) how to clean raw reads before de novo assembly; (4) how to separate the reads in RNA-Seq projects of symbiont organisms; (5) how to design and carry out a de novo transcriptome assembly that will be comprehensive and reliable; (6) how to assess transcriptome quality; (7) when and how to reduce redundancy in the transcriptome; (8) techniques and considerations in transcriptome functional annotation; (9) quantitating transcript abundance in the face of high transcriptome redundancy; and, most importantly, (10) how to achieve functional enrichment testing using available tools which either support a large range of species or enable a universal, non-species-specific analysis.Throughout the chapter, we will refer to a variety of useful software tools. For the initial analysis steps involving high-volume data, these will include Linux-based programs. For the later steps, we will describe both Linux and R packages for advanced users, as well as many user-friendly tools for nonprogrammers. Finally, we will present a full workflow for RNA-Seq analysis of nonmodel organisms using the NeatSeq-Flow platform, which can be used locally through a user-friendly interface.

Increased biological relevance of transcriptome analyses in human skeletal muscle using a model-specific pipeline

Background

Human skeletal muscle responds to weight-bearing exercise with significant inter-individual differences. Investigation of transcriptome responses could improve our understanding of this variation. However, this requires bioinformatic pipelines to be established and evaluated in study-specific contexts. Skeletal muscle subjected to mechanical stress, such as through resistance training (RT), accumulates RNA due to increased ribosomal biogenesis. When a fixed amount of total-RNA is used for RNA-seq library preparations, mRNA counts are thus assessed in different amounts of tissue, potentially invalidating subsequent conclusions. The purpose of this study was to establish a bioinformatic pipeline specific for analysis of RNA-seq data from skeletal muscles, to explore the effects of different normalization strategies and to identify genes responding to RT in a volume-dependent manner (moderate vs. low volume). To this end, we analyzed RNA-seq data derived from a twelve-week RT intervention, wherein 25 participants performed both low- and moderate-volume leg RT, allocated to the two legs in a randomized manner. Bilateral muscle biopsies were sampled from m. vastus lateralis before and after the intervention, as well as before and after the fifth training session (Week 2).

Result

Bioinformatic tools were selected based on read quality, observed gene counts, methodological variation between paired observations, and correlations between mRNA abundance and protein expression of myosin heavy chain family proteins. Different normalization strategies were compared to account for global changes in RNA to tissue ratio. After accounting for the amounts of muscle tissue used in library preparation, global mRNA expression increased by 43–53%. At Week 2, this was accompanied by dose-dependent increases for 21 genes in rested-state muscle, most of which were related to the extracellular matrix. In contrast, at Week 12, no readily explainable dose-dependencies were observed. Instead, traditional normalization and non-normalized models resulted in counterintuitive reverse dose-dependency for many genes. Overall, training led to robust transcriptome changes, with the number of differentially expressed genes ranging from 603 to 5110, varying with time point and normalization strategy.

Conclusion

Optimized selection of bioinformatic tools increases the biological relevance of transcriptome analyses from resistance-trained skeletal muscle. Moreover, normalization procedures need to account for global changes in rRNA and mRNA abundance.

Performance of gene expression analyses using de novo assembled transcripts in polyploid species

MOTIVATION

Quality of gene expression analyses using de novo assembled transcripts in species that experienced recent polyploidization remains unexplored.

RESULTS

Differential gene expression (DGE) analyses using putative genes inferred by Trinity, Corset and Grouper performed slightly differently across five plant species that experienced various polyploidy histories. In species that lack recent polyploidy events that occurred in the past several millions of years, DGE analyses using de novo assembled transcriptomes identified 54-82% of the differentially expressed genes recovered by mapping reads to the reference genes. However, in species that experienced more recent polyploidy events, the percentage decreased to 21-65%. Gene co-expression network analyses using de novo assemblies versus mapping to the reference genes recovered the same module that significantly correlated with treatment in one species that lacks recent polyploidization.

AVAILABILITY AND IMPLEMENTATION

Commands and scripts used in this study are available at https://bitbucket.org/lychen83/chen_et_al_2018_benchmark_dge/; Analysis files are available at Dryad doi: 10.5061/dryad.4p6n481.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Visual pigment evolution in Characiformes: The dynamic interplay of teleost whole-genome duplication, surviving opsins and spectral tuning

Vision represents an excellent model for studying adaptation, given the genotype-to-phenotype map that has been characterized in a number of taxa. Fish possess a diverse range of visual sensitivities and adaptations to underwater light, making them an excellent group to study visual system evolution. In particular, some speciose but understudied lineages can provide a unique opportunity to better understand aspects of visual system evolution such as opsin gene duplication and neofunctionalization. In this study, we showcase the visual system evolution of neotropical Characiformes and the spectral tuning mechanisms they exhibit to modulate their visual sensitivities. Such mechanisms include gene duplications and losses, gene conversion, opsin amino acid sequence and expression variation, and A₁ /A₂ -chromophore shifts. The Characiforms we studied utilize three cone opsin classes (SWS2, RH2, LWS) and a rod opsin (RH1). However, the characiform's entire opsin gene repertoire is a product of dynamic evolution by opsin gene loss (SWS1, RH2) and duplication (LWS, RH1). The LWS- and RH1-duplicates originated from a teleost specific whole-genome duplication as well as characiform-specific duplication events. Both LWS-opsins exhibit gene conversion and, through substitutions in key tuning sites, one of the LWS-paralogues has acquired spectral sensitivity to green light. These sequence changes suggest reversion and parallel evolution of key tuning sites. Furthermore, characiforms' colour vision is based on the expression of both LWS-paralogues and SWS2. Finally, we found interspecific and intraspecific variation in A₁ /A₂ -chromophores proportions, correlating with the light environment. These multiple mechanisms may be a result of the diverse visual environments where Characiformes have evolved.

The Utility of Genomic and Transcriptomic Data in the Construction of Proxy Protein Sequence Databases for Unsequenced Tree Nuts

As the apparent incidence of tree nut allergies rises, the development of MS methods that accurately identify tree nuts in food is critical. However, analyses are limited by few available tree nut protein sequences. We assess the utility of translated genomic and transcriptomic data for library construction with Juglans regia, walnut, as a model. Extracted walnuts were subjected to nano-liquid chromatography-mass spectrometry (n-LC-MS/MS), and spectra were searched against databases made from a six-frame translation of the genome (6FT), a transcriptome, and three proteomes. Searches against proteomic databases yielded a variable number of peptides (1156-1275), and only ten additional unique peptides were identified in the 6FT database. Searches against a transcriptomic database yielded results similar to those of the National Center for Biotechnology Information (NCBI) proteome (1200 and 1275 peptides, respectively). Performance of the transcriptomic database was improved via the adjustment of RNA-Seq read processing methods, which increased the number of identified peptides which align to seed allergen proteins by ~20%. Together, these findings establish a path towards the construction of robust proxy protein databases for tree nut species and other non-model organisms.

Homoeologous gene expression and co-expression network analyses and evolutionary inference in allopolyploids

Polyploidy is a widespread phenomenon throughout eukaryotes. Due to the coexistence of duplicated genomes, polyploids offer unique challenges for estimating gene expression levels, which is essential for understanding the massive and various forms of transcriptomic responses accompanying polyploidy. Although previous studies have explored the bioinformatics of polyploid transcriptomic profiling, the causes and consequences of inaccurate quantification of transcripts from duplicated gene copies have not been addressed. Using transcriptomic data from the cotton genus (Gossypium) as an example, we present an analytical workflow to evaluate a variety of bioinformatic method choices at different stages of RNA-seq analysis, from homoeolog expression quantification to downstream analysis used to infer key phenomena of polyploid expression evolution. In general, EAGLE-RC and GSNAP-PolyCat outperform other quantification pipelines tested, and their derived expression dataset best represents the expected homoeolog expression and co-expression divergence. The performance of co-expression network analysis was less affected by homoeolog quantification than by network construction methods, where weighted networks outperformed binary networks. By examining the extent and consequences of homoeolog read ambiguity, we illuminate the potential artifacts that may affect our understanding of duplicate gene expression, including an overestimation of homoeolog co-regulation and the incorrect inference of subgenome asymmetry in network topology. Taken together, our work points to a set of reasonable practices that we hope are broadly applicable to the evolutionary exploration of polyploids.

Comprehensive evaluation of RNA-seq analysis pipelines in diploid and polyploid species

Background

The usual analysis of RNA sequencing (RNA-seq) reads is based on an existing reference genome and annotated gene models. However, when a reference for the sequenced species is not available, alternatives include using a reference genome from a related species or reconstructing transcript sequences with de novo assembly. In addition, researchers are faced with many options for RNA-seq data processing and limited information on how their decisions will impact the final outcome. Using both a diploid and polyploid species with a distant reference genome, we have tested the influence of different tools at various steps of a typical RNA-seq analysis workflow on the recovery of useful processed data available for downstream analysis.

Findings

At the preprocessing step, we found error correction has a strong influence on de novo assembly but not on mapping results. After trimming, a greater percentage of reads could be used in downstream analysis by selecting gentle quality trimming performed with Skewer instead of strict quality trimming with Trimmomatic. This availability of reads correlated with size, quality, and completeness of de novo assemblies and with number of mapped reads. When selecting a reference genome from a related species to map reads, outcome was significantly improved when using mapping software tolerant of greater sequence divergence, such as Stampy or GSNAP.

Conclusions

The selection of bioinformatic software tools for RNA-seq data analysis can maximize quality parameters on de novo assemblies and availability of reads in downstream analysis.