The importance of data transformation in RNA-Seq preprocessing for bladder cancer subtyping

Abstract

Objective

RNA-Seq provides an accurate quantification of gene expression levels and it is widely used for molecular subtype classification in cancer, with special importance in prognosis. However, the reliability and validity of these analyses can significantly be influenced by how data are processed. In this study we evaluate how RNA-Seq preprocessing methods influence molecular subtype classification in bladder cancer. By benchmarking various aligners, quantifiers and methods of normalization and transformation, we stress the importance of preprocessing choices for accurate and consistent subtype classification.

Results

Our findings highlight that log-transformation plays a crucial role in centroid-based classifiers such as consensusMIBC and TCGAclas, while distribution-free algorithms like LundTax offer robustness to preprocessing variations. Non log-transformed data resulted in low classification rates and poor agreement with reference classifications in consensusMIBC and TCGAclas classifiers. Additionally, LundTax consistently demonstrated better separation among subtypes, compared to consensusMIBC and TCGAclas, regardless of preprocessing methods. Nonetheless, the study is limited by the lack of a true reference for objective assessment of the accuracy of the assigned subtypes. Hence, future work will be necessary to determine the robustness and scalability of the obtained results.

Introduction

While recent advancements in technologies, such as single-cell RNA sequencing and multi-omics approaches, have facilitated novel cancer classifications [1,2,3], RNA sequencing (RNA-Seq) remains a powerful and cost-effective tool that enables comprehensive transcriptome analysis by providing an accurate quantification of gene expression. RNA-Seq data are largely used for molecular subtype classification of diseases, crucial for improving prognostic accuracy and enhancing understanding of cancer biology [4,5,6,7,8,9]. In bladder cancer (BC), subtypes have shown to correlate with prognosis and treatment response [10,11,12]. Pioneering efforts by Robertson et al. defined five molecular subtypes and developed the TCGA classification (TCGAclas) [10], followed by the LundTax classifier developed by Sjöhdal et al. which focused on tumor differentiation and immune response [11]. Building on these efforts, Kamoun et al. unified previously reported molecular taxonomies and developed ConsensusMIBC, a robust classification for BC [12]. However, RNA-Seq preprocessing choices significantly impact the reliability of downstream analyses, including molecular subtype classification [13].

RNA-Seq data preprocessing workflow encompassess several steps to transform raw reads into a data matrix with the gene expression levels. The first step is the alignment of the reads to a reference genome or transcriptome using tools like STAR [14] or Hisat2 [15] among others. STAR typically offers higher accuracy in aligning complex and repetitive regions of the genome, while Hisat2 tends to be more precise in detecting single nucleotide polymorphisms. The next step is gene expression quantification, where the number of reads aligned are quantified to determine the gene expression levels. Widely used methods such as featureCounts [16] and HTSeq [17] perform straightforward read-counting methods, whereas StringTie [18] employs an EM algorithm to estimate transcript abundances, which helps in resolving ambiguities when the reads map to multiple transcripts. Although alignment and quantification are typically performed as separate steps, they can be performed together using pseudoaligners like Kallisto [19] and Salmon [20].

Gene expression levels are affected by biological and technical variability (e.g. sample quality, sequencing depth, batch effect). Normalization and transformation, such as transcript per million (TPM) or Trimmed Mean of M-values (TMM) [21], adjust raw counts to make expression values comparable. Log transformation is also essential as it balances skewed data and stabilizes the variance reducing the impact of outliers.

The main objective of this study is to evaluate the impact of different RNA-Seq preprocessing tools methods on molecular subtype classification in BC. We foresee that by ensuring the best choice for every preprocessing step, we can maximize the reliability of subtype classification [22].

Methods

We evaluated twelve combinations of preprocessing methods on three molecular subtype classifiers using four bladder cancer datasets: GSE148079 (n = 8) [23], GSE179440 (n = 12) [24], GSE186610 (n = 100) [25] and TCGA (n = 378) [10] (Fig. 1A). The preprocessing workflow included quality control (Supplementary Fig. 1), alignment (STAR, Hisat2), quantification (FeatureCounts, StringTie, HTSeq), and pseudoalignment (Salmon, Kallisto), using hg38 annotations. Different normalization and/or transformation methods (TMM, TPM and log2TPM), together with rawData and log2rawData were studied. ConsensusMIBC, LundTax, and TCGAclas classifiers were assessed using separation and coincidence metrics. (Supplementary Methods).

Fig. 1

Overview of methodology and count summaries of benchmarked tools. A Combination of datasets and preprocessing methods benchmarked. B Total counts for each dataset. Colors indicate distinct preprocessing methods for alignment and quantification. X axis shows the number of detected genes and Y axis shows the number of total counts

Results

