TrAnnoScope: A Modular Snakemake Pipeline for Comprehensive Full-Length Transcriptome Analysis and Functional Annotation

1. Introduction

RNA sequencing (RNA-Seq) has become a powerful tool for detecting novel transcripts, understanding gene expression, cataloging protein-coding genes, and revealing the biological functions of genes [1,2,3]. Additionally, RNA-Seq has unlocked the study of non-model organisms without the need for a reference genome, through de novo transcriptome analysis [4]. However, there are several challenges associated with RNA-Seq analysis, such as sequencing errors and fragmentation resulting from technological limitations, in addition to issues, such as repetitive regions and overlapping genes due to transcriptome complexity [2,4].

Short-read technologies have lower error rates and provide a higher coverage than long-read sequencing technologies. Nevertheless, transcriptomes generated exclusively from short reads often suffer from fragmentation and incomplete transcript reconstructions due to the erroneous computational predictions of isoforms. In contrast, long-read sequencing technologies can capture full-length (FL) transcripts and resolve isoform complexity; however, they retain a higher error rate and lower throughput. Accurate transcriptome assembly is crucial for downstream analyses, including functional genomics, gene discovery, and the elucidation of complex biological processes [4]. Hybrid approaches that combine short- and long-read technologies can overcome the weaknesses of each technology and improve transcriptome coverage and accuracy to obtain known and novel transcripts [2,5].

Transcriptome generation and annotation are challenging because of the complexity of the procedures, the need to select appropriate tools, and the significant computational resources required [3]. Several RNA-Seq pipelines offer an interconnected collection of tools designed to automate the process, such as RNAflow [6] and RASflow [7], which primarily focus on differential expression analysis, while others, such as TransXpress [8], TransPi [9], and Pincho [10], focus on de novo transcriptome assembly and functional annotation. However, these tools depend on short-read sequencing for analysis.

Several toolkits, such as Functional IsoTranscriptomics Analysis (FIT) [11], IsoTools [12], TAGET [13], and nf-core/isoseq [14] utilize the properties of long-read sequencing technologies for transcriptome analysis. However, they are primarily designed to function with reference annotations, and currently, Trans2Express [15] is the only reproducible protocol for non-model organisms. It enables de novo hybrid transcriptome assembly using both the Illumina and Oxford Nanopore Technologies (ONT) platforms, aiming to recover a single transcript per gene for transcriptome characterization and gene expression analysis. However, this approach may lead to the loss of important information relating to alternative splicing and isoform diversity and limit the detection of novel transcripts or isoforms, especially in the less-studied regions of a transcriptome; a more comprehensive approach that captures multiple isoforms is essential for fully elucidating the functional potential of genes [16]. Furthermore, Trans2Express offers limited flexibility, restricting users from selecting and combining tools within the pipeline to meet their specific research objectives.

Here, we present TrAnnoScope, a comprehensive FL transcriptome and annotation pipeline that integrates Illumina short-and long-read data through a number of key steps, including data preprocessing, long-read error correction, contamination removal, quality assessment, and comprehensive functional annotation. The pipeline was designed to improve transcriptome accuracy and completeness by leveraging the strengths of long reads using short reads and thorough preprocessing without the need for reference annotation. This integration can resolve complex isoform structures, identify novel genes, and provide a complete transcriptome representation. TrAnnoScope is highly modular, allowing users to customize workflows and integrate different components to suit their research goals. In addition, it supports parallel execution and cluster computing, enabling the faster processing of larger datasets. It is sufficiently versatile to be used in many research areas, from gene discovery and transcriptome profiling to comparative genomics and the study of complex biological systems. As a result, TrAnnoScope provides a powerful, reproducible approach for conducting efficient bioinformatics analyses for large-scale transcriptome analyses.

2. Materials and Methods

2.1. Components of the TrAnnoScope Pipeline

