The battle for reads: evaluating strategies to tackle multi-mapping in RNA-seq quantification in highly repetitive genomes

The battle for reads: evaluating strategies to tackle multi-mapping in RNA-seq quantification in highly repetitive genomes | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article The battle for reads: evaluating strategies to tackle multi-mapping in RNA-seq quantification in highly repetitive genomes Aldana A Cepeda Dean, Virginia Balouz, Carlos A Buscaglia, Natalia Rego, and 1 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-7888056/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract Background: RNA sequencing (RNA-seq) enables transcript quantification and isoform analysis in diverse biological contexts, but accurately measuring expression from highly related genomic regions remains challenging. Multi-mapped reads—those aligning equally well to multiple loci—pose a major computational hurdle and compromise the overall accuracy of transcriptome resolution. Results: We herein evaluated five RNA-seq pipelines—Bowtie2 + featureCounts, STAR + featureCounts, STAR + Salmon, Salmon, and Kallisto—on their ability to quantify gene expression in Trypanosoma cruzi , a parasitic protozoan with a highly repetitive genome characterized by the abundance of large multigene families.Using real RNA-seq data, we first compared gene-level outputs, with emphasis on multigene family representation. Simulated transcriptomes were used to benchmark quantification accuracy under controlled conditions. Among the best-performing strategies (Salmon, Kallisto, and STAR + Salmon), we further tested whether including untranslated regions (UTRs) in gene annotations improved the assignment of ambiguous reads. Conclusions: Overall, the alignment-free transcriptome quantifiers Salmon and Kallisto achieved the most accurate performance, closely matching simulated values. Incorporating UTR annotations improved read assignment accuracy, particularly for STAR + Salmon. These tools not only enable global expression quantification but also facilitate precise read allocation between members of the same gene family, with up to 98% sequence identity. Our results highlight the critical role of annotation quality and quantification strategy in improving gene expression estimates. Multi-mapping RNA-seq Salmon STAR Bowtie2 Kallisto Trypanosoma cruzi multigene family Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 BACKGROUND RNA sequencing (RNA-seq) has become an essential tool in transcriptomics, enabling the accurate assessment of gene expression, discovery of novel transcripts, and analysis of transcript isoforms across diverse biological contexts. However, a persistent challenge in RNA-seq data analysis is the reliable quantification of transcripts derived from duplicated loci or repetitive genomic regions. In particular, the assignment of sequencing reads that map equally well to multiple loci–commonly called multi-mappers–poses a major computational hurdle ( 1 , 2 ). Multi-mappers can arise from various sources, including gene families with high sequence similarity among members, pseudogenes, repetitive sequences, and transposable elements. In organisms with highly repetitive genomes, this can account for a substantial fraction of reads–up to 40%, depending on the sample and sequencing strategy ( 1 ). The severity of the problem depends on the read length, the nature and extent of genomic repetition, and the degree of sequence identity among paralogous regions. Discarding multi-mapped reads, a common step in many standard pipelines, can lead to biased expression estimates, particularly for multigene families and repetitive elements ( 3 – 5 ). The issue is further exacerbated in single-cell RNA-seq studies, where short reads and sparse coverage amplify the uncertainty associated with ambiguous mappings ( 6 , 7 ). A related challenge arises in complex eukaryotes, where genes produce multiple isoforms that share most of their sequence, and in certain pathogens that evolved large multigenic families as part of their strategy to achieve persistent infection in the host. This extensive sequence overlap makes it difficult to assign reads uniquely to specific transcripts, thereby complicating transcript-level quantification ( 8 , 9 ). To address this issue, various computational strategies have been developed, encompassing both alignment and quantification steps. These methods differ in how they handle ambiguity and in the trade-offs they make between accuracy, interpretability, and computational efficiency. At the alignment step, tools like Bowtie and BWA–commonly used for Illumina data–perform fast, base-level mapping and are typically employed in organisms that lack introns, where spliced alignment is not required ( 10 , 11 ). In contrast, STAR is a splice-aware aligner that retains multiple mapping positions per read, providing greater flexibility for transcriptomic analyses, similar to other aligners such as HISAT2 and Subread ( 12 , 13 ). At the quantification step, traditional pipelines often rely on tools such as featureCounts to generate gene-level counts from aligned reads ( 14 ). However, these approaches may not fully account for the uncertainty introduced by multi-mapping reads. More recent methods, such as Salmon and Kallisto, employ alignment-free strategies based on k-mer indexing and probabilistic models ( 9 , 15 ). Designed primarily for transcriptome quantification, these tools bypass the need for complete base-level alignment and incorporate statistical models to assign reads to transcripts while explicitly modeling mapping uncertainty. Their speed and ability to handle repetitive or paralogous sequences have made them increasingly popular, especially for large-scale or complex transcriptomes ( 1 , 2 ). Despite these advances, the challenge of accurately assigning multi-mapped reads remains. An underexplored issue concerns the alignment and gene quantification in eukaryotes whose genomes are not intron-rich but instead characterized by extensive repetition and structural complexity. In such contexts, where a large fraction of transcripts originates from highly similar loci, even the most advanced quantification tools may struggle to resolve ambiguous mappings. Trypanosoma cruzi , the etiological agent of Chagas disease, exemplifies this problem. More than 50% of its genome consists of repetitive sequences, which comprise a wide variety of transposable elements and large multigene families with hundreds of copies, including functional genes and pseudogenes. Prominent among these are the trans -sialidases (TS, ~ 1400 genes), mucins (~ 900 genes), mucin-associated surface proteins (MASP, ~ 800 genes), dispersed gene family-1 (DGF-1, ~ 80 genes) and GP63 proteases (~ 200 genes) ( 16 – 20 ). These gene families play critical roles in host-parasite interactions, immune evasion, and virulence; and are characterized by high sequence similarity, extensive though usually focalized allelic diversity, coordinated expression, and remarkable structural plasticity—all of which compound the difficulty of resolving transcript-level expression from RNA-seq data ( 21 , 22 ). In this study, we use T. cruzi as a model to evaluate the performance of widely used RNA-seq analysis strategies under the most challenging conditions. We first analyze RNA-seq data, comparing outputs across pipelines with special focus on parasite multigenic families. We then simulate transcriptomic datasets to further assess the capacity of the best-performing tools to resolve ambiguous read assignments. Finally, we explore whether incorporating improved gene annotations, such as including untranslated regions (UTRs), can enhance the assignment of multi-mapped reads. Through this systematic approach, we aim to provide practical guidance for RNA-seq quantification in T. cruzi and other organisms with challenging genomic architectures. RESULTS 1. Distribution of Multi-Mapping Reads in T. cruzi . Multi-mapping reads—those aligning equally well to multiple annotated features—represent a major obstacle in RNA-seq quantification, particularly in genomes with high sequence redundancy. Trypanosoma cruzi , characterized by an abundance of large multigene families and widespread repetitiveness, provides a compelling model to examine how different computational strategies handle this challenge. To quantify the extent of multi-mapping in T. cruzi , we analyzed RNA-seq data using three representative tools: the aligners Bowtie2 and STAR, and the alignment-free quantifier Salmon. STAR and Salmon are among the most widely adopted in current RNA-seq workflows, while Bowtie2, though less common today, remains frequently used in studies of T. cruzi and other intron-poor parasitic genomes (Additional file 1: Table.S1) ( 23 – 27 ). The distribution of multi-mapping percentages across genes revealed a strikingly bimodal pattern for all three methods (Fig. 1 a). One mode centered near 0% reflects genes covered almost exclusively by uniquely mapped reads, while a second mode—ranging from ~ 60% to ~ 100%—corresponds to genes with high levels of ambiguous alignments. However, the shape and intensity of these modes differed markedly across methods, reflecting differences in how each tool defines and processes multi-mapping. Salmon and STAR both showed a dominant peak at 0%, representing ~ 42% and ~ 36% of all genes, respectively, and a sharp secondary peak between 96–97%, encompassing ~ 20% of the genes. In contrast, Bowtie2 showed a flatter distribution with far fewer genes supported solely by uniquely mapped reads, and a broader peak centered around 82%, involving ~ 54% of the genes. These results underscore the influence of mapping strategy on the detection and distribution of multi-mapping reads. A summary of total multi-mapping content per method is shown in Figure S1 , while detailed mapping statistics (including the percentage of aligned reads) are presented in Additional file 1: Table.S2. For comparative purposes, we performed similar analyses on evolutionarily related pathogens such as the trypanosomatids Trypanosoma brucei and Leishmania major , and the apicomplexans Toxoplasma gondii and Plasmodium falciparum (Additional file 2: Fig. S2 ). T. brucei exhibited an overall level of multi-mapping comparable to that of T. cruzi , but with a more sharply bimodal distribution: most genes showed either 0% or nearly 100% multi-mapping, with almost no genes falling in between. In contrast, T. cruzi included a subset of genes with intermediate multi-mapping levels, resulting in a less polarized distribution. In the same line, L. major showed a similar bimodal pattern, although with generally lower levels of ambiguity. In contrast, Toxoplasma and Plasmodium presented minimal multi-mapping, with nearly all reads mapping unambiguously. To further understand how different strategies quantify ambiguous reads in T. cruzi , we focused on five of its most expanded multigene families (MASP, Mucins, GP63, TS, and DGF-1). Each one of these families comprise between 100 and 1,000 copies showing varying extent of polymorphism (Additional file 1: Table.S3 and Additional file 2: Fig.S3), and collectively account for over 34% of multi-mapping reads in this organism (Fig. 1 b). Bowtie2 yielded the highest overall density of multi-mapping, with a sharp peak above 60% across all families, reflecting its tendency to report multiple alignments without probabilistic resolution (Fig. 1 c). In contrast, STAR and Salmon produced broader, more gradual distributions, suggesting a more nuanced handling of ambiguous reads (Fig. 1 c and 1 d). Interestingly, distribution profiles varied between multigene families and the methods used. TS, for instance, consistently showed a bimodal pattern across strategies, whereas GP63 and mucins exhibited such bimodality only with STAR and Salmon (Fig. 1 c and 1 d). The GP63 family, in particular, showed more elevated levels of multi-mapping (Additional file 1: Table.S3), consistent with its internal structure (Additional file 2: Fig.S3), which includes subgroups of nearly identical genes that complicate read assignment (Berná et al., 2025). In contrast, DGF-1 exhibited markedly lower levels of multi-mapping, especially with STAR and Salmon, likely due to greater sequence divergence among its members (Fig. 1 c and 1 d). Together, these results reveal the high prevalence of multi-mapping reads in T. cruzi , and their biased distribution across multi-copy gene families. They also demonstrate that the strategy used for the genomic anchoring of sequencing reads significantly impacts on the definition of the multi-mapping landscape in this parasite. 2. Impact of Multi-Mapping on T. cruzi Gene Expression Quantification To evaluate how mapping and quantification strategies affect transcript abundance estimates in T. cruzi , we applied a standardized analysis pipeline to both experimental and simulated RNA-seq datasets (Fig. 2 ). We tested five widely used tool combinations—Bowtie2 + featureCounts, STAR + featureCounts, STAR + Salmon, Salmon, and Kallisto—to assess their performance in quantifying gene expression, particularly for highly redundant gene families. We first compared gene-level read counts across the three alignment-based strategies using experimental RNA-seq data. Overall, ~ 99% of T. cruzi genes received fewer than 1,000 reads across pipelines (Fig. 3 a), and highly expressed genes (> 1,000 reads) were consistently detected. No correlation was observed between expression level and multi-mapping percentage, although increased variability—consistent with stochastic noise—was evident among lowly expressed genes (Additional file 2: Fig.S4). When focusing on genes with up to 300 reads (covering ~ 98% of the dataset), differences between pipelines became more apparent. STAR + Salmon and Bowtie2 + featureCounts yielded similar distributions (Fig. 3 b), whereas STAR + featureCounts resulted in markedly lower counts for multi-mapping-prone genes. Interestingly, featureCounts assigned multi-mapping reads when combined with Bowtie2 but not when combined with STAR. This highlights how differences in the aligner output can affect read assignment and downstream quantification. As shown, featureCounts relies heavily on the flags and tags present in the BAM file, which are aligner-dependent, and although the performance of this tool can be modified by adjusting alignment parameters, this is rarely done in practice. 3. Performance Assessment Using Simulated Data To assess the ability of each quantification strategy to recover actual gene expression levels, we simulated RNA-seq datasets based on predefined expression profiles. Briefly, gene-level count values were specified manually from T. cruzi experimental data and used as input for Polyester, which generated synthetic RNA-seq reads according to these target values. This approach allowed us to establish a known ground truth and systematically evaluate the performance of each pipeline in approximating simulated expression levels. To verify the stability of the simulation process, we compared the expression estimates across three replicates generated by Polyester. As no variability was observed between replicates for the three tested pipelines (Additional file 2: Fig.S5), we show a single replicate for all subsequent analyses. Genome-wide analyses revealed that alignment-free transcript quantification approaches (Salmon and Kallisto) outperformed traditional alignment-based methods, producing expression estimates closer to the expected simulated values, with Kallisto showing the lowest overall error metrics (Fig. 4 ). Among the aligner-based strategies, both STAR + featureCounts and Bowtie2 + featureCounts tended to overestimate transcript abundances, particularly for ‘highly-expressed’ genes, whereas STAR + Salmon displayed a more consistent pattern, closer to that of the alignment-free tools, albeit with some dispersion (Fig. 4 ). These differences were further reflected when we assessed the number of genes with perfect quantification (i.e., identical observed and expected counts). As shown in Additional file 2: Fig.S6, Salmon achieved perfect counts for 8,417 genes (47% of the T. cruzi genome), Kallisto for 7,180 (40%), while STAR + Salmon correctly quantified as low as 243 genes (1.4%). Quite similar results were observed when we focused the analysis on parasite multigene families (Additional file 2: Fig.S7) hence reinforcing the robustness of these observations. Overall, these data indicate that RNA-seq analysis using alignment-free transcript quantification approaches (Salmon, Kallisto, and STAR + Salmon) provides for a more accurate gene expression quantification in T. cruzi . 4. Assessing the Individual Expression of Multigene Family Members: The Battle Royale for the Reads. Depending on their mode of evolution, multigene families in T. cruzi differ in their internal sequence similarity, which may range from moderate to near-identical. This variation, illustrated in Additional file 2: Fig.S3, affects the resolution capacity of mapping tools, as high similarity can hinder the unambiguous assignment of reads. To examine how the best-performing strategies (STAR + Salmon, Salmon and Kallisto) handled this issue, we designed a simulation-based analysis focused on TS, GP63, mucins, and MASP genes. Specifically, we compared the expected and observed read counts for selected subsets of these gene families comprising genes sharing > 95% sequence identity among them. In all cases, and consistent with our previous results, alignment-free methods such as Salmon and Kallisto performed better than the alignment-based STAR + Salmon approach (Fig. 5). This trend translated into much lower MedAE values for Salmon and Kallisto as compared to STAR + Salmon. As expected, MedAE values for each gene family and quantification method increased as sequence identity approached 100% (Fig. 5). To further explore the resolution capacity of each method, we designed a complementary simulated scenario in which overexpression was artificially introduced in one gene of selected paralogous subgroups. In this setup, a single gene within a group of highly similar paralogs was assigned with 1,000 reads, while the remaining members retained their original simulated expression levels. This allowed us to test the performance of the quantification methods under the most stringent conditions. (Additional file 2: Fig.S8). In the TS family (97.3–97.75% identity), all strategies accurately recovered the simulated overexpression. For GP63 genes (99.23–99.65% identity), however, Salmon and Kallisto captured the overexpression reliably, while STAR + Salmon tended to distribute part of the signal across other family members. In the Mucin family (98.9–99.39% identity), only Salmon approximated the expected values, whereas the other methods largely failed to detect the overexpressed gene (Additional file 2: Fig.S8). The MASP family, which included a pair of genes with 100% identity, presented a limit case. Kallisto split the reads nearly equally between both copies, consistent with its expectation-maximization framework (EM-based). In contrast, Salmon and STAR + Salmon assigned all reads to a single gene, completely ignoring the identical paralog (Additional file 2: Fig.S8). These results underscore the fundamental challenge of quantifying genes or transcripts with identical or near-identical sequences, where some methods retain multi-mapping uncertainty while others introduce assignment bias by artificially resolving the ambiguity. More importantly, these findings highlight the superior accuracy and resolution of alignment-free methods in quantifying expression in multigene families with high sequence similarity. 5. Effect of Including UTRs in the Quantification Building on the three selected strategies, we further explored the impact of transcript annotation quality on quantification accuracy. To that end, we compared simulated and observed read counts under two gene annotation schemes: one including only the CDS, and another one incorporating both the CDS and the UTRs (Fig. 2 ). Because T. cruzi genome annotations typically lack UTRs, we predicted these regions using the UTRme tool ( 28 ), based on RNA-Seq data from the YC6 strain. This exercise yielded a mean length of ~ 200 bp and a median length of 82 bp for 5′ UTRs and a mean length of ~ 500 bp and a median of ~ 300 bp for 3′ UTRs (Additional file 2: Fig.S9), thus consistent with previous reports ( 29 ). The resulting UTRs were integrated into the original annotation to generate an extended GFF file. As shown, the inclusion of UTRs consistently improved concordance between observed and expected expression levels across all methods (Fig. 6 and Additional file 2: Fig.S10a-b). Both the RMSE and the MAE decreased when extended annotations were used, underscoring the importance of complete transcript models for accurate quantification. For STAR + Salmon (Fig. 6 , upper panel), CDS-only annotations resulted in substantial dispersion around the identity line, with both over- and under-estimation of gene expression levels. In contrast, annotations including UTRs led to a tighter distribution of observed and expected values and thereby to lower error metrics. A similar trend was observed for Salmon alone (Fig. 6 , middle panel). However, it should be noted that although Salmon more closely followed expected values than STAR + Salmon for many genes (especially those with moderate expression), it yielded a higher overall RMSE (515.39 vs. 406.21). This increase was driven by a small number of multi-mapping-rich transcripts, where Salmon assigned all reads to a single representative transcript, leaving others with zero or near-zero counts (Additional file 2: Fig.S10 c). Most importantly, Kallisto (Fig. 6 , bottom panel) yielded the best overall performance in terms of both RMSE and MAE. Its estimates closely followed the identity line with minimal dispersion, indicating high accuracy. The positive effect of UTR inclusion was especially marked, further underscoring the value of full-length transcript definitions. Together, these results confirm that extended annotations improve the accuracy and consistency of gene expression estimates, particularly when using transcript-level quantification tools. DISCUSSION Accurate RNA-seq quantification is largely hindered by multi-mapping reads, especially in genomes with high redundancy. In mammals, multi-mapping affects ~ 5–10% of reads ( 3 ), and can reach up to 30% depending on factors such as read length, sequencing layout, and the mapping algorithm used ( 30 ). In T. cruzi , this issue is greatly amplified: depending on the tool, between 50–80% of reads map ambiguously. A major contributor to this phenomenon is the outstanding expansion of multigene families coding for virulence factors that play critical roles in parasite infection, immune evasion, and pathogenesis ( 21 , 22 ). Although multigene families represent ~ 14% of T. cruzi annotated genes, they account for over 34% of multi-mapped reads in our dataset—consistent with their high copy number and sequence similarity. Moreover, the contribution of these gene families to the overall T. cruzi multi-mapping is most likely higher than our gene-based estimations. As extensively reported, up to 80% of the total sequence dosage for these gene families is composed by transcriptionally active pseudogenes ( 21 , 22 ). In addition to multi-copy gene families, the T. cruzi genome is enriched in other repetitive sequences such as transposons and satellites that may further complicate accurate transcript expression quantification. Compared to other unicellular parasites, T. cruzi displays a uniquely complex multi-mapping landscape, highlighting how genome architecture shapes mapping ambiguity and positioning this organism as a stringent benchmark for RNA-seq quantification tools. Our comparative evaluation shows that alignment-free, transcriptome-aware quantifiers—Salmon ( 9 ) and Kallisto ( 15 )— provide the most accurate expression estimates in T. cruzi . These methods probabilistically assign reads, outperforming alignment-based workflows such as Bowtie2 + featureCounts, which has been the most widely implemented strategy in T. cruzi RNA-seq studies. Notably, they can resolve reads between paralogs sharing up to ~ 98% sequence identity, although perfect or near-perfect duplicates remain unresolved. For use cases requiring alignment (e.g., visualization, variant calling), STAR combined with Salmon quantification offers a practical compromise with accuracy comparable to alignment-free methods. Annotation quality critically affects multi-mapping resolution. Incomplete gene models—particularly those lacking UTRs—inflate ambiguity by omitting sequence features that could distinguish otherwise similar transcripts. In T. cruzi , where UTRs are not included in standard annotations, extending gene models substantially improved quantification accuracy ( 31 , 32 ). Salmon and STAR + Salmon, in particular, better approximated the simulated values when extended annotations were used, underscoring their sensitivity to UTR-derived sequence diversity. These results are consistent with benchmarking studies showing that annotation and quantification are among the primary sources of variation in RNA-seq workflows ( 33 – 36 ). Kallisto also consistently achieved low MAE and RMSE across all conditions, though it benefited less from annotation refinement. This may reflect its robust transcript-level quantification strategy, which avoids zero estimates even when transcripts are identical. In contrast, Salmon tends to assign all reads to a single transcript when encountering exact duplicates, effectively silencing the others. Together, these findings underscore the central role of both quantification strategy and annotation quality in shaping RNA-seq outcomes. Beyond our benchmarking, a broader range of methods has been proposed to explicitly address multi-mapping. Some tools aim to redistribute rather than discard ambiguous reads ( 3 , 37 ). MMquant ( 38 ), for example, assigns multi-mapped reads proportionally across all compatible features. This maximizes data usage and avoids underestimating repetitive gene families, but may overinflate expression in the presence of pseudogenes or paralogs. Other methods propose alternative allocation strategies: ShortStack ( 39 ), originally designed for small RNA-seq, assigns reads to the locus with the highest local density of uniquely mapped reads, an approach that can improve resolution in highly redundant regions, though it is limited to single-end data. Allo ( 40 ) a method developed for ChIP-seq, models the ambiguity explicitly using allocation likelihoods, and could potentially be adapted for RNA-seq contexts involving repetitive elements. MGcount offers a graph-based approach to quantify ambiguous reads across overlapping features and diverse RNA biotypes ( 41 ); while developed primarily for total-RNA and small RNA data, its underlying logic may hold promise for complex transcriptomes such as T. cruzi .Despite their potential, these tools remain underutilized, due to their specific input requirements or implementation complexity in typical RNA-seq workflows. Long-read RNA-seq technologies offer a promising alternative for resolving isoforms and paralogous transcripts that remain indistinguishable with short-read data. Despite technical challenges such as base-calling errors and incomplete transcript coverage, recent studies have shown that several dedicated tools perform well across different long-read platforms. Although short-read methods still outperform long-read-based approaches in terms of quantification accuracy under most conditions, combining both technologies may enhance transcript quantification in complex transcriptomes. As tools and protocols continue to improve, long-read RNA-seq is expected to provide increasingly accurate abundance estimates in the near future ( 42 , 43 ). Single-cell RNA-seq adds further complexity. Many pipelines, such as Cell Ranger ( 44 ), discard multi-mapped reads by default to avoid false positives. This conservative approach can underrepresent highly repetitive gene families. In contrast, tools like STARsolo ( 45 ), Salmon-Alevin ( 46 ) and Kallisto|bustools ( 47 ) use probabilistic mapping strategies that better handle ambiguity and reduce memory usage. Benchmarking studies ( 36 , 45 ) have shown that the choice of quantification method significantly impacts marker gene detection and cell-type classification. Our findings in T. cruzi mirror these observations. Conservative pipelines minimize false positives but risk overlooking biologically relevant signals, while inclusive strategies may inflate expression in repetitive regions. This trade-off is evident in recent single-cell studies. The T. cruzi cell atlas ( 32 ) employed Cell Ranger, likely underestimating expression of surface protein families. In contrast, Inchausti et al. (2025) used Kallisto|bustools with EM-based redistribution and empirically defined 3′ UTR annotations, enabling better assignment of reads and resolution of paralogous transcripts. These examples reinforce how analytical choices around multi-mapping and annotation directly shape the detection of gene family expression in T. cruzi. Another emerging challenge is posed by the increasing availability of haplotype-resolved assemblies ( 48 , 49 ), as reads may map equally well to alternative haplotypes or to true paralogs, confounding quantification. Current linear-reference genome approaches cannot resolve these cases. Although haplotype-aware tools like STARconsensus offer partial solutions, systematic allele-aware quantification methods are still in development. Graph-based reference genomes and pangenomes ( 50 , 51 ) represent a promising shift. By modeling genetic variation directly, these approaches may improve read assignment and reduce reference bias. Their integration into RNA-seq workflows remains limited, but for organisms like T. cruzi , which combine repeat-rich genomes and high haplotype divergence, they could prove transformative. Our results argue against the practice of masking repetitive regions in transcriptomic analyses. With improved annotations and modern quantification methods, expression from these regions can be recovered with reasonable accuracy—unlocking biologically meaningful information. At the same time, our study highlights remaining limitations: recent duplicates remain indistinguishable, annotations are still incomplete, and systematic biases persist depending on the aligner–quantifier combination used. Advancing toward accurate transcript quantification in complex genomes will require continued development of haplotype-aware algorithms and integration of pangenomic frameworks. In this context, T. cruzi emerges not only as a pathogen of biomedical relevance but also as a valuable system for testing the limits of current RNA-seq methods. Our study complements earlier benchmarking work by Love and Soneson ( 35 , 52 ) extending their findings to a genome characterized by extreme redundancy and structural plasticity. As transcriptomics moves into the pangenome era, organisms like T. cruzi offer critical opportunities to stress-test the next generation of analytical tools. CONCLUSIONS Our study demonstrates that T. cruzi provides a stringent benchmark for evaluating RNA-seq quantification in repetitive genomes. Multi-mapping reads, which affect only ~ 5–10% of polyadenylated transcripts in mammals, impact over half of T. cruzi genes and are disproportionately concentrated in large multigene families. Among tested strategies, alignment-free methods (Salmon, Kallisto) most closely recapitulated simulated values, while STAR + Salmon also performed robustly when supported by extended annotations. Notably, the inclusion of UTRs substantially improved quantification across pipelines, underscoring annotation quality as a key determinant of accuracy. Beyond T. cruzi , these findings highlight the broader need for probabilistic, annotation-rich, and ultimately haplotype-aware frameworks to resolve expression in complex transcriptomes. As graph-based and pangenome approaches mature, parasites with highly repetitive genomes such as T. cruzi will serve as critical models for testing and refining next-generation RNA-seq quantification algorithms. METHODS This study focuses on evaluating the impact of multi-mapping in RNA-seq analyses within the context of Illumina sequencing, which remains the predominant platform for quantifying gene expression. Accordingly, the described methodologies are either explicitly tailored to, or inherently compatible with, the characteristics of Illumina short-read data. Datasets and Quality Control RNA-Seq data from epimastigote and trypomastigote forms of T. cruzi (YC6 strain) were obtained from Li et al. (2016) and are available in the NCBI SRA under projects PRJNA251582 and PRJNA251583. Quality assessment of each dataset was performed using FastQC ( 53 ), version 0.11.9. Raw fastq files were trimmed using Sickle ( 54 ) with parameters: -t Illumina -q 30 -l 70, removing low-quality bases (Phred score < 30) and discarding reads shorter than 70 nucleotides. This ensured high-quality data for downstream analysis. Quantification of Multi-Mapping Reads Across Pathogen Genomes To evaluate the extent of multi-mapping across different eukaryotic pathogens, we used publicly available paired-end RNA-seq datasets, which were aligned to the corresponding reference genomes and annotations. Genome assemblies, gene annotations (GFF), and RNA-seq data were obtained for Trypanosoma brucei TREU927 (TriTrypDB release 68; ERR13925148), Leishmania major Friedlin2021 (TriTrypDB release 68; ERR2604479), Toxoplasma gondii RH (ToxoDB release 68; SRR28421832) (Gajria et al., 2008), and Plasmodium falciparum 3D7 (PlasmoDB GCA_000002765; SRR3455777). For T. cruzi , we used the YC6 genome assembly and RNA-seq data as described previously. Reads were aligned using STAR (v2.7.11b) with up to 50 allowed mappings per read (--outFilterMultimapNmax 50). Genome indices were built with the respective genome assemblies and annotations. BAM files were then separated into uniquely and multi-mapped reads using SAMtools (v1.15.1), and gene-level quantification was performed with featureCounts (v2.0.1). RNA-seq Read Alignment and Quantification Pipelines To evaluate the impact of mapping and quantification strategies on gene expression estimates in T. cruzi , we implemented a standardized RNA-seq analysis pipeline integrating both experimental and simulated data. Paired-end RNA-seq libraries from T. cruzi YC6 strain ( 4 ) were used as experimental input. Simulated datasets were generated based on these profiles and included scenarios with altered expression levels and annotation refinements. Overall, we tested five widely used tool combinations: Bowtie2 (v.2.5.3) + featureCounts (v.2.0.3) STAR (v.2.7.11) + featureCounts STAR + Salmon (v.1.10.1, alignment-based mode) Salmon (v.1.10.1, quasi-mapping mode) Kallisto (v.0.50.1, pseudo-alignment mode) In all cases, reads were mapped to the T. cruzi YC6 reference genome and transcriptome ( 55 ). For the Bowtie2 pipeline, reads were aligned in end-to-end mode, using the --sensitive preset, which balances speed and alignment accuracy in the presence of mismatches or indels. Additional options included --no-unal to exclude unaligned reads from the output. Aligned reads were quantified at the gene level using featureCounts. Transcript-level annotation was enabled via the -g transcript_id option and a modified GTF file including the UTRs (see below). Multi-mapping reads were not discarded (-M) and included in the quantification using --fraction, a minimum fractional overlap to the annotation feature of 10% was required (--fracOverlap 0.1). For the STAR pipeline, reads were aligned, allowing up to 50 multiple mappings per read. To tailor mammalian-trained STAR splice junction detection functions to T. cruzi genome, the following options were used: --outSJfilterReads Unique, --outSJfilterOverhangMin 150 150 150 150, and --outFilterType BySJout. Mismatch filtering was set with --outFilterMismatchNoverReadLmax 0.01, and reads were aligned in end-to-end mode (--alignEndsType EndToEnd). The output included sorted BAM files (--outSAMtype BAM SortedByCoordinate) and transcriptome alignments (--quantMode TranscriptomeSAM). These alignments were either used directly with featureCounts or passed to Salmon for quantification in alignment-based mode. For Salmon in quasi-mapping mode, transcript quantification was performed using the genome and GTF annotation (both indexed) with default parameters. Paired-end reads were processed with automatic library type detection (-l A), and the --validateMappings option was enabled to improve mapping accuracy. Finally, Kallisto was used to quantify transcript abundance via pseudo-alignment with default parameters, using the same input fasta file with the T. cruzi transcripts. To calculate the number of mapped reads for each methodology, we employed samtools flagstat (v.1.20), which provides global alignment statistics from the corresponding BAM files. Gene-level count matrices were generated for all strategies. Genes were classified into major multigene families (MASP, mucins, TS, GP63) or "Other" for comparative analysis, particularly focusing on the impact of multi-mapping reads on quantification accuracy. Genes belonging to multigene families were retrieved from CDS and GFF annotations. Pairwise sequence identities were calculated using the Biopython pairwise2 module. Global alignments were performed with the globalxx scoring scheme, and sequence identity was defined as the proportion of identical positions relative to the maximum sequence length. UTR Refinement of T. cruzi Gene Annotation The gene annotation file for T. cruzi YC6 (TriTrypDB release 68) was refined to include UTRs, which are not present in the standard annotation. UTRs were predicted using UTRme ( 28 ) with paired-end mode and default parameters, based on YC6 RNA-seq data. The highest-scoring UTRs were assigned to their corresponding genes using custom Python scripts and merged into the original GFF file. Multi-Mapping Analysis We estimated the proportion of multi-mapping reads for each mapping strategy. Only STAR, Bowtie2, and Salmon were included in this analysis, as all generate BAM or SAM files that enable tracking of read alignments. Multi-mapping reads were identified based on alignment flags and mapping quality scores. For STAR, secondary alignments were used to identify multi-mapped reads. For Bowtie2, the presence of the XS:i: tag was used to infer alternative mapping locations. Unique and multi-mapped reads were counted separately using a feature-based counting approach. Final multi-mapping percentages were calculated as the ratio of multi-mapped reads to total mapped reads, both globally and at the gene level. For Salmon, transcript-level alignments were obtained using the --writeMappings option to produce SAM files. These files were filtered to remove unmapped reads, and transcript-level counts were extracted by summarizing mappings per transcript using custom awk commands in bash. Kernel density estimation (KDE) was applied to the distribution of multi-mapping percentages for each strategy. Peaks in the KDE curves were identified as local maxima, and the area under the curve (AUC) corresponding to each peak was calculated using trapezoidal numerical integration to assess their relative contribution. All downstream processing, statistical analyses, and visualizations were performed in Python using the pandas, scipy, and numpy libraries for data manipulation and numerical operations, and seaborn and matplotlib for plotting and figure customization. Benchmark of Quantification Strategies To compare the impact of alignment and quantification methods on gene expression profiles, we evaluated the number of detected genes, transcript per million (TPM) distributions, expression correlations, and gene family representation across pipelines. Additionally, a public GitHub repository (Benchmarking-RNAseq, https://github.com/AldanaCepedaDean/Benchmarking-RNAseq ) was created to provide an overview of the methodologies discussed, together with example scripts and workflows to facilitate reproducibility and practical implementation. Simulation of RNA-seq Data Simulated RNA-seq reads were derived from real data for T. cruzi (strain YC6, trypomastigote stage; SRR1346054). Baseline expression values were obtained by averaging gene-level counts across three quantification strategies: Bowtie2 + featureCounts, STAR + featureCounts, and STAR + Salmon. These average counts were used to generate the reference expression table for simulations. Three simulation scenarios were implemented: Reference simulation : Reads were simulated according to the averaged expression values described above. Multigene family simulation : One representative gene was manually assigned with 1,000 reads, while the remaining family members retained their original expression values from the reference count table. This level of upregulation (from baseline expression levels ranging between 200 and 500 reads) was chosen to simulate a moderate and biologically plausible increase in transcript abundance. UTR-expanded simulation : Same as the reference simulation, but using transcript annotations that include both CDS and UTRs during indexing and alignment. All simulations were carried out using Polyester ( 56 ), generating 100-nucleotide paired-end reads with a 0.5% sequencing error rate, based on the error rate obtained from the YC6 reference transcriptome ( 4 ). Root Mean Square Error (RMSE) and Mean Absolute Error (MAE) were calculated to compare between simulated and observed values. For the simulation of multigene families, to assess the accuracy of expression estimates, we calculated the Median Absolute Error (MedAE) defined as: $$\:MedAE\:=\:median({|Observed}_{1}-{Expected}_{1}|,\:...,\:{|Observed}_{n}-{Expected}_{n}\left|\right)$$ where Observed n denotes the expression estimate from each quantification method for gene n , and Expected n ​ denotes the corresponding ground-truth value from the simulation. All computations were performed in Python using the NumPy and scikit-learn libraries for metric calculations. Data handling, statistical analysis, and visualization were conducted with pandas, numpy, seaborn, and matplotlib. Abbreviations RNA-seq, RNA sequencing; UTRs, Untranslated Regions; TS, trans -sialidase from T. cruzi ; MASP, mucin-associated surface protein from T. cruzi ; DGF-1, Dispersed Gene Family-1 from T. cruzi Declarations ACKNOWLEDGMENTS We are grateful to Lucas Inchausti for his help with running simulated data and for his kindness, and to Carlos Robello for his involvement in coordinating and facilitating Aldana A. Cepeda-Dean’s stay during her UNU-BIOLAC scholarship. FUNDING AACD holds a CONICET fellowship, whereas CAB and VB are career investigators from the same institution. AACD also received a CIN (Consejo Interuniversitario Nacional, Argentina) fellowship and a UNU-Biolac fellowship. LB and NR are PEDECIBA (Programa de Desarrollo de Ciencias Básicas, Uruguay) researchers and members of the Sistema Nacional de Investigadores (SNI, ANII, Uruguay). This investigation received financial support from the ANPCyT (PICT-2021-0284 to CAB) and by the Institut Pasteur de Montevideo and FOCEM - Fondo para la Convergencia Estructural del Mercosur (COF 03/11). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. AUTHOR’S CONTRIBUTION Conceptualization: AACD, LB, NR, VB, CAB. Design: AACD, LB, NR. Formal analysis: AACD, LB, NR. Methodology: AACD, LB, NR. Funding acquisition and resources: LB, NR, CAB, VB, AACD. Writing - original draft: AACD, LB, NR. Writing - review and editing: AACD, VB, CAB, LB, NR. Final approval of the version to be submitted: All the co-authors. All authors read and approved the final manuscript. ETHICS APPROVAL AND CONSENT TO PARTICIPATE Not applicable. CONSENT TO PUBLICATION Not applicable. COMPETING INTEREST The authors declare no competing interests. AVAILABILITY OF DATA AND MATERIALS. All data generated and analysed during this study are included in this published article and its supplementary information files. Additionally, a public GitHub repository has been created to provide an overview of the RNA-seq methodologies discussed in this review, along with example scripts and workflows to facilitate reproducibility and practical use by the community: Project name: Benchmarking-RNAseq Project home page: https://github.com/AldanaCepedaDean/Benchmarking-RNAseq Operating system(s): Linux, MacOS, Windows (with WSL) Programming language: Shell, Bash, R Other requirements: None This repository is intended as a resource for researchers to explore and implement the discussed strategies for handling multi-mapping reads in RNA-seq experiments. References Deschamps-Francoeur G, Simoneau J, Scott MS. Handling multi-mapped reads in RNA-seq. Comput Struct Biotechnol J. 2020 Jan 1;18:1569–76. Srivastava A, Malik L, Sarkar H, Zakeri M, Almodaresi F, Soneson C, et al. Alignment and mapping methodology influence transcript abundance estimation. Genome Biol. 2020 Sept 7;21(1):239. Almeida da Paz M, Warger S, Taher L. Disregarding multimappers leads to biases in the functional assessment of NGS data. BMC Genomics. 2024 May 8;25:455. Li Y, Shah-Simpson S, Okrah K, Belew AT, Choi J, Caradonna KL, et al. Transcriptome Remodeling in Trypanosoma cruzi and Human Cells during Intracellular Infection. PLOS Pathog. 2016 Apr 5;12(4):e1005511. Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008 July;5(7):621–8. Hicks SC, Townes FW, Teng M, Irizarry RA. Missing data and technical variability in single-cell RNA-sequencing experiments. Biostatistics. 2018 Oct 1;19(4):562–78. Ning L, Liu G, Li G, Hou Y, Tong Y, He J. Current Challenges in the Bioinformatics of Single Cell Genomics. Front Oncol [Internet]. 2014 Jan 27 [cited 2025 Aug 27];4. Available from: https://www.frontiersin.org/journals/oncology/articles/10.3389/fonc.2014.00007/full Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011 Dec;12(1):1–16. Patro R, Duggal G, Love MI, Irizarry RA, Kingsford C. Salmon provides fast and bias-aware quantification of transcript expression. Nat Methods. 2017 Apr;14(4):417–9. Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009 Mar 4;10(3):R25. Li H, Durbin R. Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics. 2009 July 15;25(14):1754–60. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015 Apr;12(4):357–60. Liao Y, Smyth GK, Shi W. The Subread aligner: fast, accurate and scalable read mapping by seed-and-vote. Nucleic Acids Res. 2013 May 1;41(10):e108. Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014 Apr 1;30(7):923–30. Bray NL, Pimentel H, Melsted P, Pachter L. Near-optimal probabilistic RNA-seq quantification. Nat Biotechnol. 2016 May;34(5):525–7. Berná L, Chiribao ML, Pita S, Alvarez-Valin F, Parodi-Talice A. Exploring the genomic landscape of the GP63 family in Trypanosoma cruzi: Evolutionary dynamics and functional peculiarities. PLoS Negl Trop Dis. 2025 Mar 17;19(3):e0012950. Berná L, Rodriguez M, Chiribao ML, Parodi-Talice A, Pita S, Rijo G, et al. Expanding an expanded genome: long-read sequencing of Trypanosoma cruzi. Microb Genomics. 2018;4(5):e000177. Dean AAC, Berná L, Robello C, Buscaglia CA, Balouz V. An algorithm for annotation and classification of T. cruzi MASP sequences: towards a better understanding of the parasite genetic variability. BMC Genomics. 2025 Dec;26(1):1–17. El-Sayed NM, Myler PJ, Blandin G, Berriman M, Crabtree J, Aggarwal G, et al. Comparative genomics of trypanosomatid parasitic protozoa. Science. 2005 July 15;309(5733):404–9. Talavera-López C, Messenger LA, Lewis MD, Yeo M, Reis-Cunha JL, Matos GM, et al. Repeat-Driven Generation of Antigenic Diversity in a Major Human Pathogen, Trypanosoma cruzi. Front Cell Infect Microbiol [Internet]. 2021 Mar 3 [cited 2025 Aug 20];11. Available from: https://www.frontiersin.org/journals/cellular-and-infection-microbiology/articles/10.3389/fcimb.2021.614665/full Albanaz ATS, Bartholomeu DC. Trypanosoma cruzi surface components: Why so many? Why so polymorphic? Curr Top Membr. 2025;95:175–214. Campetella O, Buscaglia CA, Mucci J, Leguizamón MS. Parasite-host glycan interactions during Trypanosoma cruzi infection: trans -Sialidase rides the show. Biochim Biophys Acta BBA - Mol Basis Dis. 2020 May 1;1866(5):165692. Choi YJ, Aliota MT, Mayhew GF, Erickson SM, Christensen BM. Dual RNA-seq of Parasite and Host Reveals Gene Expression Dynamics during Filarial Worm–Mosquito Interactions. PLoS Negl Trop Dis. 2014 May 22;8(5):e2905. Ehret T, Spork S, Dieterich C, Lucius R, Heitlinger E. Dual RNA-seq reveals no plastic transcriptional response of the coccidian parasite Eimeria falciformis to host immune defenses. BMC Genomics. 2017 Sept 5;18:686. Rastrojo A, Corvo L, Lombraña R, Solana JC, Aguado B, Requena JM. Analysis by RNA-seq of transcriptomic changes elicited by heat shock in Leishmania major. Sci Rep. 2019 May 6;9(1):6919. Shadab M, Das S, Banerjee A, Sinha R, Asad M, Kamran M, et al. RNA-Seq Revealed Expression of Many Novel Genes Associated With Leishmania donovani Persistence and Clearance in the Host Macrophage. Front Cell Infect Microbiol [Internet]. 2019 Feb 5 [cited 2025 Sept 30];9. Available from: https://www.frontiersin.org/journals/cellular-and-infection-microbiology/articles/10.3389/fcimb.2019.00017/full Watson AK, Williams TA, Williams BAP, Moore KA, Hirt RP, Embley TM. Transcriptomic profiling of host-parasite interactions in the microsporidian Trachipleistophora hominis. BMC Genomics. 2015 Dec;16(1):1–20. Radío S, Fort RS, Garat B, Sotelo-Silveira J, Smircich P. UTRme: A Scoring-Based Tool to Annotate Untranslated Regions in Trypanosomatid Genomes. Front Genet [Internet]. 2018 Dec 18 [cited 2025 Aug 20];9. Available from: https://www.frontiersin.org/journals/genetics/articles/10.3389/fgene.2018.00671/full Brandão A, Jiang T. The composition of untranslated regions in Trypanosoma cruzi genes. Parasitol Int. 2009 Sept 1;58(3):215–9. McDermaid A, Monier B, Zhao J, Liu B, Ma Q. Interpretation of differential gene expression results of RNA-seq data: review and integration. Brief Bioinform. 2019 Nov 27;20(6):2044–54. Inchausti L, Bilbao L, Campo VA, Garat J, Sotelo-Silveira J, Rinaldi G, et al. Single-cell RNA-seq reveals trans-sialidase-like superfamily gene expression heterogeneity in Trypanosoma cruzi populations. eLife [Internet]. 2025 Apr 2 [cited 2025 Sept 30];14. Available from: https://elifesciences.org/reviewed-preprints/105822 Laidlaw R, Sanchez MG, Pacheco JDS, Paradela LDS, Otto TD, Rycker MD. The Trypanosoma cruzi cell atlas; a single-cell resource for understanding parasite population heterogeneity and differentiation [Internet]. bioRxiv; 2025 [cited 2025 Sept 12]. p. 2024.10.01.616042. Available from: https://www.biorxiv.org/content/10.1101/2024.10.01.616042v2 Hamaguchi Y, Zeng C, Hamada M. Impact of human gene annotations on RNA-seq differential expression analysis. BMC Genomics. 2021 Oct 8;22(1):730. Lawson ND, Li R, Shin M, Grosse A, Yukselen O, Stone OA, et al. An improved zebrafish transcriptome annotation for sensitive and comprehensive detection of cell type-specific genes. Stainier DY, Busch-Nentwich E, Crump G, Busch-Nentwich E, editors. eLife. 2020 Aug 24;9:e55792. Soneson C, Love MI, Robinson MD. Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences [Internet]. F1000Research; 2016 [cited 2025 Sept 26]. Available from: https://f1000research.com/articles/4-1521 Wang D, Liu Y, Zhang Y, Chen Q, Han Y, Hou W, et al. A real-world multi-center RNA-seq benchmarking study using the Quartet and MAQC reference materials. Nat Commun. 2024 July 22;15(1):6167. Hashimoto T, de Hoon MJL, Grimmond SM, Daub CO, Hayashizaki Y, Faulkner GJ. Probabilistic resolution of multi-mapping reads in massively parallel sequencing data using MuMRescueLite. Bioinforma Oxf Engl. 2009 Oct 1;25(19):2613–4. Zytnicki M. mmquant: how to count multi-mapping reads? BMC Bioinformatics. 2017 Sept 15;18(1):411. Axtell MJ. ShortStack: comprehensive annotation and quantification of small RNA genes. RNA N Y N. 2013 June;19(6):740–51. Morrissey A, Shi J, James DQ, Mahony S. Accurate allocation of multimapped reads enables regulatory element analysis at repeats. Genome Res. 2024 July 23;34(6):937–51. Hita A, Brocart G, Fernandez A, Rehmsmeier M, Alemany A, Schvartzman S. MGcount: a total RNA-seq quantification tool to address multi-mapping and multi-overlapping alignments ambiguity in non-coding transcripts. BMC Bioinformatics. 2022 Dec;23(1):1–21. Pardo-Palacios FJ, Wang D, Reese F, Diekhans M, Carbonell-Sala S, Williams B, et al. Systematic assessment of long-read RNA-seq methods for transcript identification and quantification. Nat Methods. 2024 July;21(7):1349–63. Ji HJ, Pertea M. Enhancing transcriptome expression quantification through accurate assignment of long RNA sequencing reads with TranSigner. Genome Biol. 2025 Dec;26(1):1–26. Zheng GXY, Terry JM, Belgrader P, Ryvkin P, Bent ZW, Wilson R, et al. Massively parallel digital transcriptional profiling of single cells. Nat Commun. 2017 Jan 16;8(1):14049. Kaminow B, Yunusov D, Dobin A. STARsolo: accurate, fast and versatile mapping/quantification of single-cell and single-nucleus RNA-seq data [Internet]. bioRxiv; 2021 [cited 2025 Sept 30]. p. 2021.05.05.442755. Available from: https://www.biorxiv.org/content/10.1101/2021.05.05.442755v1 Srivastava A, Malik L, Smith T, Sudbery I, Patro R. Alevin efficiently estimates accurate gene abundances from dscRNA-seq data. Genome Biol. 2019 Dec;20(1):1–16. Melsted P, Booeshaghi AS, Liu L, Gao F, Lu L, Min KH (Joseph), et al. Modular, efficient and constant-memory single-cell RNA-seq preprocessing. Nat Biotechnol. 2021 July;39(7):813–8. Hoyos Sanchez MC, Ospina Zapata HS, Suarez BD, Ospina C, Barbosa HJ, Carranza Martinez JC, et al. A phased genome assembly of a Colombian Trypanosoma cruzi TcI strain and the evolution of gene families. Sci Rep. 2024 Jan 24;14(1):2054. Rhie A, McCarthy SA, Fedrigo O, Damas J, Formenti G, Koren S, et al. Towards complete and error-free genome assemblies of all vertebrate species. Nature. 2021 Apr;592(7856):737–46. Hickey JW, Becker WR, Nevins SA, Horning A, Perez AE, Zhu C, et al. Organization of the human intestine at single-cell resolution. Nature. 2023 July;619(7970):572–84. The Computational Pan-Genomics Consortium. Computational pan-genomics: status, promises and challenges. Brief Bioinform. 2018 Jan 1;19(1):118–35. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014 Dec;15(12):1–21. Andrews S. FastQC: A Quality Control Tool for High Throughput Sequence Data. 2010; Available from: [Online]. Available online at: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ Joshi NA, Fass JN. Sickle: a sliding-window, adaptive, quality-based trimming tool for FastQ files. 2011. Wang W, Peng D, Baptista RP, Li Y, Kissinger JC, Tarleton RL. Strain-specific genome evolution in Trypanosoma cruzi, the agent of Chagas disease. PLOS Pathog. 2021 Jan 28;17(1):e1009254. Frazee AC, Jaffe AE, Langmead B, Leek JT. Polyester: simulating RNA-seq datasets with differential transcript expression. Bioinformatics. 2015 Sept 1;31(17):2778–84. Additional Declarations No competing interests reported. Supplementary Files SupplementaryFile1Benchmarking.pdf Additional File 1 (PDF format) Supplementary Tables Contains Supplementary Tables 1–3, including: ● Summary of RNA-seq quantification strategies used in T. cruzi studies (Table S1) ● Mapping statistics for each quantification strategy (Table S2) ● Multi-mapping levels across T. cruzi multigene families (Table S3) SupplementaryFile2Benchmarking.pdf Additional File 2 (PDF format) Supplementary Figures Contains Supplementary Figures 1–10, supporting the main results presented in the manuscript, including visualizations of mapping distributions, quantification accuracy, and benchmarking comparisons. Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-7888056","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":538082392,"identity":"21018fb7-ad1c-4689-9fb6-aff97ce58c65","order_by":0,"name":"Aldana A Cepeda Dean","email":"","orcid":"","institution":"Universidad Nacional de San Martín","correspondingAuthor":false,"prefix":"","firstName":"Aldana","middleName":"A Cepeda","lastName":"Dean","suffix":""},{"id":538082393,"identity":"17854607-930d-4ec7-a0af-4a03409e2687","order_by":1,"name":"Virginia Balouz","email":"","orcid":"","institution":"Universidad Nacional de San Martín","correspondingAuthor":false,"prefix":"","firstName":"Virginia","middleName":"","lastName":"Balouz","suffix":""},{"id":538082394,"identity":"e5514a07-b8c0-415d-bd6d-d89890324eb6","order_by":2,"name":"Carlos A Buscaglia","email":"","orcid":"","institution":"Universidad Nacional de San Martín","correspondingAuthor":false,"prefix":"","firstName":"Carlos","middleName":"A","lastName":"Buscaglia","suffix":""},{"id":538082395,"identity":"544690e6-7788-454b-9a3c-4ffbd17c65d4","order_by":3,"name":"Natalia Rego","email":"","orcid":"","institution":"Institut Pasteur de Montevideo","correspondingAuthor":false,"prefix":"","firstName":"Natalia","middleName":"","lastName":"Rego","suffix":""},{"id":538082396,"identity":"6726f2e3-918b-409d-93e6-55369d5b7280","order_by":4,"name":"Luisa Berná","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA30lEQVRIiWNgGAWjYBACPgkwZYMkdICAFjaIljQgZiZNy2FStEg3H/tcUHM+j7///MHHvHsY5PhuJLBJ8+DTInMsefaMY7eLJW4kMxvzPGMwliSoRSLHmJmH7XZiww1moMoDDIkbCGvJ/8zM8+9c4vzzh9l/A7XUE6Elh5mZt+1A4oYDyWzMQC0JBkT4xZiZty+52PBGsrHknAMShjPPPGy2nINHC79082Nmnm92eXLnDz788OaAjTzf8eSDN97g0QIDCVAaFE2MDUz4HIauBQIYfxChZRSMglEwCkYMAACP+EcnVIOTFAAAAABJRU5ErkJggg==","orcid":"","institution":"Institut Pasteur de Montevideo","correspondingAuthor":true,"prefix":"","firstName":"Luisa","middleName":"","lastName":"Berná","suffix":""}],"badges":[],"createdAt":"2025-10-17 15:23:17","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-7888056/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-7888056/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":95291860,"identity":"9c9a7c09-60a3-461d-8d11-4366ed6217cb","added_by":"auto","created_at":"2025-11-06 11:08:51","extension":"docx","order_by":0,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":5335317,"visible":true,"origin":"","legend":"","description":"","filename":"ManuscritoBenchmarkingRNAseq.docx","url":"https://assets-eu.researchsquare.com/files/rs-7888056/v1/679bbf087804bf28baccef74.docx"},{"id":95291942,"identity":"ea4fefc4-0f26-49c6-b1ed-f747661d4805","added_by":"auto","created_at":"2025-11-06 11:08:54","extension":"json","order_by":1,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":7920,"visible":true,"origin":"","legend":"","description":"","filename":"ca430803c38f4217a1f009123c614724.json","url":"https://assets-eu.researchsquare.com/files/rs-7888056/v1/e302dfce4dd831db9801ad56.json"},{"id":95291943,"identity":"795861b4-ac48-4303-a1f2-eb3fbd84dafd","added_by":"auto","created_at":"2025-11-06 11:08:54","extension":"pdf","order_by":2,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":137685,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementaryFile1Benchmarking.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7888056/v1/f81de3078be8ddeaad03ef30.pdf"},{"id":95291852,"identity":"db5118c4-3652-4f28-a145-0131016bfda2","added_by":"auto","created_at":"2025-11-06 11:08:49","extension":"pdf","order_by":3,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":3200301,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementaryFile2Benchmarking.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7888056/v1/b02fc8548c58b75d0a759797.pdf"},{"id":95291868,"identity":"a8a3c567-88a5-4521-b9f8-78b58354a0d1","added_by":"auto","created_at":"2025-11-06 11:08:54","extension":"xml","order_by":4,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":140372,"visible":true,"origin":"","legend":"","description":"","filename":"ca430803c38f4217a1f009123c6147241enriched.xml","url":"https://assets-eu.researchsquare.com/files/rs-7888056/v1/43b9120fbb750432dca773d4.xml"},{"id":95291866,"identity":"f9ef4264-d80b-415a-b035-9978e9107656","added_by":"auto","created_at":"2025-11-06 11:08:53","extension":"png","order_by":5,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":544664,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage1.png","url":"https://assets-eu.researchsquare.com/files/rs-7888056/v1/2ab513424ca2d24626ce3d96.png"},{"id":95314724,"identity":"10023250-6a07-4042-985f-26a6b441129b","added_by":"auto","created_at":"2025-11-06 15:53:15","extension":"png","order_by":6,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":128092,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage2.png","url":"https://assets-eu.researchsquare.com/files/rs-7888056/v1/46157c6087cb61c84ba843ad.png"},{"id":95291951,"identity":"7e6419a1-1fe3-471a-9457-70b44d0cf15d","added_by":"auto","created_at":"2025-11-06 11:08:56","extension":"png","order_by":7,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":994539,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage3.png","url":"https://assets-eu.researchsquare.com/files/rs-7888056/v1/5dcee574fd828f49ea652432.png"},{"id":95291853,"identity":"c027e95b-ebf9-415a-818c-61f3df9cf1af","added_by":"auto","created_at":"2025-11-06 11:08:49","extension":"png","order_by":8,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":1506818,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage4.png","url":"https://assets-eu.researchsquare.com/files/rs-7888056/v1/d110f9181514c9b97a7025a0.png"},{"id":95291948,"identity":"fc0abb4c-bb48-4ccc-9461-2ced7c36dd4e","added_by":"auto","created_at":"2025-11-06 11:08:55","extension":"png","order_by":9,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":689994,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage5.png","url":"https://assets-eu.researchsquare.com/files/rs-7888056/v1/08debb6d7aa99256c7220a85.png"},{"id":95291946,"identity":"64663416-da3e-4fe2-8774-04e8631a21ed","added_by":"auto","created_at":"2025-11-06 11:08:55","extension":"png","order_by":10,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":1339909,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage6.png","url":"https://assets-eu.researchsquare.com/files/rs-7888056/v1/67a3838e09901bf52cf82083.png"},{"id":95291857,"identity":"df91dead-929b-47c7-83ba-ad5c629de3a3","added_by":"auto","created_at":"2025-11-06 11:08:50","extension":"png","order_by":11,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":72094,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage1.png","url":"https://assets-eu.researchsquare.com/files/rs-7888056/v1/45a90e4e264b8dcca7291452.png"},{"id":95291952,"identity":"8a1a12bd-6e30-4083-8a11-90abe05d0422","added_by":"auto","created_at":"2025-11-06 11:08:56","extension":"png","order_by":12,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":21984,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage2.png","url":"https://assets-eu.researchsquare.com/files/rs-7888056/v1/6ae8d6e05b871c11e8b273ab.png"},{"id":95291867,"identity":"106fc6bd-dff4-42fc-bba9-c643f1be24b6","added_by":"auto","created_at":"2025-11-06 11:08:53","extension":"png","order_by":13,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":202437,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage3.png","url":"https://assets-eu.researchsquare.com/files/rs-7888056/v1/4ffaf46e26932bc0dee1c948.png"},{"id":95291950,"identity":"51425802-0993-4d57-8739-cec6667f2c70","added_by":"auto","created_at":"2025-11-06 11:08:56","extension":"png","order_by":14,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":338156,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage4.png","url":"https://assets-eu.researchsquare.com/files/rs-7888056/v1/7ebad269e8c114a7fb5d34dc.png"},{"id":95313673,"identity":"f6de45b9-0a71-4d58-b648-48196f086591","added_by":"auto","created_at":"2025-11-06 15:51:50","extension":"png","order_by":15,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":139366,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage5.png","url":"https://assets-eu.researchsquare.com/files/rs-7888056/v1/23c982bd3981db5237073742.png"},{"id":95291849,"identity":"605e6eea-e7ef-48be-ac0a-7ba0a890f53e","added_by":"auto","created_at":"2025-11-06 11:08:49","extension":"png","order_by":16,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":230659,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage6.png","url":"https://assets-eu.researchsquare.com/files/rs-7888056/v1/6468c501bce35c86a55595d4.png"},{"id":95291944,"identity":"5d79497a-75cc-4ac0-9602-aac5ed008a9a","added_by":"auto","created_at":"2025-11-06 11:08:54","extension":"xml","order_by":17,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":138295,"visible":true,"origin":"","legend":"","description":"","filename":"ca430803c38f4217a1f009123c6147241structuring.xml","url":"https://assets-eu.researchsquare.com/files/rs-7888056/v1/60f4e101a2c9ee63c5c6a089.xml"},{"id":95291856,"identity":"b4e8c492-14f4-4c14-bdf9-d7e22f78b101","added_by":"auto","created_at":"2025-11-06 11:08:50","extension":"html","order_by":18,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":150670,"visible":true,"origin":"","legend":"","description":"","filename":"earlyproof.html","url":"https://assets-eu.researchsquare.com/files/rs-7888056/v1/797cdc46c4cace8a55c8dc7c.html"},{"id":95291851,"identity":"eef4a2c1-34e6-4657-9fa7-cc307b0d0d88","added_by":"auto","created_at":"2025-11-06 11:08:49","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":118984,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eDistribution of multi-mapping reads in \u003c/strong\u003e\u003cem\u003e\u003cstrong\u003eT. cruzi\u003c/strong\u003e\u003c/em\u003e\u003cstrong\u003e.\u003c/strong\u003e (a) Density distribution of the percentage of multi-mapping reads across all \u003cem\u003eT. cruzi\u003c/em\u003e genes, comparing different alignment and quasi-mapping strategies. Dashed vertical lines indicate the two most prominent density peaks for each strategy. Panels summarize the position of these peaks (in % multi-mapping) and their corresponding AUC, calculated by numerical integration of the KDE. (b) Pie chart showing the contribution of selected high-copy ‘multigene families’ (MASP, mucins, GP63, DGF-1 and TS) to \u003cem\u003eT. cruzi\u003c/em\u003e multi-mapping. (c) Density plots of multi-mapping percentages per strategy for each multigene family. (d) Violin plots depicting the distribution of multi-mapping percentages for multigene families with \u0026gt;0% multi-mapping (top) and \u0026gt;60% multi-mapping (bottom). Colors correspond to the strategies shown in (a).\u003c/p\u003e","description":"","filename":"1.png","url":"https://assets-eu.researchsquare.com/files/rs-7888056/v1/e6e404a3927b5cac1787526f.png"},{"id":95291858,"identity":"6418299b-1192-405f-b80c-e0af8f3f5a2f","added_by":"auto","created_at":"2025-11-06 11:08:51","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":134163,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eSchematic representation of the pipeline used to evaluate RNA-seq quantification in \u003c/strong\u003e\u003cem\u003e\u003cstrong\u003eT. cruzi\u003c/strong\u003e\u003c/em\u003e\u003cstrong\u003e. \u003c/strong\u003eExperimental RNA-seq data from \u003cem\u003eT. cruzi\u003c/em\u003e were processed using five combinations of aligners/mappers and quantifiers. Simulated RNA-seq datasets were generated based on expression profiles derived from real data and included additional scenarios with altered expression levels in selected genes or modified transcript annotations. Reads that map to a single or multiple transcripts are denoted in blue and grey, respectively. All datasets were processed through the same pipeline structure, and gene-level expression matrices were generated for comparative analysis.\u003c/p\u003e","description":"","filename":"2.png","url":"https://assets-eu.researchsquare.com/files/rs-7888056/v1/295278c6462f5f02ecfd7d62.png"},{"id":95291854,"identity":"a9071b69-59cf-4034-98f5-e808d6d27c5a","added_by":"auto","created_at":"2025-11-06 11:08:50","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":170871,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eComparative handling of multi-mapping reads across alignment strategies. \u003c/strong\u003eIn the scatter plots, each point represents a single gene quantified by the two strategies being compared. The color scheme corresponds to per-gene multi-mapping bins previously estimated. (a) Global view of read counts per gene, restricted to genes with read counts in the 0–5,000 range. (b) Zoomed view focusing on genes with 0–300 reads. The black dashed line represents the identity line (\u003cem\u003ey = x\u003c/em\u003e), while the grey area corresponds to the 95% confidence interval (CI) of the trend line, indicating the expected variability in the comparison. BOWTIE_FC, STAR_FC, and STAR_SALMON refer to the BOWTIE+featureCounts, STAR+featureCounts, and STAR+Salmon strategies, respectively.\u003c/p\u003e","description":"","filename":"3.png","url":"https://assets-eu.researchsquare.com/files/rs-7888056/v1/d495d24af20bcfb94708c6e1.png"},{"id":95291949,"identity":"e4a18d37-bdf7-49a5-b49a-92956f90a6ed","added_by":"auto","created_at":"2025-11-06 11:08:55","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":344683,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003ePerformance of quantification strategies in the simulated dataset. \u003c/strong\u003eComparison between simulated and observed gene-level counts across all \u003cem\u003eT. cruzi\u003c/em\u003e genes (n=17,650) using different alignment-based and alignment-free approaches. The left column displays the full range of simulated counts (0–5,000), while the right column zooms in on lowly expressed genes (0–300), corresponding to the gray-shaded region in the left panels. The black dashed line represents the identity line.\u003c/p\u003e","description":"","filename":"4.png","url":"https://assets-eu.researchsquare.com/files/rs-7888056/v1/8d04af9a853b54bbc26fa3a1.png"},{"id":95291945,"identity":"68ed866d-61a8-4dfd-a9ce-0eb684639961","added_by":"auto","created_at":"2025-11-06 11:08:55","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":194730,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eQuantification accuracy across highly paralogous genomic regions. \u003c/strong\u003eScatter plots of the difference between observed and simulated read counts across TS, GP63, mucins, and MASP genes sharing \u0026gt;95% sequence identity with at least one other family member of their corresponding family. Each point represents a gene, with the x-axis indicating the maximum sequence identity to its closest homolog within the same gene family. Shown are the median absolute errors (MedAE) of observed versus simulated expression values for each gene family and quantification method.\u003c/p\u003e","description":"","filename":"5.png","url":"https://assets-eu.researchsquare.com/files/rs-7888056/v1/a314787c7f355ba3f3ba66af.png"},{"id":95291928,"identity":"d5fc0d10-afde-4bb9-8504-c971d1334da5","added_by":"auto","created_at":"2025-11-06 11:08:54","extension":"png","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":266004,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eInclusion of UTR sequence in gene models improves detection and quantification of gene expression in \u003c/strong\u003e\u003cem\u003e\u003cstrong\u003eT. cruzi\u003c/strong\u003e\u003c/em\u003e\u003cstrong\u003e. \u003c/strong\u003eComparison between simulated and estimated counts for \u003cem\u003eT. cruzi\u003c/em\u003e genes using the pre-selected methodologies (STAR+Salmon, Salmon and Kallisto), on gene annotations including only CDS (left panels) or CDS+UTR (right panels). Black line represents the identity line, and the grey dashed line indicates the trend line. RMSE and MAE are shown.\u003c/p\u003e","description":"","filename":"6.png","url":"https://assets-eu.researchsquare.com/files/rs-7888056/v1/7bdbf11df2bcab08d49f08e7.png"},{"id":97895006,"identity":"4af8d092-390c-4ae5-bfef-895709721a3c","added_by":"auto","created_at":"2025-12-10 15:33:22","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":1976140,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7888056/v1/5e65fe86-640b-4c0a-bac0-49e4ac7015e6.pdf"},{"id":95291859,"identity":"173db6b7-cc6f-40e2-a103-c906502a931e","added_by":"auto","created_at":"2025-11-06 11:08:51","extension":"pdf","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":137685,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eAdditional File 1 \u003c/strong\u003e(PDF format)\u003cbr\u003e\n \u0026nbsp;Supplementary Tables\u003cstrong\u003e\u003cbr\u003e\n \u003c/strong\u003e\u0026nbsp;Contains Supplementary Tables 1–3, including:\u003c/p\u003e\n\u003cp\u003e● Summary of RNA-seq quantification strategies used in T. cruzi studies (Table S1)\u003c/p\u003e\n\u003cp\u003e● Mapping statistics for each quantification strategy (Table S2)\u003c/p\u003e\n\u003cp\u003e● Multi-mapping levels across T. cruzi multigene families (Table S3)\u003c/p\u003e","description":"","filename":"SupplementaryFile1Benchmarking.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7888056/v1/cb8fc391ecf4cb9b7cd8be96.pdf"},{"id":95291865,"identity":"59ab1cca-3fd1-4534-b5ef-18c0c5e3d064","added_by":"auto","created_at":"2025-11-06 11:08:53","extension":"pdf","order_by":2,"title":"","display":"","copyAsset":false,"role":"supplement","size":3200301,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eAdditional File 2 \u003c/strong\u003e(PDF format)\u003cbr\u003e\n \u0026nbsp;Supplementary Figures\u003cstrong\u003e\u003cbr\u003e\n \u003c/strong\u003e\u0026nbsp;Contains Supplementary Figures 1–10, supporting the main results presented in the manuscript, including visualizations of mapping distributions, quantification accuracy, and benchmarking comparisons.\u003c/p\u003e","description":"","filename":"SupplementaryFile2Benchmarking.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7888056/v1/5baf41897459952d3cdb283b.pdf"}],"financialInterests":"No competing interests reported.","formattedTitle":"The battle for reads: evaluating strategies to tackle multi-mapping in RNA-seq quantification in highly repetitive genomes","fulltext":[{"header":"BACKGROUND","content":"\u003cp\u003eRNA sequencing (RNA-seq) has become an essential tool in transcriptomics, enabling the accurate assessment of gene expression, discovery of novel transcripts, and analysis of transcript isoforms across diverse biological contexts. However, a persistent challenge in RNA-seq data analysis is the reliable quantification of transcripts derived from duplicated loci or repetitive genomic regions. In particular, the assignment of sequencing reads that map equally well to multiple loci\u0026ndash;commonly called multi-mappers\u0026ndash;poses a major computational hurdle (\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e, \u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eMulti-mappers can arise from various sources, including gene families with high sequence similarity among members, pseudogenes, repetitive sequences, and transposable elements. In organisms with highly repetitive genomes, this can account for a substantial fraction of reads\u0026ndash;up to 40%, depending on the sample and sequencing strategy (\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e). The severity of the problem depends on the read length, the nature and extent of genomic repetition, and the degree of sequence identity among paralogous regions. Discarding multi-mapped reads, a common step in many standard pipelines, can lead to biased expression estimates, particularly for multigene families and repetitive elements (\u003cspan additionalcitationids=\"CR4\" citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e). The issue is further exacerbated in single-cell RNA-seq studies, where short reads and sparse coverage amplify the uncertainty associated with ambiguous mappings (\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e, \u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e). A related challenge arises in complex eukaryotes, where genes produce multiple isoforms that share most of their sequence, and in certain pathogens that evolved large multigenic families as part of their strategy to achieve persistent infection in the host. This extensive sequence overlap makes it difficult to assign reads uniquely to specific transcripts, thereby complicating transcript-level quantification (\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e, \u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eTo address this issue, various computational strategies have been developed, encompassing both alignment and quantification steps. These methods differ in how they handle ambiguity and in the trade-offs they make between accuracy, interpretability, and computational efficiency. At the alignment step, tools like Bowtie and BWA\u0026ndash;commonly used for Illumina data\u0026ndash;perform fast, base-level mapping and are typically employed in organisms that lack introns, where spliced alignment is not required (\u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e, \u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e). In contrast, STAR is a splice-aware aligner that retains multiple mapping positions per read, providing greater flexibility for transcriptomic analyses, similar to other aligners such as HISAT2 and Subread (\u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e, \u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eAt the quantification step, traditional pipelines often rely on tools such as featureCounts to generate gene-level counts from aligned reads (\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e). However, these approaches may not fully account for the uncertainty introduced by multi-mapping reads. More recent methods, such as Salmon and Kallisto, employ alignment-free strategies based on k-mer indexing and probabilistic models (\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e, \u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e). Designed primarily for transcriptome quantification, these tools bypass the need for complete base-level alignment and incorporate statistical models to assign reads to transcripts while explicitly modeling mapping uncertainty. Their speed and ability to handle repetitive or paralogous sequences have made them increasingly popular, especially for large-scale or complex transcriptomes (\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e, \u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eDespite these advances, the challenge of accurately assigning multi-mapped reads remains. An underexplored issue concerns the alignment and gene quantification in eukaryotes whose genomes are not intron-rich but instead characterized by extensive repetition and structural complexity. In such contexts, where a large fraction of transcripts originates from highly similar loci, even the most advanced quantification tools may struggle to resolve ambiguous mappings. \u003cem\u003eTrypanosoma cruzi\u003c/em\u003e, the etiological agent of Chagas disease, exemplifies this problem. More than 50% of its genome consists of repetitive sequences, which comprise a wide variety of transposable elements and large multigene families with hundreds of copies, including functional genes and pseudogenes. Prominent among these are the \u003cem\u003etrans\u003c/em\u003e-sialidases (TS, ~\u0026thinsp;1400 genes), mucins (~\u0026thinsp;900 genes), mucin-associated surface proteins (MASP, ~\u0026thinsp;800 genes), dispersed gene family-1 (DGF-1, ~\u0026thinsp;80 genes) and GP63 proteases (~\u0026thinsp;200 genes) (\u003cspan additionalcitationids=\"CR17 CR18 CR19\" citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e). These gene families play critical roles in host-parasite interactions, immune evasion, and virulence; and are characterized by high sequence similarity, extensive though usually focalized allelic diversity, coordinated expression, and remarkable structural plasticity\u0026mdash;all of which compound the difficulty of resolving transcript-level expression from RNA-seq data (\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e, \u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eIn this study, we use \u003cem\u003eT. cruzi\u003c/em\u003e as a model to evaluate the performance of widely used RNA-seq analysis strategies under the most challenging conditions. We first analyze RNA-seq data, comparing outputs across pipelines with special focus on parasite multigenic families. We then simulate transcriptomic datasets to further assess the capacity of the best-performing tools to resolve ambiguous read assignments. Finally, we explore whether incorporating improved gene annotations, such as including untranslated regions (UTRs), can enhance the assignment of multi-mapped reads. Through this systematic approach, we aim to provide practical guidance for RNA-seq quantification in \u003cem\u003eT. cruzi\u003c/em\u003e and other organisms with challenging genomic architectures.\u003c/p\u003e"},{"header":"RESULTS","content":"\u003cp\u003e1. Distribution of Multi-Mapping Reads in \u003cem\u003eT. cruzi\u003c/em\u003e.\u003c/p\u003e\n\u003cp\u003eMulti-mapping reads\u0026mdash;those aligning equally well to multiple annotated features\u0026mdash;represent a major obstacle in RNA-seq quantification, particularly in genomes with high sequence redundancy. \u003cem\u003eTrypanosoma cruzi\u003c/em\u003e, characterized by an abundance of large multigene families and widespread repetitiveness, provides a compelling model to examine how different computational strategies handle this challenge.\u003c/p\u003e\n\u003cp\u003eTo quantify the extent of multi-mapping in \u003cem\u003eT. cruzi\u003c/em\u003e, we analyzed RNA-seq data using three representative tools: the aligners Bowtie2 and STAR, and the alignment-free quantifier Salmon. STAR and Salmon are among the most widely adopted in current RNA-seq workflows, while Bowtie2, though less common today, remains frequently used in studies of \u003cem\u003eT. cruzi\u003c/em\u003e and other intron-poor parasitic genomes (Additional file 1: Table.S1) (\u003cspan class=\"CitationRef\"\u003e23\u003c/span\u003e\u0026ndash;\u003cspan class=\"CitationRef\"\u003e27\u003c/span\u003e). The distribution of multi-mapping percentages across genes revealed a strikingly bimodal pattern for all three methods (Fig. \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003ea). One mode centered near 0% reflects genes covered almost exclusively by uniquely mapped reads, while a second mode\u0026mdash;ranging from ~\u0026thinsp;60% to ~\u0026thinsp;100%\u0026mdash;corresponds to genes with high levels of ambiguous alignments. However, the shape and intensity of these modes differed markedly across methods, reflecting differences in how each tool defines and processes multi-mapping. Salmon and STAR both showed a dominant peak at 0%, representing\u0026thinsp;~\u0026thinsp;42% and ~\u0026thinsp;36% of all genes, respectively, and a sharp secondary peak between 96\u0026ndash;97%, encompassing\u0026thinsp;~\u0026thinsp;20% of the genes. In contrast, Bowtie2 showed a flatter distribution with far fewer genes supported solely by uniquely mapped reads, and a broader peak centered around 82%, involving\u0026thinsp;~\u0026thinsp;54% of the genes. These results underscore the influence of mapping strategy on the detection and distribution of multi-mapping reads. A summary of total multi-mapping content per method is shown in Figure \u003cspan class=\"InternalRef\"\u003eS1\u003c/span\u003e, while detailed mapping statistics (including the percentage of aligned reads) are presented in Additional file 1: Table.S2.\u003c/p\u003e\n\u003cp\u003eFor comparative purposes, we performed similar analyses on evolutionarily related pathogens such as the trypanosomatids \u003cem\u003eTrypanosoma brucei\u003c/em\u003e and \u003cem\u003eLeishmania major\u003c/em\u003e, and the apicomplexans \u003cem\u003eToxoplasma gondii\u003c/em\u003e and \u003cem\u003ePlasmodium falciparum\u003c/em\u003e (Additional file 2: Fig.\u003cspan class=\"InternalRef\"\u003eS2\u003c/span\u003e). \u003cem\u003eT. brucei\u003c/em\u003e exhibited an overall level of multi-mapping comparable to that of \u003cem\u003eT. cruzi\u003c/em\u003e, but with a more sharply bimodal distribution: most genes showed either 0% or nearly 100% multi-mapping, with almost no genes falling in between. In contrast, \u003cem\u003eT. cruzi\u003c/em\u003e included a subset of genes with intermediate multi-mapping levels, resulting in a less polarized distribution. In the same line, \u003cem\u003eL. major\u003c/em\u003e showed a similar bimodal pattern, although with generally lower levels of ambiguity. In contrast, \u003cem\u003eToxoplasma\u003c/em\u003e and \u003cem\u003ePlasmodium\u003c/em\u003e presented minimal multi-mapping, with nearly all reads mapping unambiguously.\u003c/p\u003e\n\u003cp\u003eTo further understand how different strategies quantify ambiguous reads in \u003cem\u003eT. cruzi\u003c/em\u003e, we focused on five of its most expanded multigene families (MASP, Mucins, GP63, TS, and DGF-1). Each one of these families comprise between 100 and 1,000 copies showing varying extent of polymorphism (Additional file 1: Table.S3 and Additional file 2: Fig.S3), and collectively account for over 34% of multi-mapping reads in this organism (Fig. \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003eb). Bowtie2 yielded the highest overall density of multi-mapping, with a sharp peak above 60% across all families, reflecting its tendency to report multiple alignments without probabilistic resolution (Fig. \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003ec). In contrast, STAR and Salmon produced broader, more gradual distributions, suggesting a more nuanced handling of ambiguous reads (Fig. \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003ec and \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003ed).\u003c/p\u003e\n\u003cp\u003eInterestingly, distribution profiles varied between multigene families and the methods used. TS, for instance, consistently showed a bimodal pattern across strategies, whereas GP63 and mucins exhibited such bimodality only with STAR and Salmon (Fig. \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003ec and \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003ed). The GP63 family, in particular, showed more elevated levels of multi-mapping (Additional file 1: Table.S3), consistent with its internal structure (Additional file 2: Fig.S3), which includes subgroups of nearly identical genes that complicate read assignment (Bern\u0026aacute; et al., 2025). In contrast, DGF-1 exhibited markedly lower levels of multi-mapping, especially with STAR and Salmon, likely due to greater sequence divergence among its members (Fig. \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003ec and \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003ed).\u003c/p\u003e\n\u003cp\u003eTogether, these results reveal the high prevalence of multi-mapping reads in \u003cem\u003eT. cruzi\u003c/em\u003e, and their biased distribution across multi-copy gene families. They also demonstrate that the strategy used for the genomic anchoring of sequencing reads significantly impacts on the definition of the multi-mapping landscape in this parasite.\u003c/p\u003e\n\u003cp\u003e\u003cbr\u003e\u003c/p\u003e\n\u003cp\u003e2. Impact of Multi-Mapping on \u003cem\u003eT. cruzi\u003c/em\u003e Gene Expression Quantification\u003c/p\u003e\n\u003cp\u003eTo evaluate how mapping and quantification strategies affect transcript abundance estimates in \u003cem\u003eT. cruzi\u003c/em\u003e, we applied a standardized analysis pipeline to both experimental and simulated RNA-seq datasets (Fig. \u003cspan class=\"InternalRef\"\u003e2\u003c/span\u003e). We tested five widely used tool combinations\u0026mdash;Bowtie2\u0026thinsp;+\u0026thinsp;featureCounts, STAR\u0026thinsp;+\u0026thinsp;featureCounts, STAR\u0026thinsp;+\u0026thinsp;Salmon, Salmon, and Kallisto\u0026mdash;to assess their performance in quantifying gene expression, particularly for highly redundant gene families.\u003c/p\u003e\n\u003cp\u003eWe first compared gene-level read counts across the three alignment-based strategies using experimental RNA-seq data. Overall, ~\u0026thinsp;99% of \u003cem\u003eT. cruzi\u003c/em\u003e genes received fewer than 1,000 reads across pipelines (Fig. \u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003ea), and highly expressed genes (\u0026gt;\u0026thinsp;1,000 reads) were consistently detected. No correlation was observed between expression level and multi-mapping percentage, although increased variability\u0026mdash;consistent with stochastic noise\u0026mdash;was evident among lowly expressed genes (Additional file 2: Fig.S4).\u003c/p\u003e\n\u003cp\u003eWhen focusing on genes with up to 300 reads (covering\u0026thinsp;~\u0026thinsp;98% of the dataset), differences between pipelines became more apparent. STAR\u0026thinsp;+\u0026thinsp;Salmon and Bowtie2\u0026thinsp;+\u0026thinsp;featureCounts yielded similar distributions (Fig. \u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003eb), whereas STAR\u0026thinsp;+\u0026thinsp;featureCounts resulted in markedly lower counts for multi-mapping-prone genes. Interestingly, featureCounts assigned multi-mapping reads when combined with Bowtie2 but not when combined with STAR. This highlights how differences in the aligner output can affect read assignment and downstream quantification. As shown, featureCounts relies heavily on the flags and tags present in the BAM file, which are aligner-dependent, and although the performance of this tool can be modified by adjusting alignment parameters, this is rarely done in practice.\u003c/p\u003e\n\u003cp\u003e3. Performance Assessment Using Simulated Data\u003c/p\u003e\n\u003cp\u003eTo assess the ability of each quantification strategy to recover actual gene expression levels, we simulated RNA-seq datasets based on predefined expression profiles. Briefly, gene-level count values were specified manually from \u003cem\u003eT. cruzi\u003c/em\u003e experimental data and used as input for Polyester, which generated synthetic RNA-seq reads according to these target values. This approach allowed us to establish a known ground truth and systematically evaluate the performance of each pipeline in approximating simulated expression levels.\u003c/p\u003e\n\u003cp\u003eTo verify the stability of the simulation process, we compared the expression estimates across three replicates generated by Polyester. As no variability was observed between replicates for the three tested pipelines (Additional file 2: Fig.S5), we show a single replicate for all subsequent analyses.\u003c/p\u003e\n\u003cp\u003eGenome-wide analyses revealed that alignment-free transcript quantification approaches (Salmon and Kallisto) outperformed traditional alignment-based methods, producing expression estimates closer to the expected simulated values, with Kallisto showing the lowest overall error metrics (Fig. \u003cspan class=\"InternalRef\"\u003e4\u003c/span\u003e). Among the aligner-based strategies, both STAR\u0026thinsp;+\u0026thinsp;featureCounts and Bowtie2\u0026thinsp;+\u0026thinsp;featureCounts tended to overestimate transcript abundances, particularly for \u0026lsquo;highly-expressed\u0026rsquo; genes, whereas STAR\u0026thinsp;+\u0026thinsp;Salmon displayed a more consistent pattern, closer to that of the alignment-free tools, albeit with some dispersion (Fig. \u003cspan class=\"InternalRef\"\u003e4\u003c/span\u003e). These differences were further reflected when we assessed the number of genes with perfect quantification (i.e., identical observed and expected counts). As shown in Additional file 2: Fig.S6, Salmon achieved perfect counts for 8,417 genes (47% of the \u003cem\u003eT. cruzi\u003c/em\u003e genome), Kallisto for 7,180 (40%), while STAR\u0026thinsp;+\u0026thinsp;Salmon correctly quantified as low as 243 genes (1.4%).\u003c/p\u003e\n\u003cp\u003eQuite similar results were observed when we focused the analysis on parasite multigene families (Additional file 2: Fig.S7) hence reinforcing the robustness of these observations. Overall, these data indicate that RNA-seq analysis using alignment-free transcript quantification approaches (Salmon, Kallisto, and STAR\u0026thinsp;+\u0026thinsp;Salmon) provides for a more accurate gene expression quantification in \u003cem\u003eT. cruzi\u003c/em\u003e.\u003c/p\u003e\n\u003cp\u003e4. Assessing the Individual Expression of Multigene Family Members: The Battle Royale for the Reads.\u003c/p\u003e\n\u003cp\u003eDepending on their mode of evolution, multigene families in \u003cem\u003eT. cruzi\u003c/em\u003e differ in their internal sequence similarity, which may range from moderate to near-identical. This variation, illustrated in Additional file 2: Fig.S3, affects the resolution capacity of mapping tools, as high similarity can hinder the unambiguous assignment of reads. To examine how the best-performing strategies (STAR\u0026thinsp;+\u0026thinsp;Salmon, Salmon and Kallisto) handled this issue, we designed a simulation-based analysis focused on TS, GP63, mucins, and MASP genes. Specifically, we compared the expected and observed read counts for selected subsets of these gene families comprising genes sharing\u0026thinsp;\u0026gt;\u0026thinsp;95% sequence identity among them.\u003c/p\u003e\n\u003cp\u003eIn all cases, and consistent with our previous results, alignment-free methods such as Salmon and Kallisto performed better than the alignment-based STAR\u0026thinsp;+\u0026thinsp;Salmon approach (Fig.\u0026nbsp;5). This trend translated into much lower MedAE values for Salmon and Kallisto as compared to STAR\u0026thinsp;+\u0026thinsp;Salmon. As expected, MedAE values for each gene family and quantification method increased as sequence identity approached 100% (Fig.\u0026nbsp;5).\u003c/p\u003e\n\u003cp\u003eTo further explore the resolution capacity of each method, we designed a complementary simulated scenario in which overexpression was artificially introduced in one gene of selected paralogous subgroups. In this setup, a single gene within a group of highly similar paralogs was assigned with 1,000 reads, while the remaining members retained their original simulated expression levels. This allowed us to test the performance of the quantification methods under the most stringent conditions. (Additional file 2: Fig.S8). In the TS family (97.3\u0026ndash;97.75% identity), all strategies accurately recovered the simulated overexpression. For GP63 genes (99.23\u0026ndash;99.65% identity), however, Salmon and Kallisto captured the overexpression reliably, while STAR\u0026thinsp;+\u0026thinsp;Salmon tended to distribute part of the signal across other family members. In the Mucin family (98.9\u0026ndash;99.39% identity), only Salmon approximated the expected values, whereas the other methods largely failed to detect the overexpressed gene (Additional file 2: Fig.S8). The MASP family, which included a pair of genes with 100% identity, presented a limit case. Kallisto split the reads nearly equally between both copies, consistent with its expectation-maximization framework (EM-based). In contrast, Salmon and STAR\u0026thinsp;+\u0026thinsp;Salmon assigned all reads to a single gene, completely ignoring the identical paralog (Additional file 2: Fig.S8). These results underscore the fundamental challenge of quantifying genes or transcripts with identical or near-identical sequences, where some methods retain multi-mapping uncertainty while others introduce assignment bias by artificially resolving the ambiguity. More importantly, these findings highlight the superior accuracy and resolution of alignment-free methods in quantifying expression in multigene families with high sequence similarity.\u003c/p\u003e\n\u003cp\u003e5. Effect of Including UTRs in the Quantification\u003c/p\u003e\n\u003cp\u003eBuilding on the three selected strategies, we further explored the impact of transcript annotation quality on quantification accuracy. To that end, we compared simulated and observed read counts under two gene annotation schemes: one including only the CDS, and another one incorporating both the CDS and the UTRs (Fig. \u003cspan class=\"InternalRef\"\u003e2\u003c/span\u003e).\u003c/p\u003e\n\u003cp\u003eBecause \u003cem\u003eT. cruzi\u003c/em\u003e genome annotations typically lack UTRs, we predicted these regions using the UTRme tool (\u003cspan class=\"CitationRef\"\u003e28\u003c/span\u003e), based on RNA-Seq data from the YC6 strain. This exercise yielded a mean length of ~\u0026thinsp;200 bp and a median length of 82 bp for 5\u0026prime; UTRs and a mean length of ~\u0026thinsp;500 bp and a median of ~\u0026thinsp;300 bp for 3\u0026prime; UTRs (Additional file 2: Fig.S9), thus consistent with previous reports (\u003cspan class=\"CitationRef\"\u003e29\u003c/span\u003e). The resulting UTRs were integrated into the original annotation to generate an extended GFF file. As shown, the inclusion of UTRs consistently improved concordance between observed and expected expression levels across all methods (Fig. \u003cspan class=\"InternalRef\"\u003e6\u003c/span\u003e and Additional file 2: Fig.S10a-b). Both the RMSE and the MAE decreased when extended annotations were used, underscoring the importance of complete transcript models for accurate quantification. For STAR\u0026thinsp;+\u0026thinsp;Salmon (Fig. \u003cspan class=\"InternalRef\"\u003e6\u003c/span\u003e, upper panel), CDS-only annotations resulted in substantial dispersion around the identity line, with both over- and under-estimation of gene expression levels. In contrast, annotations including UTRs led to a tighter distribution of observed and expected values and thereby to lower error metrics. A similar trend was observed for Salmon alone (Fig. \u003cspan class=\"InternalRef\"\u003e6\u003c/span\u003e, middle panel). However, it should be noted that although Salmon more closely followed expected values than STAR\u0026thinsp;+\u0026thinsp;Salmon for many genes (especially those with moderate expression), it yielded a higher overall RMSE (515.39 vs. 406.21). This increase was driven by a small number of multi-mapping-rich transcripts, where Salmon assigned all reads to a single representative transcript, leaving others with zero or near-zero counts (Additional file 2: Fig.S10 c).\u003c/p\u003e\n\u003cp\u003eMost importantly, Kallisto (Fig. \u003cspan class=\"InternalRef\"\u003e6\u003c/span\u003e, bottom panel) yielded the best overall performance in terms of both RMSE and MAE. Its estimates closely followed the identity line with minimal dispersion, indicating high accuracy. The positive effect of UTR inclusion was especially marked, further underscoring the value of full-length transcript definitions.\u003c/p\u003e\n\u003cp\u003eTogether, these results confirm that extended annotations improve the accuracy and consistency of gene expression estimates, particularly when using transcript-level quantification tools.\u003c/p\u003e"},{"header":"DISCUSSION","content":"\u003cp\u003eAccurate RNA-seq quantification is largely hindered by multi-mapping reads, especially in genomes with high redundancy. In mammals, multi-mapping affects\u0026thinsp;~\u0026thinsp;5\u0026ndash;10% of reads (\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e), and can reach up to 30% depending on factors such as read length, sequencing layout, and the mapping algorithm used (\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e). In \u003cem\u003eT. cruzi\u003c/em\u003e, this issue is greatly amplified: depending on the tool, between 50\u0026ndash;80% of reads map ambiguously. A major contributor to this phenomenon is the outstanding expansion of multigene families coding for virulence factors that play critical roles in parasite infection, immune evasion, and pathogenesis (\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e, \u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e). Although multigene families represent\u0026thinsp;~\u0026thinsp;14% of \u003cem\u003eT. cruzi\u003c/em\u003e annotated genes, they account for over 34% of multi-mapped reads in our dataset\u0026mdash;consistent with their high copy number and sequence similarity. Moreover, the contribution of these gene families to the overall \u003cem\u003eT. cruzi\u003c/em\u003e multi-mapping is most likely higher than our gene-based estimations. As extensively reported, up to 80% of the total sequence dosage for these gene families is composed by transcriptionally active pseudogenes (\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e, \u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e). In addition to multi-copy gene families, the \u003cem\u003eT. cruzi\u003c/em\u003e genome is enriched in other repetitive sequences such as transposons and satellites that may further complicate accurate transcript expression quantification. Compared to other unicellular parasites, \u003cem\u003eT. cruzi\u003c/em\u003e displays a uniquely complex multi-mapping landscape, highlighting how genome architecture shapes mapping ambiguity and positioning this organism as a stringent benchmark for RNA-seq quantification tools.\u003c/p\u003e\u003cp\u003eOur comparative evaluation shows that alignment-free, transcriptome-aware quantifiers\u0026mdash;Salmon (\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e) and Kallisto (\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e)\u0026mdash; provide the most accurate expression estimates in \u003cem\u003eT. cruzi\u003c/em\u003e. These methods probabilistically assign reads, outperforming alignment-based workflows such as Bowtie2\u0026thinsp;+\u0026thinsp;featureCounts, which has been the most widely implemented strategy in \u003cem\u003eT. cruzi\u003c/em\u003e RNA-seq studies. Notably, they can resolve reads between paralogs sharing up to ~\u0026thinsp;98% sequence identity, although perfect or near-perfect duplicates remain unresolved. For use cases requiring alignment (e.g., visualization, variant calling), STAR combined with Salmon quantification offers a practical compromise with accuracy comparable to alignment-free methods.\u003c/p\u003e\u003cp\u003eAnnotation quality critically affects multi-mapping resolution. Incomplete gene models\u0026mdash;particularly those lacking UTRs\u0026mdash;inflate ambiguity by omitting sequence features that could distinguish otherwise similar transcripts. In \u003cem\u003eT. cruzi\u003c/em\u003e, where UTRs are not included in standard annotations, extending gene models substantially improved quantification accuracy (\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e, \u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e). Salmon and STAR\u0026thinsp;+\u0026thinsp;Salmon, in particular, better approximated the simulated values when extended annotations were used, underscoring their sensitivity to UTR-derived sequence diversity. These results are consistent with benchmarking studies showing that annotation and quantification are among the primary sources of variation in RNA-seq workflows (\u003cspan additionalcitationids=\"CR34 CR35\" citationid=\"CR33\" class=\"CitationRef\"\u003e33\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e36\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eKallisto also consistently achieved low MAE and RMSE across all conditions, though it benefited less from annotation refinement. This may reflect its robust transcript-level quantification strategy, which avoids zero estimates even when transcripts are identical. In contrast, Salmon tends to assign all reads to a single transcript when encountering exact duplicates, effectively silencing the others.\u003c/p\u003e\u003cp\u003eTogether, these findings underscore the central role of both quantification strategy and annotation quality in shaping RNA-seq outcomes. Beyond our benchmarking, a broader range of methods has been proposed to explicitly address multi-mapping.\u003c/p\u003e\u003cp\u003eSome tools aim to redistribute rather than discard ambiguous reads (\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e, \u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e37\u003c/span\u003e). MMquant (\u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e38\u003c/span\u003e), for example, assigns multi-mapped reads proportionally across all compatible features. This maximizes data usage and avoids underestimating repetitive gene families, but may overinflate expression in the presence of pseudogenes or paralogs. Other methods propose alternative allocation strategies: ShortStack (\u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e39\u003c/span\u003e), originally designed for small RNA-seq, assigns reads to the locus with the highest local density of uniquely mapped reads, an approach that can improve resolution in highly redundant regions, though it is limited to single-end data. Allo (\u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e40\u003c/span\u003e) a method developed for ChIP-seq, models the ambiguity explicitly using allocation likelihoods, and could potentially be adapted for RNA-seq contexts involving repetitive elements. MGcount offers a graph-based approach to quantify ambiguous reads across overlapping features and diverse RNA biotypes (\u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e41\u003c/span\u003e); while developed primarily for total-RNA and small RNA data, its underlying logic may hold promise for complex transcriptomes such as \u003cem\u003eT. cruzi\u003c/em\u003e.Despite their potential, these tools remain underutilized, due to their specific input requirements or implementation complexity in typical RNA-seq workflows.\u003c/p\u003e\u003cp\u003eLong-read RNA-seq technologies offer a promising alternative for resolving isoforms and paralogous transcripts that remain indistinguishable with short-read data. Despite technical challenges such as base-calling errors and incomplete transcript coverage, recent studies have shown that several dedicated tools perform well across different long-read platforms. Although short-read methods still outperform long-read-based approaches in terms of quantification accuracy under most conditions, combining both technologies may enhance transcript quantification in complex transcriptomes. As tools and protocols continue to improve, long-read RNA-seq is expected to provide increasingly accurate abundance estimates in the near future (\u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e42\u003c/span\u003e, \u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e43\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eSingle-cell RNA-seq adds further complexity. Many pipelines, such as Cell Ranger (\u003cspan citationid=\"CR44\" class=\"CitationRef\"\u003e44\u003c/span\u003e), discard multi-mapped reads by default to avoid false positives. This conservative approach can underrepresent highly repetitive gene families. In contrast, tools like STARsolo (\u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e45\u003c/span\u003e), Salmon-Alevin (\u003cspan citationid=\"CR46\" class=\"CitationRef\"\u003e46\u003c/span\u003e) and Kallisto|bustools (\u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e47\u003c/span\u003e) use probabilistic mapping strategies that better handle ambiguity and reduce memory usage. Benchmarking studies (\u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e36\u003c/span\u003e, \u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e45\u003c/span\u003e) have shown that the choice of quantification method significantly impacts marker gene detection and cell-type classification.\u003c/p\u003e\u003cp\u003eOur findings in \u003cem\u003eT. cruzi\u003c/em\u003e mirror these observations. Conservative pipelines minimize false positives but risk overlooking biologically relevant signals, while inclusive strategies may inflate expression in repetitive regions. This trade-off is evident in recent single-cell studies. The \u003cem\u003eT. cruzi\u003c/em\u003e cell atlas (\u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e) employed Cell Ranger, likely underestimating expression of surface protein families. In contrast, Inchausti et al. (2025) used Kallisto|bustools with EM-based redistribution and empirically defined 3\u0026prime; UTR annotations, enabling better assignment of reads and resolution of paralogous transcripts. These examples reinforce how analytical choices around multi-mapping and annotation directly shape the detection of gene family expression in \u003cem\u003eT. cruzi.\u003c/em\u003e\u003c/p\u003e\u003cp\u003eAnother emerging challenge is posed by the increasing availability of haplotype-resolved assemblies (\u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e48\u003c/span\u003e, \u003cspan citationid=\"CR49\" class=\"CitationRef\"\u003e49\u003c/span\u003e), as reads may map equally well to alternative haplotypes or to true paralogs, confounding quantification. Current linear-reference genome approaches cannot resolve these cases. Although haplotype-aware tools like STARconsensus offer partial solutions, systematic allele-aware quantification methods are still in development. Graph-based reference genomes and pangenomes (\u003cspan citationid=\"CR50\" class=\"CitationRef\"\u003e50\u003c/span\u003e, \u003cspan citationid=\"CR51\" class=\"CitationRef\"\u003e51\u003c/span\u003e) represent a promising shift. By modeling genetic variation directly, these approaches may improve read assignment and reduce reference bias. Their integration into RNA-seq workflows remains limited, but for organisms like \u003cem\u003eT. cruzi\u003c/em\u003e, which combine repeat-rich genomes and high haplotype divergence, they could prove transformative.\u003c/p\u003e\u003cp\u003eOur results argue against the practice of masking repetitive regions in transcriptomic analyses. With improved annotations and modern quantification methods, expression from these regions can be recovered with reasonable accuracy\u0026mdash;unlocking biologically meaningful information. At the same time, our study highlights remaining limitations: recent duplicates remain indistinguishable, annotations are still incomplete, and systematic biases persist depending on the aligner\u0026ndash;quantifier combination used. Advancing toward accurate transcript quantification in complex genomes will require continued development of haplotype-aware algorithms and integration of pangenomic frameworks.\u003c/p\u003e\u003cp\u003eIn this context, \u003cem\u003eT. cruzi\u003c/em\u003e emerges not only as a pathogen of biomedical relevance but also as a valuable system for testing the limits of current RNA-seq methods. Our study complements earlier benchmarking work by Love and Soneson (\u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e35\u003c/span\u003e, \u003cspan citationid=\"CR52\" class=\"CitationRef\"\u003e52\u003c/span\u003e) extending their findings to a genome characterized by extreme redundancy and structural plasticity. As transcriptomics moves into the pangenome era, organisms like \u003cem\u003eT. cruzi\u003c/em\u003e offer critical opportunities to stress-test the next generation of analytical tools.\u003c/p\u003e"},{"header":"CONCLUSIONS","content":"\u003cp\u003eOur study demonstrates that \u003cem\u003eT. cruzi\u003c/em\u003e provides a stringent benchmark for evaluating RNA-seq quantification in repetitive genomes. Multi-mapping reads, which affect only\u0026thinsp;~\u0026thinsp;5\u0026ndash;10% of polyadenylated transcripts in mammals, impact over half of \u003cem\u003eT. cruzi\u003c/em\u003e genes and are disproportionately concentrated in large multigene families. Among tested strategies, alignment-free methods (Salmon, Kallisto) most closely recapitulated simulated values, while STAR\u0026thinsp;+\u0026thinsp;Salmon also performed robustly when supported by extended annotations. Notably, the inclusion of UTRs substantially improved quantification across pipelines, underscoring annotation quality as a key determinant of accuracy. Beyond \u003cem\u003eT. cruzi\u003c/em\u003e, these findings highlight the broader need for probabilistic, annotation-rich, and ultimately haplotype-aware frameworks to resolve expression in complex transcriptomes. As graph-based and pangenome approaches mature, parasites with highly repetitive genomes such as \u003cem\u003eT. cruzi\u003c/em\u003e will serve as critical models for testing and refining next-generation RNA-seq quantification algorithms.\u003c/p\u003e"},{"header":"METHODS","content":"\u003cp\u003eThis study focuses on evaluating the impact of multi-mapping in RNA-seq analyses within the context of Illumina sequencing, which remains the predominant platform for quantifying gene expression. Accordingly, the described methodologies are either explicitly tailored to, or inherently compatible with, the characteristics of Illumina short-read data.\u003c/p\u003e\n\u003ch3\u003eDatasets and Quality Control\u003c/h3\u003e\n\u003cp\u003eRNA-Seq data from epimastigote and trypomastigote forms of \u003cem\u003eT. cruzi\u003c/em\u003e (YC6 strain) were obtained from Li et al. (2016) and are available in the NCBI SRA under projects PRJNA251582 and PRJNA251583. Quality assessment of each dataset was performed using FastQC (\u003cspan citationid=\"CR53\" class=\"CitationRef\"\u003e53\u003c/span\u003e), version 0.11.9. Raw fastq files were trimmed using Sickle (\u003cspan citationid=\"CR54\" class=\"CitationRef\"\u003e54\u003c/span\u003e) with parameters: -t Illumina -q 30 -l 70, removing low-quality bases (Phred score\u0026thinsp;\u0026lt;\u0026thinsp;30) and discarding reads shorter than 70 nucleotides. This ensured high-quality data for downstream analysis.\u003c/p\u003e\n\u003ch3\u003eQuantification of Multi-Mapping Reads Across Pathogen Genomes\u003c/h3\u003e\n\u003cp\u003eTo evaluate the extent of multi-mapping across different eukaryotic pathogens, we used publicly available paired-end RNA-seq datasets, which were aligned to the corresponding reference genomes and annotations. Genome assemblies, gene annotations (GFF), and RNA-seq data were obtained for \u003cem\u003eTrypanosoma brucei\u003c/em\u003e TREU927 (TriTrypDB release 68; ERR13925148), \u003cem\u003eLeishmania major\u003c/em\u003e Friedlin2021 (TriTrypDB release 68; ERR2604479), \u003cem\u003eToxoplasma gondii\u003c/em\u003e RH (ToxoDB release 68; SRR28421832) (Gajria et al., 2008), and \u003cem\u003ePlasmodium falciparum\u003c/em\u003e 3D7 (PlasmoDB GCA_000002765; SRR3455777). For \u003cem\u003eT. cruzi\u003c/em\u003e, we used the YC6 genome assembly and RNA-seq data as described previously. Reads were aligned using STAR (v2.7.11b) with up to 50 allowed mappings per read (--outFilterMultimapNmax 50). Genome indices were built with the respective genome assemblies and annotations. BAM files were then separated into uniquely and multi-mapped reads using SAMtools (v1.15.1), and gene-level quantification was performed with featureCounts (v2.0.1).\u003c/p\u003e\u003cdiv id=\"Sec8\" class=\"Section2\"\u003e\u003ch2\u003eRNA-seq Read Alignment and Quantification Pipelines\u003c/h2\u003e\u003cp\u003eTo evaluate the impact of mapping and quantification strategies on gene expression estimates in \u003cem\u003eT. cruzi\u003c/em\u003e, we implemented a standardized RNA-seq analysis pipeline integrating both experimental and simulated data. Paired-end RNA-seq libraries from \u003cem\u003eT. cruzi\u003c/em\u003e YC6 strain (\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e) were used as experimental input. Simulated datasets were generated based on these profiles and included scenarios with altered expression levels and annotation refinements. Overall, we tested five widely used tool combinations:\u003c/p\u003e\u003cp\u003e\u003cul\u003e\u003cli\u003e\u003cp\u003eBowtie2 (v.2.5.3)\u0026thinsp;+\u0026thinsp;featureCounts (v.2.0.3)\u003c/p\u003e\u003c/li\u003e\u003cli\u003e\u003cp\u003eSTAR (v.2.7.11)\u0026thinsp;+\u0026thinsp;featureCounts\u003c/p\u003e\u003c/li\u003e\u003cli\u003e\u003cp\u003eSTAR\u0026thinsp;+\u0026thinsp;Salmon (v.1.10.1, alignment-based mode)\u003c/p\u003e\u003c/li\u003e\u003cli\u003e\u003cp\u003eSalmon (v.1.10.1, quasi-mapping mode)\u003c/p\u003e\u003c/li\u003e\u003cli\u003e\u003cp\u003eKallisto (v.0.50.1, pseudo-alignment mode)\u003c/p\u003e\u003c/li\u003e\u003c/ul\u003e\u003c/p\u003e\u003cp\u003eIn all cases, reads were mapped to the \u003cem\u003eT. cruzi\u003c/em\u003e YC6 reference genome and transcriptome (\u003cspan citationid=\"CR55\" class=\"CitationRef\"\u003e55\u003c/span\u003e). For the Bowtie2 pipeline, reads were aligned in end-to-end mode, using the --sensitive preset, which balances speed and alignment accuracy in the presence of mismatches or indels. Additional options included --no-unal to exclude unaligned reads from the output. Aligned reads were quantified at the gene level using featureCounts. Transcript-level annotation was enabled via the -g transcript_id option and a modified GTF file including the UTRs (see below). Multi-mapping reads were not discarded (-M) and included in the quantification using --fraction, a minimum fractional overlap to the annotation feature of 10% was required (--fracOverlap 0.1). For the STAR pipeline, reads were aligned, allowing up to 50 multiple mappings per read. To tailor mammalian-trained STAR splice junction detection functions to \u003cem\u003eT. cruzi\u003c/em\u003e genome, the following options were used: --outSJfilterReads Unique, --outSJfilterOverhangMin 150 150 150 150, and --outFilterType BySJout. Mismatch filtering was set with --outFilterMismatchNoverReadLmax 0.01, and reads were aligned in end-to-end mode (--alignEndsType EndToEnd). The output included sorted BAM files (--outSAMtype BAM SortedByCoordinate) and transcriptome alignments (--quantMode TranscriptomeSAM). These alignments were either used directly with featureCounts or passed to Salmon for quantification in alignment-based mode. For Salmon in quasi-mapping mode, transcript quantification was performed using the genome and GTF annotation (both indexed) with default parameters. Paired-end reads were processed with automatic library type detection (-l A), and the --validateMappings option was enabled to improve mapping accuracy. Finally, Kallisto was used to quantify transcript abundance via pseudo-alignment with default parameters, using the same input fasta file with the \u003cem\u003eT. cruzi\u003c/em\u003e transcripts. To calculate the number of mapped reads for each methodology, we employed samtools flagstat (v.1.20), which provides global alignment statistics from the corresponding BAM files.\u003c/p\u003e\u003cp\u003eGene-level count matrices were generated for all strategies. Genes were classified into major multigene families (MASP, mucins, TS, GP63) or \"Other\" for comparative analysis, particularly focusing on the impact of multi-mapping reads on quantification accuracy. Genes belonging to multigene families were retrieved from CDS and GFF annotations. Pairwise sequence identities were calculated using the Biopython pairwise2 module. Global alignments were performed with the globalxx scoring scheme, and sequence identity was defined as the proportion of identical positions relative to the maximum sequence length.\u003c/p\u003e\u003cp\u003e\u003cb\u003eUTR Refinement of\u003c/b\u003e \u003cb\u003eT. cruzi\u003c/b\u003e \u003cb\u003eGene Annotation\u003c/b\u003e\u003c/p\u003e\u003cp\u003eThe gene annotation file for \u003cem\u003eT. cruzi\u003c/em\u003e YC6 (TriTrypDB release 68) was refined to include UTRs, which are not present in the standard annotation. UTRs were predicted using UTRme (\u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e) with paired-end mode and default parameters, based on YC6 RNA-seq data. The highest-scoring UTRs were assigned to their corresponding genes using custom Python scripts and merged into the original GFF file.\u003c/p\u003e\u003c/div\u003e\n\u003ch3\u003eMulti-Mapping Analysis\u003c/h3\u003e\n\u003cp\u003eWe estimated the proportion of multi-mapping reads for each mapping strategy. Only STAR, Bowtie2, and Salmon were included in this analysis, as all generate BAM or SAM files that enable tracking of read alignments. Multi-mapping reads were identified based on alignment flags and mapping quality scores. For STAR, secondary alignments were used to identify multi-mapped reads. For Bowtie2, the presence of the XS:i: tag was used to infer alternative mapping locations. Unique and multi-mapped reads were counted separately using a feature-based counting approach. Final multi-mapping percentages were calculated as the ratio of multi-mapped reads to total mapped reads, both globally and at the gene level. For Salmon, transcript-level alignments were obtained using the --writeMappings option to produce SAM files. These files were filtered to remove unmapped reads, and transcript-level counts were extracted by summarizing mappings per transcript using custom awk commands in bash. Kernel density estimation (KDE) was applied to the distribution of multi-mapping percentages for each strategy. Peaks in the KDE curves were identified as local maxima, and the area under the curve (AUC) corresponding to each peak was calculated using trapezoidal numerical integration to assess their relative contribution. All downstream processing, statistical analyses, and visualizations were performed in Python using the pandas, scipy, and numpy libraries for data manipulation and numerical operations, and seaborn and matplotlib for plotting and figure customization.\u003c/p\u003e\n\u003ch3\u003eBenchmark of Quantification Strategies\u003c/h3\u003e\n\u003cp\u003eTo compare the impact of alignment and quantification methods on gene expression profiles, we evaluated the number of detected genes, transcript per million (TPM) distributions, expression correlations, and gene family representation across pipelines. Additionally, a public GitHub repository (Benchmarking-RNAseq, \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://github.com/AldanaCepedaDean/Benchmarking-RNAseq\u003c/span\u003e\u003cspan address=\"https://github.com/AldanaCepedaDean/Benchmarking-RNAseq\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003cspan type=\"Underline\" class=\"Underline\" name=\"Emphasis\"\u003e)\u003c/span\u003e was created to provide an overview of the methodologies discussed, together with example scripts and workflows to facilitate reproducibility and practical implementation.\u003c/p\u003e\u003cdiv id=\"Sec11\" class=\"Section2\"\u003e\u003ch2\u003eSimulation of RNA-seq Data\u003c/h2\u003e\u003cp\u003eSimulated RNA-seq reads were derived from real data for \u003cem\u003eT. cruzi\u003c/em\u003e (strain YC6, trypomastigote stage; SRR1346054). Baseline expression values were obtained by averaging gene-level counts across three quantification strategies: Bowtie2\u0026thinsp;+\u0026thinsp;featureCounts, STAR\u0026thinsp;+\u0026thinsp;featureCounts, and STAR\u0026thinsp;+\u0026thinsp;Salmon. These average counts were used to generate the reference expression table for simulations. Three simulation scenarios were implemented:\u003c/p\u003e\u003cp\u003e\u003col\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003e\u003cb\u003eReference simulation\u003c/b\u003e: Reads were simulated according to the averaged expression values described above.\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003e\u003cb\u003eMultigene family simulation\u003c/b\u003e: One representative gene was manually assigned with 1,000 reads, while the remaining family members retained their original expression values from the reference count table. This level of upregulation (from baseline expression levels ranging between 200 and 500 reads) was chosen to simulate a moderate and biologically plausible increase in transcript abundance.\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003e\u003cb\u003eUTR-expanded simulation\u003c/b\u003e: Same as the reference simulation, but using transcript annotations that include both CDS and UTRs during indexing and alignment.\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003c/ol\u003e\u003c/p\u003e\u003cp\u003eAll simulations were carried out using Polyester (\u003cspan citationid=\"CR56\" class=\"CitationRef\"\u003e56\u003c/span\u003e), generating 100-nucleotide paired-end reads with a 0.5% sequencing error rate, based on the error rate obtained from the YC6 reference transcriptome (\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e). Root Mean Square Error (RMSE) and Mean Absolute Error (MAE) were calculated to compare between simulated and observed values. For the simulation of multigene families, to assess the accuracy of expression estimates, we calculated the Median Absolute Error (MedAE) defined as:\u003cdiv id=\"Equa\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equa\" name=\"EquationSource\"\u003e\n$$\\:MedAE\\:=\\:median({|Observed}_{1}-{Expected}_{1}|,\\:...,\\:{|Observed}_{n}-{Expected}_{n}\\left|\\right)$$\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e\u003cp\u003ewhere \u003cem\u003eObserved\u003c/em\u003e\u003csub\u003e\u003cem\u003en\u003c/em\u003e\u003c/sub\u003e denotes the expression estimate from each quantification method for gene \u003cem\u003en\u003c/em\u003e, and \u003cem\u003eExpected\u003c/em\u003e\u003csub\u003e\u003cem\u003en\u003c/em\u003e​\u003c/sub\u003e denotes the corresponding ground-truth value from the simulation. All computations were performed in Python using the NumPy and scikit-learn libraries for metric calculations. Data handling, statistical analysis, and visualization were conducted with pandas, numpy, seaborn, and matplotlib.\u003c/p\u003e\u003c/div\u003e"},{"header":"Abbreviations","content":"\u003cp\u003eRNA-seq, RNA sequencing; UTRs, Untranslated Regions; TS, \u003cem\u003etrans\u003c/em\u003e-sialidase from \u003cem\u003eT. cruzi\u003c/em\u003e; MASP, mucin-associated surface protein from \u003cem\u003eT. cruzi\u003c/em\u003e; DGF-1, Dispersed Gene Family-1 from \u003cem\u003eT. cruzi\u003c/em\u003e\u003c/p\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eACKNOWLEDGMENTS\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eWe are grateful to Lucas Inchausti for his help with running simulated data and for his kindness, and to Carlos Robello for his involvement in coordinating and facilitating Aldana A. Cepeda-Dean\u0026rsquo;s stay during her UNU-BIOLAC scholarship.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eFUNDING\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eAACD holds a CONICET fellowship, whereas CAB and VB are career investigators from the same institution. AACD also received a CIN (Consejo Interuniversitario Nacional, Argentina) fellowship and a UNU-Biolac fellowship. LB and NR are PEDECIBA (Programa de Desarrollo de Ciencias B\u0026aacute;sicas, Uruguay) researchers and members of the Sistema Nacional de Investigadores (SNI, ANII, Uruguay). This investigation received financial support from the ANPCyT (PICT-2021-0284 to CAB) and by the Institut Pasteur de Montevideo and FOCEM - Fondo para la Convergencia Estructural del Mercosur (COF 03/11). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAUTHOR\u0026rsquo;S CONTRIBUTION\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eConceptualization: AACD, LB, NR, VB, CAB.\u003c/p\u003e\n\u003cp\u003eDesign: AACD, LB, NR.\u003c/p\u003e\n\u003cp\u003eFormal analysis: AACD, LB, NR.\u003c/p\u003e\n\u003cp\u003eMethodology: AACD, LB, NR.\u003c/p\u003e\n\u003cp\u003eFunding acquisition and resources: LB, NR, CAB, VB, AACD.\u003c/p\u003e\n\u003cp\u003eWriting - original draft: AACD, LB, NR.\u003c/p\u003e\n\u003cp\u003eWriting - review and editing: AACD, VB, CAB, LB, NR.\u003c/p\u003e\n\u003cp\u003eFinal approval of the version to be submitted: All the co-authors. All authors read and approved the final manuscript.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eETHICS APPROVAL AND CONSENT TO PARTICIPATE\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eNot applicable.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCONSENT TO PUBLICATION\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eNot applicable.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCOMPETING INTEREST\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe authors declare no competing interests.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAVAILABILITY OF DATA AND MATERIALS.\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eAll data generated and analysed during this study are included in this published article and its supplementary information files.\u003c/p\u003e\n\u003cp\u003eAdditionally, a public GitHub repository has been created to provide an overview of the RNA-seq methodologies discussed in this review, along with example scripts and workflows to facilitate reproducibility and practical use by the community:\u003c/p\u003e\n\u003cul\u003e\n \u003cli\u003e\u003cstrong\u003eProject name:\u003c/strong\u003e Benchmarking-RNAseq\u003c/li\u003e\n \u003cli\u003e\u003cstrong\u003eProject home page:\u003c/strong\u003ehttps://github.com/AldanaCepedaDean/Benchmarking-RNAseq\u003c/li\u003e\n \u003cli\u003e\u003cstrong\u003eOperating system(s):\u003c/strong\u003e Linux, MacOS, Windows (with WSL)\u003c/li\u003e\n \u003cli\u003e\u003cstrong\u003eProgramming language:\u003c/strong\u003e Shell, Bash, R\u003c/li\u003e\n \u003cli\u003e\u003cstrong\u003eOther requirements:\u003c/strong\u003e None\u003c/li\u003e\n\u003c/ul\u003e\n\u003cp\u003eThis repository is intended as a resource for researchers to explore and implement the discussed strategies for handling multi-mapping reads in RNA-seq experiments.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eDeschamps-Francoeur G, Simoneau J, Scott MS. Handling multi-mapped reads in RNA-seq. Comput Struct Biotechnol J. 2020 Jan 1;18:1569\u0026ndash;76.\u003c/li\u003e\n\u003cli\u003eSrivastava A, Malik L, Sarkar H, Zakeri M, Almodaresi F, Soneson C, et al. Alignment and mapping methodology influence transcript abundance estimation. Genome Biol. 2020 Sept 7;21(1):239.\u003c/li\u003e\n\u003cli\u003eAlmeida da Paz M, Warger S, Taher L. Disregarding multimappers leads to biases in the functional assessment of NGS data. BMC Genomics. 2024 May 8;25:455.\u003c/li\u003e\n\u003cli\u003eLi Y, Shah-Simpson S, Okrah K, Belew AT, Choi J, Caradonna KL, et al. Transcriptome Remodeling in Trypanosoma cruzi and Human Cells during Intracellular Infection. PLOS Pathog. 2016 Apr 5;12(4):e1005511.\u003c/li\u003e\n\u003cli\u003eMortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008 July;5(7):621\u0026ndash;8.\u003c/li\u003e\n\u003cli\u003eHicks SC, Townes FW, Teng M, Irizarry RA. Missing data and technical variability in single-cell RNA-sequencing experiments. Biostatistics. 2018 Oct 1;19(4):562\u0026ndash;78.\u003c/li\u003e\n\u003cli\u003eNing L, Liu G, Li G, Hou Y, Tong Y, He J. Current Challenges in the Bioinformatics of Single Cell Genomics. Front Oncol [Internet]. 2014 Jan 27 [cited 2025 Aug 27];4. Available from: https://www.frontiersin.org/journals/oncology/articles/10.3389/fonc.2014.00007/full\u003c/li\u003e\n\u003cli\u003eLi B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011 Dec;12(1):1\u0026ndash;16.\u003c/li\u003e\n\u003cli\u003ePatro R, Duggal G, Love MI, Irizarry RA, Kingsford C. Salmon provides fast and bias-aware quantification of transcript expression. Nat Methods. 2017 Apr;14(4):417\u0026ndash;9.\u003c/li\u003e\n\u003cli\u003eLangmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009 Mar 4;10(3):R25.\u003c/li\u003e\n\u003cli\u003eLi H, Durbin R. Fast and accurate short read alignment with Burrows\u0026ndash;Wheeler transform. Bioinformatics. 2009 July 15;25(14):1754\u0026ndash;60.\u003c/li\u003e\n\u003cli\u003eKim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015 Apr;12(4):357\u0026ndash;60.\u003c/li\u003e\n\u003cli\u003eLiao Y, Smyth GK, Shi W. The Subread aligner: fast, accurate and scalable read mapping by seed-and-vote. Nucleic Acids Res. 2013 May 1;41(10):e108.\u003c/li\u003e\n\u003cli\u003eLiao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014 Apr 1;30(7):923\u0026ndash;30.\u003c/li\u003e\n\u003cli\u003eBray NL, Pimentel H, Melsted P, Pachter L. Near-optimal probabilistic RNA-seq quantification. Nat Biotechnol. 2016 May;34(5):525\u0026ndash;7.\u003c/li\u003e\n\u003cli\u003eBern\u0026aacute; L, Chiribao ML, Pita S, Alvarez-Valin F, Parodi-Talice A. Exploring the genomic landscape of the GP63 family in Trypanosoma cruzi: Evolutionary dynamics and functional peculiarities. PLoS Negl Trop Dis. 2025 Mar 17;19(3):e0012950.\u003c/li\u003e\n\u003cli\u003eBern\u0026aacute; L, Rodriguez M, Chiribao ML, Parodi-Talice A, Pita S, Rijo G, et al. Expanding an expanded genome: long-read sequencing of Trypanosoma cruzi. Microb Genomics. 2018;4(5):e000177.\u003c/li\u003e\n\u003cli\u003eDean AAC, Bern\u0026aacute; L, Robello C, Buscaglia CA, Balouz V. An algorithm for annotation and classification of T. cruzi MASP sequences: towards a better understanding of the parasite genetic variability. BMC Genomics. 2025 Dec;26(1):1\u0026ndash;17.\u003c/li\u003e\n\u003cli\u003eEl-Sayed NM, Myler PJ, Blandin G, Berriman M, Crabtree J, Aggarwal G, et al. Comparative genomics of trypanosomatid parasitic protozoa. Science. 2005 July 15;309(5733):404\u0026ndash;9.\u003c/li\u003e\n\u003cli\u003eTalavera-L\u0026oacute;pez C, Messenger LA, Lewis MD, Yeo M, Reis-Cunha JL, Matos GM, et al. Repeat-Driven Generation of Antigenic Diversity in a Major Human Pathogen, Trypanosoma cruzi. Front Cell Infect Microbiol [Internet]. 2021 Mar 3 [cited 2025 Aug 20];11. Available from: https://www.frontiersin.org/journals/cellular-and-infection-microbiology/articles/10.3389/fcimb.2021.614665/full\u003c/li\u003e\n\u003cli\u003eAlbanaz ATS, Bartholomeu DC. Trypanosoma cruzi surface components: Why so many? Why so polymorphic? Curr Top Membr. 2025;95:175\u0026ndash;214.\u003c/li\u003e\n\u003cli\u003eCampetella O, Buscaglia CA, Mucci J, Leguizam\u0026oacute;n MS. Parasite-host glycan interactions during \u003cem\u003eTrypanosoma cruzi\u003c/em\u003e infection: \u003cem\u003etrans\u003c/em\u003e-Sialidase rides the show. Biochim Biophys Acta BBA - Mol Basis Dis. 2020 May 1;1866(5):165692.\u003c/li\u003e\n\u003cli\u003eChoi YJ, Aliota MT, Mayhew GF, Erickson SM, Christensen BM. Dual RNA-seq of Parasite and Host Reveals Gene Expression Dynamics during Filarial Worm\u0026ndash;Mosquito Interactions. PLoS Negl Trop Dis. 2014 May 22;8(5):e2905.\u003c/li\u003e\n\u003cli\u003eEhret T, Spork S, Dieterich C, Lucius R, Heitlinger E. Dual RNA-seq reveals no plastic transcriptional response of the coccidian parasite Eimeria falciformis to host immune defenses. BMC Genomics. 2017 Sept 5;18:686.\u003c/li\u003e\n\u003cli\u003eRastrojo A, Corvo L, Lombra\u0026ntilde;a R, Solana JC, Aguado B, Requena JM. Analysis by RNA-seq of transcriptomic changes elicited by heat shock in Leishmania major. Sci Rep. 2019 May 6;9(1):6919.\u003c/li\u003e\n\u003cli\u003eShadab M, Das S, Banerjee A, Sinha R, Asad M, Kamran M, et al. RNA-Seq Revealed Expression of Many Novel Genes Associated With Leishmania donovani Persistence and Clearance in the Host Macrophage. Front Cell Infect Microbiol [Internet]. 2019 Feb 5 [cited 2025 Sept 30];9. Available from: https://www.frontiersin.org/journals/cellular-and-infection-microbiology/articles/10.3389/fcimb.2019.00017/full\u003c/li\u003e\n\u003cli\u003eWatson AK, Williams TA, Williams BAP, Moore KA, Hirt RP, Embley TM. Transcriptomic profiling of host-parasite interactions in the microsporidian Trachipleistophora hominis. BMC Genomics. 2015 Dec;16(1):1\u0026ndash;20.\u003c/li\u003e\n\u003cli\u003eRad\u0026iacute;o S, Fort RS, Garat B, Sotelo-Silveira J, Smircich P. UTRme: A Scoring-Based Tool to Annotate Untranslated Regions in Trypanosomatid Genomes. Front Genet [Internet]. 2018 Dec 18 [cited 2025 Aug 20];9. Available from: https://www.frontiersin.org/journals/genetics/articles/10.3389/fgene.2018.00671/full\u003c/li\u003e\n\u003cli\u003eBrand\u0026atilde;o A, Jiang T. The composition of untranslated regions in \u003cem\u003eTrypanosoma cruzi\u003c/em\u003e genes. Parasitol Int. 2009 Sept 1;58(3):215\u0026ndash;9.\u003c/li\u003e\n\u003cli\u003eMcDermaid A, Monier B, Zhao J, Liu B, Ma Q. Interpretation of differential gene expression results of RNA-seq data: review and integration. Brief Bioinform. 2019 Nov 27;20(6):2044\u0026ndash;54.\u003c/li\u003e\n\u003cli\u003eInchausti L, Bilbao L, Campo VA, Garat J, Sotelo-Silveira J, Rinaldi G, et al. Single-cell RNA-seq reveals trans-sialidase-like superfamily gene expression heterogeneity in Trypanosoma cruzi populations. eLife [Internet]. 2025 Apr 2 [cited 2025 Sept 30];14. Available from: https://elifesciences.org/reviewed-preprints/105822\u003c/li\u003e\n\u003cli\u003eLaidlaw R, Sanchez MG, Pacheco JDS, Paradela LDS, Otto TD, Rycker MD. The Trypanosoma cruzi cell atlas; a single-cell resource for understanding parasite population heterogeneity and differentiation [Internet]. bioRxiv; 2025 [cited 2025 Sept 12]. p. 2024.10.01.616042. Available from: https://www.biorxiv.org/content/10.1101/2024.10.01.616042v2\u003c/li\u003e\n\u003cli\u003eHamaguchi Y, Zeng C, Hamada M. Impact of human gene annotations on RNA-seq differential expression analysis. BMC Genomics. 2021 Oct 8;22(1):730.\u003c/li\u003e\n\u003cli\u003eLawson ND, Li R, Shin M, Grosse A, Yukselen O, Stone OA, et al. An improved zebrafish transcriptome annotation for sensitive and comprehensive detection of cell type-specific genes. Stainier DY, Busch-Nentwich E, Crump G, Busch-Nentwich E, editors. eLife. 2020 Aug 24;9:e55792.\u003c/li\u003e\n\u003cli\u003eSoneson C, Love MI, Robinson MD. Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences [Internet]. F1000Research; 2016 [cited 2025 Sept 26]. Available from: https://f1000research.com/articles/4-1521\u003c/li\u003e\n\u003cli\u003eWang D, Liu Y, Zhang Y, Chen Q, Han Y, Hou W, et al. A real-world multi-center RNA-seq benchmarking study using the Quartet and MAQC reference materials. Nat Commun. 2024 July 22;15(1):6167.\u003c/li\u003e\n\u003cli\u003eHashimoto T, de Hoon MJL, Grimmond SM, Daub CO, Hayashizaki Y, Faulkner GJ. Probabilistic resolution of multi-mapping reads in massively parallel sequencing data using MuMRescueLite. Bioinforma Oxf Engl. 2009 Oct 1;25(19):2613\u0026ndash;4.\u003c/li\u003e\n\u003cli\u003eZytnicki M. mmquant: how to count multi-mapping reads? BMC Bioinformatics. 2017 Sept 15;18(1):411.\u003c/li\u003e\n\u003cli\u003eAxtell MJ. ShortStack: comprehensive annotation and quantification of small RNA genes. RNA N Y N. 2013 June;19(6):740\u0026ndash;51.\u003c/li\u003e\n\u003cli\u003eMorrissey A, Shi J, James DQ, Mahony S. Accurate allocation of multimapped reads enables regulatory element analysis at repeats. Genome Res. 2024 July 23;34(6):937\u0026ndash;51.\u003c/li\u003e\n\u003cli\u003eHita A, Brocart G, Fernandez A, Rehmsmeier M, Alemany A, Schvartzman S. MGcount: a total RNA-seq quantification tool to address multi-mapping and multi-overlapping alignments ambiguity in non-coding transcripts. BMC Bioinformatics. 2022 Dec;23(1):1\u0026ndash;21.\u003c/li\u003e\n\u003cli\u003ePardo-Palacios FJ, Wang D, Reese F, Diekhans M, Carbonell-Sala S, Williams B, et al. Systematic assessment of long-read RNA-seq methods for transcript identification and quantification. Nat Methods. 2024 July;21(7):1349\u0026ndash;63.\u003c/li\u003e\n\u003cli\u003eJi HJ, Pertea M. Enhancing transcriptome expression quantification through accurate assignment of long RNA sequencing reads with TranSigner. Genome Biol. 2025 Dec;26(1):1\u0026ndash;26.\u003c/li\u003e\n\u003cli\u003eZheng GXY, Terry JM, Belgrader P, Ryvkin P, Bent ZW, Wilson R, et al. Massively parallel digital transcriptional profiling of single cells. Nat Commun. 2017 Jan 16;8(1):14049.\u003c/li\u003e\n\u003cli\u003eKaminow B, Yunusov D, Dobin A. STARsolo: accurate, fast and versatile mapping/quantification of single-cell and single-nucleus RNA-seq data [Internet]. bioRxiv; 2021 [cited 2025 Sept 30]. p. 2021.05.05.442755. Available from: https://www.biorxiv.org/content/10.1101/2021.05.05.442755v1\u003c/li\u003e\n\u003cli\u003eSrivastava A, Malik L, Smith T, Sudbery I, Patro R. Alevin efficiently estimates accurate gene abundances from dscRNA-seq data. Genome Biol. 2019 Dec;20(1):1\u0026ndash;16.\u003c/li\u003e\n\u003cli\u003eMelsted P, Booeshaghi AS, Liu L, Gao F, Lu L, Min KH (Joseph), et al. Modular, efficient and constant-memory single-cell RNA-seq preprocessing. Nat Biotechnol. 2021 July;39(7):813\u0026ndash;8.\u003c/li\u003e\n\u003cli\u003eHoyos Sanchez MC, Ospina Zapata HS, Suarez BD, Ospina C, Barbosa HJ, Carranza Martinez JC, et al. A phased genome assembly of a Colombian Trypanosoma cruzi TcI strain and the evolution of gene families. Sci Rep. 2024 Jan 24;14(1):2054.\u003c/li\u003e\n\u003cli\u003eRhie A, McCarthy SA, Fedrigo O, Damas J, Formenti G, Koren S, et al. Towards complete and error-free genome assemblies of all vertebrate species. Nature. 2021 Apr;592(7856):737\u0026ndash;46.\u003c/li\u003e\n\u003cli\u003eHickey JW, Becker WR, Nevins SA, Horning A, Perez AE, Zhu C, et al. Organization of the human intestine at single-cell resolution. Nature. 2023 July;619(7970):572\u0026ndash;84.\u003c/li\u003e\n\u003cli\u003eThe Computational Pan-Genomics Consortium. Computational pan-genomics: status, promises and challenges. Brief Bioinform. 2018 Jan 1;19(1):118\u0026ndash;35.\u003c/li\u003e\n\u003cli\u003eLove MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014 Dec;15(12):1\u0026ndash;21.\u003c/li\u003e\n\u003cli\u003eAndrews S. FastQC: A Quality Control Tool for High Throughput Sequence Data. 2010; Available from: [Online]. Available online at: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/\u003c/li\u003e\n\u003cli\u003eJoshi NA, Fass JN. Sickle: a sliding-window, adaptive, quality-based trimming tool for FastQ files. 2011.\u003c/li\u003e\n\u003cli\u003eWang W, Peng D, Baptista RP, Li Y, Kissinger JC, Tarleton RL. Strain-specific genome evolution in Trypanosoma cruzi, the agent of Chagas disease. PLOS Pathog. 2021 Jan 28;17(1):e1009254.\u003c/li\u003e\n\u003cli\u003eFrazee AC, Jaffe AE, Langmead B, Leek JT. Polyester: simulating RNA-seq datasets with differential transcript expression. Bioinformatics. 2015 Sept 1;31(17):2778\u0026ndash;84.\u003c/li\u003e\n\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":true,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"Multi-mapping, RNA-seq, Salmon, STAR, Bowtie2, Kallisto, Trypanosoma cruzi, multigene family","lastPublishedDoi":"10.21203/rs.3.rs-7888056/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-7888056/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003e\u003cb\u003eBackground:\u003c/b\u003e\u003c/p\u003e\u003cp\u003eRNA sequencing (RNA-seq) enables transcript quantification and isoform analysis in diverse biological contexts, but accurately measuring expression from highly related genomic regions remains challenging. Multi-mapped reads\u0026mdash;those aligning equally well to multiple loci\u0026mdash;pose a major computational hurdle and compromise the overall accuracy of transcriptome resolution.\u003c/p\u003e\u003cp\u003e\u003cb\u003eResults:\u003c/b\u003e\u003c/p\u003e\u003cp\u003eWe herein evaluated five RNA-seq pipelines\u0026mdash;Bowtie2\u0026thinsp;+\u0026thinsp;featureCounts, STAR\u0026thinsp;+\u0026thinsp;featureCounts, STAR\u0026thinsp;+\u0026thinsp;Salmon, Salmon, and Kallisto\u0026mdash;on their ability to quantify gene expression in \u003cem\u003eTrypanosoma cruzi\u003c/em\u003e, a parasitic protozoan with a highly repetitive genome characterized by the abundance of large multigene families.Using real RNA-seq data, we first compared gene-level outputs, with emphasis on multigene family representation. Simulated transcriptomes were used to benchmark quantification accuracy under controlled conditions. Among the best-performing strategies (Salmon, Kallisto, and STAR\u0026thinsp;+\u0026thinsp;Salmon), we further tested whether including untranslated regions (UTRs) in gene annotations improved the assignment of ambiguous reads.\u003c/p\u003e\u003cp\u003e\u003cb\u003eConclusions:\u003c/b\u003e\u003c/p\u003e\u003cp\u003eOverall, the alignment-free transcriptome quantifiers Salmon and Kallisto achieved the most accurate performance, closely matching simulated values. Incorporating UTR annotations improved read assignment accuracy, particularly for STAR\u0026thinsp;+\u0026thinsp;Salmon. These tools not only enable global expression quantification but also facilitate precise read allocation between members of the same gene family, with up to 98% sequence identity. Our results highlight the critical role of annotation quality and quantification strategy in improving gene expression estimates.\u003c/p\u003e","manuscriptTitle":"The battle for reads: evaluating strategies to tackle multi-mapping in RNA-seq quantification in highly repetitive genomes","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2025-11-06 11:08:39","doi":"10.21203/rs.3.rs-7888056/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"0b8f64cd-0647-4131-8fa3-76f0ae987b75","owner":[],"postedDate":"November 6th, 2025","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[],"tags":[],"updatedAt":"2025-12-08T16:23:25+00:00","versionOfRecord":[],"versionCreatedAt":"2025-11-06 11:08:39","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-7888056","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-7888056","identity":"rs-7888056","version":["v1"]},"buildId":"8U1c8b4HqxoKbykW_rLl7","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}