The number of total counts and detected genes were consistent among the different alignment and quantification methods across datasets. STAR and Hisat2 outperformed the two pseudoaligners (Kallisto and Salmon) in the number of counts. Even though both pseudoaligners yielded the lowest number of counts, they detected an equivalent number of genes as the StringTie quantifier, regardless of the aligner. In contrast, featureCounts, followed by HTSeq, consistently detected the highest number of genes, except for the GSE179440 dataset, in which StringTie detected more genes than HTSeq. The number of counts was proportional to the number of sequenced reads, being the highest for the TCGA (average of 62 million reads) (Fig. 1B, Supplementary Table 1 and Supplementary Methods). Overall, the best performing methods across datasets were STAR or Hisat2 combined with featureCounts since they retrieve the highest number of genes (Fig. 1B, Supplementary Table 1).

We next evaluated if the differences observed in the number of counts and number of detected genes have an impact on BC subtyping. Our results on the consensusMIBC classifier [12] showed that, across all methods and datasets, using non log-transformed data (rawData or TPM) resulted in low correlation values and many unclassified samples (up to 87.5–100% in the two smallest datasets and 34.4%-64% in the largest). Even when few or no samples were unclassified (eg. featureCounts), correlation values were consistently higher for log transformed data, being log2TPM and TMM the highest. Of note, HTSeq and StringTie quantifiers were more affected by this issue, regardless of the aligner and specially in combination with TMM normalization (e.g. 0% vs 1.06–34.4% of unclassified samples, using log-transformed vs non log-transformed data in TCGA dataset). Similar results were observed with the pseudoaligners (0% vs 0–16.7%) (Fig. 2, Supplementary Table 2).

Fig. 2

Heatmap of consensusMIBC classification using benchmarked tools. Each column shows the results for the different RNA-Seq combinations of aligner, quantifier, of normalization and transformation method. Each row represents a sample from several public datasets (GSE148079, GSE179440, GSE186610 and TCGA). Colors of the heatmap indicate the predicted molecular subtype. Column on the left represents the reference for each sample (i.e. the most commonly predicted subtype). Luminal papillary (LumP), luminal unstable (LumU), luminal non-specified (LumNS), stroma rich (Stroma-rich), basal squamous (Ba/Sq) and neuroendocrine-like (NE-like). Kallisto, Salmon, STAR + featureCounts (STAR + featC), STAR + HTSeq, STAR + StringTie (STAR + StringT), Hisat2 + featureCounts (Hisat2 + featC), Hisat2 + HTSeq, Hisat2 + StringTie (Hisat2 + StringT) and Reference. Barplots show the mean and standard deviation of the correlation values for each methodology

Additionally, we assessed two other BC classifiers: TCGAclas [10] and LundTax [26]. The results for the TCGAclas classifier showed a high variability across methods, especially when comparing log-transformed and non log-transformed data (Supplementary Fig. 2). In contrast, LundTax appeared very stable, even across log and non-log data (Supplementary Fig. 3). For instance, the TCGAclas on TPM compared to log2TPM predicted 53% vs 18% of basal squamous, 37% vs 68% of luminal papillary (GSE186610, STAR + featureCounts). In contrast, the subtype classification of the same dataset using consensusMIBC and LundTax classifiers were more consistent across both pipelines, with variations of at most 3%. These results were observed in most of the different preprocessing pipelines that were applied, such as Hisat2 + featureCounts, Kallisto or STAR + HTSeq (Supplementary Figs. 2, 3 and Supplementary Tables 3, 4).

In an effort to fairly compare across the three classifiers, we evaluated the distribution of their scores (Supplementary Methods). The distribution of the non log-transformed data was distinctly separated from that of the log-transformed for consensusMIBC and TCGAclas classifiers. Importantly, all scores below the minimum confidence threshold belonged to non log data (Fig. 3A). In contrast, LundTax scores distribution was not influenced by log transformation (Supplementary Tables 5 and 6).

Fig. 3

Score distribution across classifiers and molecular subtypes interconnections of the consensusMIBC classifier. A Density plot of the score values distribution for each classifier, splitted by log-transformed and non log-transformed values. Shadows in each density plot represent the 95% confidence intervals. B Directed graph with reclassification percentages. Each node represents a consensusMIBC molecular subtype. The size shows the number of samples classified as that specific molecular subtype across methods. The distance between nodes and the width of the links represent the proportion of reclassified samples across methods between the subtypes. Edge colors are inherited by the original node. Edges below 1% were filtered out

To assess the performance, two metrics were evaluated: separation and coincidence. Separation shows how a sample is representative of its subtype whereas coincidence is the percentage of samples corresponding to the most frequent subtype (Supplementary Methods, Supplementary Fig. 4). The most important differences were influenced by the classifier and the normalization and transformation methods. Specifically, TPM normalization combined with HTSeq or StringTie exhibited the lowest coincidence scores across both the consensusMIBC (0.22–0.40) and TCGAclas (0.46–0.66) classifiers. LundTax was not influenced by log transformation and showed more stable metrics. Finally, pseudoaligners, when combined with log transformed data, performed similarly to other quantifiers, regardless of the aligner (Supplementary Fig. 4, Supplementary Table 7).