We implemented our pipeline using Snakemake owing to its simplicity and ability to automate complex workflows while managing dependencies [17]. The TrAnnoScope pipeline consists of several modules, starting with a Python script (Figure 1A) that simplifies the setup by installing the necessary dependencies and databases for TrAnnoScope, allowing users to install only the components essential for their needs. Following this, the pipeline includes quality control of Illumina reads (Figure 1B), preprocessing of Illumina reads (Figure 1E), preprocessing of PacBio reads (Figure 1C), contamination removal, error correction using Illumina reads, isoform clustering and classification (Figure 1F), quality assessment (Figure 1D), and annotation (Figure 1G). A configuration file is provided for users to customize the tool parameters according to their needs. To address the time-intensive steps, such as contamination removal and annotation, TrAnnoScope supports parallel execution by dividing input files for faster processing. Detailed instructions are available on the TrAnnoScope GitHub page, [https://github.com/aysevllpkts/TrAnnoScope accessed on 6 Nov 2024], and users can run the pipeline fully or selectively based on their needs.

2.2. Pipeline Implementation

2.2.1. Quality and Preprocessing for Illumina Reads

The first module of the TrAnnoScope pipeline focuses on performing the quality control of Illumina short reads, which is an essential step for ensuring data reliability. FastQC [18] is used to generate quality reports for individual samples, and MultiQC [19] combines these reports into a comprehensive summary of all Illumina datasets.

To detect potential contamination, FastQScreen [20] maps the reads against a set of reference databases. TrAnnoScope allows users to use the default FastQScreen Genome database or create a custom database to filter out undesired sequences, such as species-specific mitochondrial RNA (mtRNA) or ribosomal RNA (rRNA) sequences. For adapter trimming and quality filtering, fastp [21] is employed due to its speed, low memory usage, and detailed quality reports. Additionally, it detects and removes poly-G tails from Illumina NextSeq/NovaSeq data. Ensuring high-quality reads, such as via contamination removal and error correction, is crucial for downstream analysis, as poor-quality reads can lead to inaccurate results. A final quality control step is performed for preprocessed reads using FastQC and MultiQC to confirm that they are suitable for downstream analysis.

2.2.2. Preprocessing of PacBio Reads

The Isoseq3 package [22] is used to preprocess the PacBio reads, which is essential for converting raw sequencing data into high-quality, FL transcripts. The module begins by converting subreads into circular consensus sequences (CCSs), which are self-corrected sequences, using the ccs tool to ensure a higher accuracy in the downstream analysis. The lima tool is then used to demultiplex and trim primer sequences from the CCS reads to obtain FL reads. Further, isoseq3 refine removes chimeric sequences, which can originate from multiple transcripts, and retains only reads with poly-A tails, producing FL non-chimeric (FLNC) reads that improve read quality. Isoseq3 cluster/cluster2 groups reads based on their sequence similarity to identify unique transcripts and their isoforms, offering insights into alternative splicing and transcript diversity, both of which are crucial for understanding gene expression dynamics. This preprocessing step currently applies only to the PacBio reads. However, if the user has Nanopore reads in the FASTA format, they can still be used for the subsequent steps of the pipeline.

2.2.3. Contamination Removal

Contaminants in RNA-Seq data can significantly impact the quality and accuracy of results, leading to biased gene expression and the incorrect identification of splice variants. These contaminants can be introduced at various stages of the RNA-Seq process [23,24].

To address this, the TrAnnoScope pipeline uses Blobtools2 [25] to remove microbial and cross-contamination from long reads. Blobtools2 analyzes, visualizes, and filters assemblies based on the GC content, coverage, and taxonomic information. This tool is especially useful for de novo data, helping to identify and remove contaminants, thereby improving the quality of assemblies. The input files for Blobtools2 are generated using Bowtie2 [26] for coverage data, Benchmarking Universal Single-Copy Orthologs (BUSCO) for taxonomic classification [27], and BLAST [28] against the National Center for Biotechnology Information (NCBI) nucleotide (NT) database for taxonomic information.

2.2.4. Error Correction

Error correction can be beneficial for improving the accuracy of long reads, which typically have higher error rates than short reads [29]. In TrAnnoScope, the FMLRC tool is available as an optional step for correcting errors in long reads by leveraging complementary Illumina reads. FMLRC utilizes a multi-string Burrows–Wheeler transform and FM index to retrieve k-mer frequencies and construct de Bruijn graphs from short reads. It performs two passes with short and long k-mer values to correct unsupported regions in long reads, resulting in a more comprehensive correction process [30]. Owing to its efficiency and accuracy, FMLRC is a robust choice for error correction [29].

2.2.5. Clustering and Classification

To eliminate redundancy and reduce the complexity of the transcript data, TrAnnoScope employs two tools for clustering and classification: CD-HIT-Est [31], which is provided as an optional tool, and EvidentialGene [32]. CD-HIT-Est clusters transcripts based on their sequence similarity to remove redundancy within each sample, whereas EvidentialGene classifies the transcripts as primary and alternate forms based on their quality and potential function across the combined dataset.

2.2.6. Quality Assessment

TrAnnoScope includes a quality assessment step that utilizes several tools to evaluate the transcriptome comprehensively. NanoPlot [33] generates descriptive statistics, such as the mean, median, and N50 values, providing a clear overview of transcript continuity. BUSCO [27] further assesses transcriptome completeness against a user-defined lineage, ensuring the presence of elements for the organism of interest. Additionally, transcriptome quality is evaluated by comparing the number of FL or nearly FL transcripts against known protein databases, using an approach similar to the Trinity method of counting FL transcripts [34]. However, we implemented an in-house Bash script that calculates the percentage of high-coverage proteins present in the transcriptome compared to SwissProt by default or a user-defined custom database for closely related organisms. Together, these methods provide a detailed and accurate measure of transcriptome quality.

2.2.7. Annotation

Functional annotation of the transcriptome is performed using Trinotate [35], an advanced annotation suite created for the automated functional annotation of transcriptomes. Trinotate integrates multiple sequence databases, including Pfam [36], SwissProt [37], SignalP [38], TMHMM [39], EggNOG [40], and Infernal [41], to provide a comprehensive annotation. Additionally, TrAnnoScope provides homology searches against NCBI (non-redundant protein) NR and NT databases. Users can select the databases they wish to use for annotation in the configuration file. However, they must manually download and prepare the necessary files for the NT and NR databases prior to the annotation process. To facilitate this, Bash scripts to automate the downloading and indexing of the databases were provided.

To accelerate the annotation process, we implemented a strategy that splits the input files into user-defined chunks, enabling the parallel execution of homology searches. This approach significantly reduces the time required to annotate large datasets. After completing the homology searches, the results are parsed using modified Trinotate helper scripts to generate detailed annotation files enriched with Gene Ontology (GO) [42] from SwissProt, Pfam and EggNOG. For further functional insights, we implemented in-house R scripts to create plots, including GO, Kyoto Encyclopedia of Genes and Genomes (KEGG) [43], Eukaryotic Orthologous Groups (KOG) [44], and species distributions (based on the NCBI NR database).

2.2.8. Data Selection

To demonstrate the functionality and versatility of our pipeline, we processed publicly available RNA-Seq reads from the zebra finch (Taeniopygia guttata), an avian model used to study the neural mechanisms of local learning and social behavior [45]. This species is notable for its complex vocalizations, ease of breeding in captivity, and pronounced sexual dimorphism, making it a valuable model for understanding vocal learning and its implications for human speech and language development [45,46].

RNA-Seq data were obtained from the NCBI Sequence Read Archive [47] (SRR8551559, SRR8551563, SRR8551565, SRR8551567, SRR8551558, SRR8551562, SRR8551564, and SRR8551566). This dataset includes Illumina NextSeq 500 paired-end reads (2 × 76 bp) and long-read sequences generated via PacBio SMRT Sequel from various tissues, including the brain, ovary, and testis (Table S1). These complementary data types provided an ideal scenario for evaluating the ability of the pipeline to integrate and interpret both short and long reads for a comprehensive transcriptome assembly. TrAnnoScope was executed on the GenomeDK cluster using SLURM for all analysis steps. Detailed information about the outputs and runtime is provided in Supplementary File S3.

To fully utilize TrAnnoScope, users must provide both long and short reads. However, the pipeline is flexible: users can choose to run the quality control and preprocessing modules with only Illumina reads or execute the preprocessing, clustering and classification, quality assessment, and annotation modules with only long reads.

2.2.9. Mapping to the Zebra Finch Reference Genome

The final transcriptome generated by the TrAnnoScope pipeline was mapped to the current zebra finch reference genome (RefSeq: GCF_003957565.2) using minimap2 -ax splice –secondary=no -C5 parameters. Alignment statistics were obtained using the samtools flagstats option.

3. Results and Discussion

In this study, we applied TrAnnoScope to RNA-Seq data from the zebra finch (Taeniopygia guttata) to evaluate its effectiveness in transcriptome assembly and functional annotation. The dataset included Illumina reads and PacBio long-read sequences from brain, ovary, and testis tissues. Our primary objective was to assess the ability of the pipeline to integrate and interpret these complementary data types to achieve a comprehensive transcriptome assembly. This section details the results obtained, focusing on the key findings related to preprocessing, contamination removal, error correction, clustering and classification, quality assessment, and functional annotation.

In the preprocessing step for Illumina reads, we evaluated the read quality before and after data processing using FastQC, with the quality metrics compiled into a single report for concise visualization via MultiQC. Contaminants, including rRNA, mtRNA, and other potential contaminants, were removed using FastQScreen. This process involved hits against the LSU_Ref and SSU_Ref Silva databases v.138 [48], the zebra finch mitochondrial genome (NCBI Reference Sequence: NC_007897.1), and the FastQScreen database of vectors, adapters, and GRCm38 rRNA. For each sample, only minor hits were detected in the rRNA databases, primarily mitochondrial reads from the zebra finch (~5%); approximately 95% of the reads had no-hits (Figure S1). Reads that did not map to these databases (no-hits) were retained for downstream analyses. Adapter sequences and low-quality bases were trimmed using fastp. Table 1 presents the preprocessing statistics for each Illumina sample. Following preprocessing, the number of retained high-quality reads ranged from 27,446,223 to 33,201,568 per sample, ensuring robust data for the downstream assembly.

The initial preprocessing of the PacBio raw data (subreads) followed the Isoseq3 package for each sample. CCSs were generated using the default minimum number of subreads (default: 3). FL transcripts were identified, and primers were removed using lima with –peek-guess, while isoseq3 refine was used to remove poly-A tails and artificial concatemers to obtain FLNC reads. High-quality FL consensus sequences were obtained using the isoseq3 cluster2 with the –singletons parameter. By default, isoseq3 cluster2 retains isoforms that are represented by at least two FLNC reads. To capture rare but potentially significant isoforms, the --singletons option was employed for this analysis to include these single-read isoforms in the consensus sequences.

To eliminate potential contamination from PacBio reads, BlobTools2 was employed to retain only vertebrate sequences and other hits for further analysis. Except for the ovary_2 sample, which contained hits from the Annelia phylum, all other samples contained only vertebrate sequences and other hits (Figure S2). Additionally, rRNA and mitochondrial fragments were identified and removed by aligning the reads against the NCBI nucleotide database (NT, retrieved on 9 February 2024) using Blastn. High-quality, clean FL reads from each sample were corrected using FLMRC with default parameters. The error-correction step is provided as an option for users. With the advancements in sequencing technologies, the accuracy of long reads has been gradually increasing. The error rate of Nanopore sequencing has improved from approximately 64% for R7 to approximately 84–95% for R9.4 [49]. In contrast, the PacBio platform, utilizing the CCS approach, achieves a greater than 99% consensus accuracy [50,51]; however, systematic errors can still persist, especially in homopolymeric regions [52,53]. Depending on the accuracy of the sequencing platform utilized, users can proceed with downstream analysis without an additional error-correction step. In our study, we chose to perform the error correction because our analysis revealed that the BUSCO scores of the error-corrected clean reads were superior to those of the non-error-corrected clean reads (Figure S3). This finding aligns with the literature suggesting that error correction can still enhance the overall quality of transcriptomic data [51]. Table 2 shows the number of reads obtained at each preprocessing step. To remove redundancy, CD-HIT-Est was first employed for each sample (Table 2), and subsequently, EvidentialGene was used for the combined sample for the further classification of mRNA reads and predicted protein sequences (Table 3).

After EvidentialGene, 39,984 transcripts were obtained with an average transcript length of 3097.7 bp, a median of 2794 bp, and an N50 of 4108 bp, indicating the robustness of the assembly and the inclusion of long, high-quality transcripts (Table 3). Among the 39,984 predicted proteins, the average protein length was 398.2 amino acids with a median of 271 amino acids and an N50 of 597 amino acids. Notably, 86.7% of the transcripts were classified as complete proteins, further demonstrating the effectiveness of the pipeline for generating FL sequences. For FL representation analysis, a total of 26,141 transcripts matched the zebra finch protein sequences (GCF_003957565.2) with an e-value threshold of 1 × 10⁻²⁰. Of these, 24,560 transcripts presented as nearly FL (>70% coverage) relative to the zebra finch reference protein sequences. Among these, 10,575 transcripts were classified as FL transcripts with 100% coverage (Figure 2A, Table S2).

The BUSCO assessment results demonstrated the completeness of the transcriptome based on the presence of BUSCOs from the vertebrate lineage (Figure 2B). The analysis revealed that FL transcripts matched 83.8% of the single-copy orthologs in BUSCOs, comprising 79.1% complete orthologs and 4.7% fragmented orthologs out of a total of 3354 orthologs. Additionally, 16.2% of the orthologs were classified as missing (Table 3). While these BUSCO result indicates a reasonably high level of completeness, it also highlights potential limitations in the transcriptome assembly. The observed BUSCO completeness score can be attributed to the limited sampling of only three specific tissues. BUSCO scores assess the presence of conversed, single-copy orthologs expected in a comprehensive transcriptome. By focusing on a smaller number of tissues only a subset of actively expressed transcripts was captured, which reduces the overall BUSCO score. This result does not indicate poor data quality but reflects the targeted nature of the sampling strategy [54]. To improve the completeness of the transcriptome, future study could incorporate additional tissue samples to capture broader range of gene expression. Despite the limitations, the BUSCO results, combined with the descriptive statistics, underscore the ability of the pipeline to accurately obtain transcriptomes from Illumina and PacBio data, providing a reliable foundation for subsequent functional analyses.

The alignment of our transcriptome to current zebra finch genome yielded a high mapping rate of 99.63%, underscoring the accuracy and reliability of the TrAnnoScope pipeline for generating high-quality transcriptomic data. Despite this high mapping rate, a total of 154 transcripts did not map to the genome.

The identification of 154 transcripts that did not map to genome presents an intriguing opportunity for further exploration into the genomic landscape if the zebra finch. Among these unmapped transcripts, a significant portion (136 transcripts) exhibited notable hits in homology searches (see Supplementary File S2). It indicates that these transcripts likely represent real biological data rather than artifacts. Interestingly, 7 of these 136 transcripts aligned with sequences from previous zebra finch genome assemblies that are not present in the current genome version. This suggests that the unmapped transcripts may represent genomic regions that have been lost or altered in the latest assembly. This observation highlights the ongoing refinement of genomic resources and emphasizes the importance of considering multiple assembly versions in genomic data analysis. The rest of the unmapped 136 transcripts aligned with closely related species, particularly those within the passerine bird family. This indicates that these transcripts may possess functional relevance, potentially aligning with genes conserved across closely related passerine species and suggesting a shared heritage that may be critical for understanding evolutionary relationships in this group.

Among the remaining unmapped transcripts lacking homology search results, eight showed hits only for SignalP and TmHMM, indicating that they might encode peptides with specific targeting signals or transmembrane domains. This finding hints at their potential functional roles within cellular processes. However, the remaining ten transcripts, which yielded no significant information, raise questions about their biological significance. These transcripts could represent novel genes, warranting further exploration into their functions. Alternatively, they may be artifacts, emphasizing the need for additional validation.

The functional annotation of the assembled transcriptome was conducted using a comprehensive set of databases via Trinotate (e-value 1 × 10⁻⁵), providing valuable insights into the roles and characteristics of the predicted proteins (Table 4). The annotation revealed broad coverage across multiple databases, enhancing confidence in the functional assignments. Out of the 39,984 transcripts, 70.7% (28,274) had significant hits against the UniProt/SwissProt database using Blastx, while 62.7% (25,051) were confirmed through Blastp searches. Domain-based searches using Pfam identified conserved protein domains in 59.2% (23,689) of transcripts, highlighting their protein-coding potential.

GO terms from SwissProt, Pfam and EggNOG were assigned to 71.0% (28,399) of transcripts, categorizing them into biological processes, molecular functions, and cellular components (Figure 3). Additionally, KOG classifications were found for 66.0% (26,409) of transcripts, offering insights into their evolutionary relationships and potential functional roles (Figure 4). The KEGG pathway analysis annotated 65.3% (26,097) of transcripts, linking them to various metabolic and signaling pathways (Figure 5).

Structural and localization predictions identified transmembrane domains in 17.9% (7155) of transcripts, as predicted using TMHMM, and signal peptides in 5.8% (2305), as predicted using SignalP. A small portion of the transcripts (0.4%, 169) were annotated as non-coding RNAs using Infernal.

Further annotation against the NR database using Blastx produced hits for 82.7% (33,049) of transcripts, while Blastp hits were obtained for 70.8% (28,303). Nearly all transcripts (99.6%, 39,827) matched in the NT database through Blastn, indicating the strong conservation of these sequences against the known sequences. These extensive annotation results underscore the robustness of the transcriptome, with most transcripts being functionally annotated across various databases, providing a rich resource for downstream biological analyses.

The species distribution of the transcripts obtained from the TrAnnoScope pipeline underscores the effectiveness of our approach in functional annotation of zebra finch data (Figure 6). The Blastx homology search against the NR database revealed that a significant number of hits (17,121) were assigned to the zebra finch, indicating a strong alignment between our assembled transcripts and existing annotations. This result highlights the robustness of the TrAnnoScope pipeline in processing RNA-Seq data and achieving meaningful annotations that are crucial for downstream analyses. Furthermore, the presence of additional hits to closely related species, such as the society finch (Lonchura striata domestica), Gouldian finch (Chloebia gouldiae), and canary (Serinus canaria) can be attributed to high degree of genetic similarity among these species. This observation reinforces the notion that functional conservation is common among closely related species, facilitating the identification of homologous genes and conserved biological functions. Moreover, the identification of transcripts that align with other birds in the Passeriformes order, such as starlings (Lamprotornis superbus), swallows (Hirundo rustica rustica), sparrows (Melospiza melodia maxima, Passer montanus), and Réunion grey white-eyes (Zosterops borbonicus), reflects the shared ancestry within this diverse avian group.

The successful functional annotation of these transcripts not only validates the efficiency of the TrAnnoScope pipeline but also provides valuable insights into the evolutionary relationships among these avian species. The alignment of our data with that of the zebra finch and its relatives may contribute to a deeper understanding of the genetic basis of traits relevant to adaptation and survival in varying environments. Additionally, these findings can serve as a foundation for future studies aimed at exploring gene function, expression patterns, and evolutionary dynamics within the Passeriformes order and beyond.

While the TrAnnoScope pipeline has demonstrated its effectiveness in processing and annotating transcriptomic data, there are several areas for future improvement. One significant limitation is the current support for long-read preprocessing, which is restricted to PacBio reads. Incorporating preprocessing capabilities for ONT reads would greatly enhance the versatility of the pipeline, allowing users to leverage the strengths of both sequencing platforms for comprehensive transcriptome assembly. Additionally, providing support for differential expression analysis within the pipeline would facilitate more in-depth investigations into gene expression patterns across various conditions and tissues. This enhancement could empower researchers to derive meaningful biological insights from their data, further expanding the utility of TrAnnoScope in the field of transcriptomics. Addressing these gaps will not only improve the overall functionality of the pipeline but also increase its appeal to a broader range of users conducting diverse transcriptomic studies.

4. Conclusions

In this study, we introduced TrAnnoScope, a comprehensive pipeline for transcriptome analysis and annotation that utilizes both short- and long-read data. Applying TrAnnoScope to zebra finch RNA-Seq data demonstrated its capability to generate high-quality transcripts and functional annotations across multiple databases. The pipeline efficiently processes large datasets, from quality control to final annotation, resulting in a transcriptome with significant functional insights.

TrAnnoScope is built with Snakemake, and its modular design allows for easy customization while requiring minimal programming skills, making it accessible to users with varying levels of expertise. Its parallelized steps and user-defined parameters enhance the speed and reliability of transcriptome analysis. Although some manual database preparation is necessary, the pipeline remains a valuable tool for researchers, particularly those working with non-model organisms. Future updates will focus on automating database management and expanding the preprocessing options for other platforms.

Overall, TrAnnoScope is a versatile and efficient tool for transcriptomics, providing a robust platform for transcriptome analysis.