Of note, consensusMIBC and TCGAclas classifiers demonstrated low separation values regardless of the preprocessing methods used (0.1–0.32), indicating that the samples were less representative and less distinctly separated from other molecular subtypes. Conversely, the LundTax classifier achieved consistently the highest separation values (0.45–0.63) across different methods (Supplementary Fig. 4, Supplementary Table 7).

As observed in Fig. 2, the stability of the molecular subtype classification was highly sample-dependent. We studied the changes in molecular subtype classification across datasets and methods to analyze their exchangeability. The most frequent subtypes, Ba/Sq and LumP, were the most stable ones (93.8%; 90.3%), followed by the less frequent, NE-like (89.7%). As expected, we observed high exchangeability among luminal subtypes (LumP, LumU, and LumNS). Around 11.4% of LumNS and 13.2% of LumU were also classified as LumP, the most abundant luminal subtype. In contrast, the likelihood of a luminal subtype being classified as a non-luminal subtype was low or zero, as in the case of NE-like (Fig. 3B, Supplementary Table 8). The likelihood of interchange among non-luminal subtypes (Stroma-rich, Ba/Sq, and NE-Like) was not as high as among luminal subtypes (0–7.22%), but it was still higher than being reclassified as any of the luminal subtypes. Stroma-rich showed similar reclassification proportions with most subtypes (Fig. 3B, Supplementary Table 8).

Conclusions

Evaluating the impact of RNA-Seq preprocessing tools is essential for standardizing RNA-Seq pipelines, and thus maximizing the reliability of molecular subtype classification in bladder cancer. Our findings suggest that preprocessing steps including read alignment, gene quantification, normalization and data transformation, have an impact on the downstream analysis of molecular subtypes in bladder cancer. Despite some preprocessing methodologies being more affected by normalization and transformation, results showed that the critical factor when using consensusMIBC [12] and TCGAclas [10] classifiers is the log transformation of the data. Both classifiers are centroid-based and therefore highly sensitive to the distribution of the data. Indeed, Kamoun et al. recommend using log-transformed data to reduce common issues such as outliers or skewness [9]. The LundTax classifier [26], in contrast, is less affected by log-transformation. This may be attributed to the use of a Random Forest algorithm, which relies on relative comparisons rather than absolute values [27,28,29].

We examined the interconnections among molecular subtypes within the consensusMIBC classification to assess their exchangeability. As expected and previously reported, luminal subtypes showed a higher likelihood of being reclassified within the luminal group compared to non-luminal subtypes [28]. However, reclassification between luminal and non-luminal subtypes is not likely to happen, suggesting a distinct separation between these groups greater than differences across preprocessing methods [27, 30].

Based on our comprehensive evaluation and previous studies, both pseudoaligners (Kallisto [19] and Salmon [20]) and the featureCounts [16] quantifier combined with any aligner (STAR [14] or Hisat2 [15]) are effective for molecular subtype classification. While featureCounts provides better sensitivity in detecting lowly expressed genes, Kallisto and Salmon offer notable advantages in speed and memory efficiency [31]. Our findings suggest that, although Kallisto and Salmon identify a smaller number of genes, this does not compromise the accuracy of subtype classification, likely because the classifiers rely only on a subset of genes.

Moreover, we evaluated the classification rates to compare among the results obtained by the three classifiers. Our results showed low classification rates for HTSeq [17] and StringTie [18] quantifiers, regardless of the aligner, across all datasets and classifiers, indicating a low-accuracy in subtype assignment. The highest classification rates were achieved by featureCounts combined with STAR or Hisat2 aligners, along with Kallisto and Salmon.

We have tested gene-based quantifiers (featureCounts and HTSeq), as well as transcript quantifiers (StringTie), with the objective of assessing subtyping classification for which most methods are based on gene estimations. Despite gene-level results being often more accurate, powerful and interpretable than transcript-level results, difference between transcript-based and exon-based quantifiers is expected to be relatively minor when performing molecular subtype classification compared to other analyses, such as differential expression analysis, where gene quantification plays a crucial role [32,33,34].

In conclusion, our results show that log-transformation is a required step for centroid-based classifiers such as consensusMIBC and TCGAclas, that were trained on log-transformed data. In addition, we recommend distribution-free algorithms such as LundTax, which show less sensitivity to preprocessing steps. Future research should be focused on the validation of the robustness and scalability from the findings of our study in MIBC as well as on other cancer type molecular classifiers. Having standardized workflows could improve the application of molecular subtype classification across cancers, enabling a more thorough study of their association with treatment response and promoting the development of personalized therapeutic strategies in clinical practice.

Limitations

The main limitation of this study is the lack of a true reference to compute objective evaluation metrics, which makes it challenging to assess the accuracy and reliability of assigned subtypes. Coincidence can also be misleading as it measures the agreement of an assigned subtype with the most frequent subtype, without considering subtype designation accuracy.

Additionally, relying only on limited methods and datasets might introduce bias and miss important insights, compromising the reliability of the findings. Nevertheless, our findings reveal parallel constraints in predictive models development, which are heavily reliant on the datasets and methods used, particularly when external validation is not feasible.