Optimizing bioinformatic workflows to extract clinically usable gene expression data from targeted RNA sequencing panels: comparison with total RNAseq

Optimizing bioinformatic workflows to extract clinically usable gene expression data from targeted RNA sequencing panels: comparison with total RNAseq | 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 Article Optimizing bioinformatic workflows to extract clinically usable gene expression data from targeted RNA sequencing panels: comparison with total RNAseq Xiaokang Pan, Ashley Patton, Yi Seok Chang, Ryan Stevens, Nehad Mohamed, and 7 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-8695099/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 Targeted RNA sequencing (RNAseq) is widely used to detect gene fusions in tumors but clinical diagnostic use of expression data from these panels in fusion-negative cases has been limited. To facilitate this application, we evaluated methods for sequence read counting and gene normalization to optimize them for smaller gene sets. We present comparative methods to derive differential gene expression (DGE) data using ~ 200-gene clinically validated RNAseq fusion panels and compared them to parallel full RNAseq. We compared five methods for read counting, demonstrating that featureCounts is the most rapid and robust. For normalization prior to DGE with DESeq2, we compared five different normalization strategies and showed normalization using the 5 most stably expressed genes provided optimal centralization for these smaller gene sets. DGE output was assessed by principal component analysis (PCA), t-SNE and heatmap-clustering. The final pipeline was validated using PCA and pathway analysis by comparison with full RNAseq separately performed on a common set of challenging tumors with comparable results observed with the targeted gene panel. Overall, we show using an optimized bioinformatic pipeline that usable gene expression data can be obtained from smaller targeted RNAseq panels to maximize the clinical utility of these assays. Biological sciences/Biological techniques Health sciences/Biomarkers Biological sciences/Cancer Biological sciences/Computational biology and bioinformatics Biological sciences/Genetics Health sciences/Oncology Read counting normalization clustering tumor grading cell lineage assessment Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Introduction RNA sequencing (RNAseq) by next-generation sequencing (NGS) is commonly used to detect diagnostic or therapy-related gene fusions in formalin-fixed paraffin-embedded (FFPE) tumor samples. To achieve rapid signout and maximal sensitivity through deep coverage, clinical laboratories typically perform fusion detection using limited/targeted gene sets of 50–500 genes rather than assessing the entire transcriptome. Methods for targeted enrichment include bait-probe hybridization or anchored PCR amplification [ 1 , 2 , 3 ]. However, since relevant gene fusions are present in only a minority of tumors, many targeted RNAseq studies do not advance the diagnostic workup. Combining gene fusion detection with differential gene expression (DGE) analysis can increase the value of such panels. Clinically relevant goals of DGE in tumors include distinguishing site of origin (eg. bladder versus lung), distinguishing lower grade from higher grade soft tissue tumors and particularly the classification of poorly differentiated tumors where the differential includes carcinoma, melanoma and sarcoma. The first step in RNAseq-based DGE analysis is read counting where the number of sequence reads that align to each gene segment are quantified. With bait-probe capture and bridge amplification (Illumina) sequencing, NGS read counts have been shown to correlate with expression as determined by microarray or real-time PCR analysis [ 4 , 5 ]. Although bioinformatics tools for extracting valid tumor gene expression from whole transcriptome data have been well-established [ 6 ], methods optimized for small/targeted RNAseq data sets are currently limited. Furthermore, the accuracy and reproducibility of any given read counting method needs to be established for a fully validated DGE assay to be used clinically [ 4 ]. Given the wide variation in the quality of FFPE samples, fixation artifacts and tumor dilutional effects on differential gene expression can be significant, especially as compared to fresh/frozen tumor tissues [ 7 ]. Normalization or centralization methods for RNAseq DGE are usually applied to reduce bias due to variable tumor content and/or RNA degradation. Many normalization methods are available and have been evaluated for whole transcriptome RNAseq data [ 8 , 9 , 10 , 11 ]. Whether these studies are valid or applicable to limited gene RNAseq panels has not been well-elucidated. For ease of clinical signout, data presentation using diagnostically meaningful outputs that accurately represent differences between tumors are critical for clinical utility. For example, heatmaps display comparisons of the entre sample and are most useful for classification discovery with static datasets but are difficult to apply dynamically. Principal component analysis (PCA), a linear dimensionality reduction technique that preserves large pairwise distances and variance in the data, may more easily highlight clustering of individual samples without explicitly specifying distinguishing gene sets. Similarly, t-distributed stochastic neighbor embedding (t-SNE) is another data reduction technique that can potentially easily highlight similarity of a new sample to the reference set of different tumors. In this study, we employed clinically validated ~ 200-gene RNAseq panels optimized for fusion detection in solid tumors to systemically consider optimal methods of read counting, normalization and data visualization. We compared the outputs of the optimized pipelines for a targeted panel that included a limited number of genes included cell lineage and tumor grade assessment with parallel full RNAseq and showed highly comparable results in a set of poorly differentiated tumors. Results Optimizing read counting To convert NGS output from the targeted RNAseq panels into gene expression data, raw sequence reads obtained from each gene segment were counted for each sample using 6 different tools: Samtools-view (referent clinical assay method), featureCounts, HTSeq-count, coverageBed, cuffdiff and countOverlaps. A sample output from a 16-sample RNAseq fusion assay is shown comparing usability, time-to-output and resource use. Using 8 CPU processors, featureCounts, took only 2 minutes to complete as compared to coverageBed at 9 minutes ( Table S3 ). In contrast, HTSeq-count, CuffDiff and Samtools-view produced outputs in 3, 4 and 6 hours, respectively. The total reads in the outputs of Samtools-view, featureCounts, HTSeq-count and coverageBed were similar for 10 samples (Table 1 ), whereas the read counts from cuffdiff outputs were much smaller due to larger genes producing high read counts being skipped in the calculation by the software. Excluding cuffdiff, the CV values for each of the 222 expressed genes were computed for each program. The median CV ranks of the 222 expressed genes from the four programs were also nearly identical (Fig. 1 a and Table 2 ). The correlation of each expressed gene in the panel was compared with the validated Samtools-view value using Pearson correlation coefficients and showed high correlations (Fig. 1 b); the median rank of the correlation coefficients showed only small differences (Table 2 ). Considering the performance metrics of output time, precision and accuracy, featureCounts was rated as the best read counting method for this application. Table 1 Number of read counts in 10-sample tun using different read counting methods (assay 2). Method S-107 S-143 S-153 S-175 S-359 S-394 S-458 S-504 S-525 S-614 Samtools-view 23863915 84929701 38449531 103286313 61451551 52339994 173011118 32439398 60600048 72260163 FeatureCounts 23559590 84701209 37895291 102612244 61105419 51932562 171937331 32037416 60076885 71446450 HTSeq-count 23505495 84600839 37801803 102395919 60991722 51824717 171667472 31965039 59889358 71322334 CoverageBED 24306691 85794224 39043406 104359380 62160975 53154627 174990470 32974310 61740756 73204253 CuffDiff 3863910 3428790 62443716 7726710 4301250 5106000 9124866 4543470 9254736 10259952 Table 2 Precision and accuracy performance of different read counting methods. Rank Counting method Median precision Median accuracy Overall 1 Samtools-view 110 na na 2 coverageBed 110 93 203 3 featureCounts 110 111 221 4 HTSeq 110 113 223 NA: Not applicable as Samtools-view is the referent method. Effect of expression normalization strategies on tumor grouping and data centralization Clinical biopsies/resections show a range of tumor cellularity and varying admixture of neoplastic and non-neoplastic cell types. Therefore, minimizing these effects on overall gene expression distribution through read count normalization is critical in highlighting diagnostically relevant sample differences. Especially in suboptimal FFPE samples, different normalization strategies are expected to have large effects. Using the raw read counts obtained with featureCounts as input (“NoNorm”), we compared the impact of normalization with ACTB , TraditionHK5, IdentifiedHK5, Top10H, Top10S and AllGenes (see Methods for definitions of these groups) followed by gene expression analysis in DESeq2. One readout of effective normalization assessed was a reduction in intragroup variation across tumors of similar types. This was assessed by calculating the average CV percentage values of read counts before or after normalization for groups of carcinomas, low-grade soft tissue tumors and high-grade sarcomas. The lowest CV values for each of the 3 tumor groups were achieved with the IdentifiedHK5, Top10S and AllGenes normalization methods (Table 3 ). The intragroup variation for these methods was significantly reduced compared to the NoNorm control condition (P < 0.01, Wilcoxon signed rank test). Intra-tumor group variation was higher than NoNorm for the TraditionHK5 and ACTB normalization methods, with the latter having significantly higher average CV values for both the carcinoma and sarcoma groups. Table 3 Average CV percentage values after normalization by different methods using DESeq2. Group NoNorm ACTB TraditionHK5 IdentifiedHK5 Top10H Top10S AllGenes Carcinomas 93.17 93.99 88.45* 84.99** 104.65** 86.67** 85.40** Soft Tissue Tumors 82.93 85.30* 84.20 83.35 89.09** 82.57 82.16 Low-grade Sarcoma 100.15 96.67 84.36** 82.53** 101.60 83.33** 81.54** High-grade Sarcoma 119.15 127.86** 121.65 105.67** 141.57** 105.44** 105.76** Note: Comparisons via “estimateSizeFactors” function in DESeq2 across different number of genes as control. ** and * represent 0.01 and 0.05 significantly different from the value without normalization NoNorm; red and black stars indicate significantly higher and lower, respectively. Using DESeq2 to compare gene expression in distinct tumor groups, the distributions of all the logFC values were calculated for downregulated and upregulated genes in carcinoma as compared to soft tissue tumors in assay 2 (Fig. 2 a). The distribution of the logFC values of upregulated genes was highly biased without normalization likely influenced by a small set of highly expressed set of genes related to extracellular matrix production such as collagen isoforms. Normalization using ACTB or multiple HK genes (TraditionHK5) exaggerated this effect. Normalization using Top10H better highlighted differences among down-regulated genes. Normalization using AllGenes and Top10S moved the distribution of the logFC values toward to 0 for upregulated genes, whereas the IdentifiedHK5 method centralized the distribution of the logFC values for both upregulated and downregulated genes. To assess the effects of normalization strategies on determining tumor grade by gene expression patterns, the distribution of the logFC values of the read counts was compared to a set of low-grade soft tissue tumors and high-grade sarcomas in assay 1 (Fig. 2 b). Without normalization, upregulated genes biased the distribution as did including all genes, whereas ACTB or the TraditionHK5 and Top10H methods improved the separation between the groups. IdentifiedHK5 and Top10S produced centralization of the distribution. Possibly due to effects on minimizing differences simply due to tumor cellularity, normalization using IdentifiedHK5 and Top10S produced excellent results for this assay indication. Effect of expression normalization strategies on tumor type clustering An important rationale for normalization strategies in targeted tumor RNAseq assays is to enhance separation of tumor type clusters to improve mapping of samples of unknown lineage. Here, we evaluated several clustering and visualization methods, including PCA, t-SNE and Heatmap-clustering algorithms for their ease of interpretation by several molecular pathologists among the authors. Using a comparison of carcinoma and soft tissue tumor groups, the effects of read counts normalized by different methods on clustering by PCA was assessed. The IdentifiedHK5 and Top10S methods enhanced the separation of PCA clusters with IdentifiedHK5 showing best results (Fig. 3 ). The other normalization methods reduced the separation of PCA clusters as compared to no normalization. Although a heatmap of the same data did highlight the tumor groups, the complicated display did not easily facilitate expression pattern of new samples during a quick review (not shown). With the IdentifiedHK5 normalization method, there was no obvious clustering with t-SNE for the carcinoma versus soft tissue diagnostic indication (not shown). Similar results are seen for the low-grade and high-grade tumor comparison (not shown). Using DESeq2, the most differentially expressed genes comparing carcinoma and soft tissue tumors, the IdentifiedHK5 method again showed the best results with a few more distinguishing genes (Fig. 5 a). A similar pattern was seen with the low-grade versus high-grade tumor sets as a separate distinct diagnostic indication (Fig. 5 b). Comparison of targeted RNAseq outputs with total RNAseq To validate the DGE output of the targeted panel 2, we compared the outputs with those obtained by total RNAseq on a diagnostically relevant set of 32 poorly differentiated tumors (19 carcinoma, 13 sarcoma). Although the number of expressed genes across all samples was vast differently (226 versus 20218; Table 4 ), the significantly differentially expressed genes between carcinoma and sarcoma cases (at adjusted p < 0.05 level) were very similar for the shared genes ( Supplementary Table 6 ). The mapping of these cases into carcinoma and sarcoma clusters was also similar by PCA (Fig. 5 a). Similarly, the top differentially upregulated or downregulated canonical pathways by IPA analysis were also similar (Fig. 5 b). However, most of the disease states differentially regulated in sarcomatoid carcinoma and sarcoma, as identified by IPA, were distinct when exome data was compared to the targeted panel (Fig. 5 c). Although S100 family signaling was a shared upregulated pathway in sarcomas by both analyses, the targeted panel identified more specific signaling pathways in the carcinoma group. Table 4 Comparison of expressed genes in optimized outputs for targeted RNAseq (panel 2) and full RNAseq. Datasets Total genes studied Correlation coefficient for shared genes* Expressed genes Expressed genes (%) Expressed genes shared Total RNAseq 20218 0.95 2591 12.8 32 Targeted RNAseq 226 47 14.2 Note: Expressed genes selected by cutoffs of adjust p-value 0.05 and log2fc = 0.7. Discussion Using clinical-grade ~ 200-gene RNAseq assays originally developed for gene fusion detection, we have systematically evaluated their suitability for tumor lineage assessment and grading. The best methods for read counting (for pipeline optimization) and data normalization using DGE and PCA for evaluation, were determined. The output of the optimized targeted panel pipeline was then validated for a difficult application (distinguishing sarcomatoid poorly differentiated tumors representing either carcinoma or sarcoma) using total RNAseq. These comparisons provide a model for validating these methods for routine clinical use and highlight the importance of closely matching the methods to the design and goals of each assay. As highlighted by several use cases, the purposes for performing differential gene expression by RNAseq in tumors are varied. A common goal is to highlight one or more highly overexpressed lineage-associated genes in any given sample that may assist diagnosis. These uncommon expression patterns often signal cellular differentiation characteristics of a specific tumor type (eg high MYOD1 indicating skeletal muscle differentiation in rhabdomyosarcoma) [12{]. Another goal is to superimpose the overall gene expression pattern of any given tumor against a model set of different classes of tumors (e.g. carcinoma, melanoma and soft tissue tumors) to help determine site of origin when other biomarkers and microscopic studies are not informative. However, given the clinical impact of both indications, RNAseq analysis methods need to be carefully validated, and outputs need to be relatively easy to interpret. Targeted RNAseq often includes targets whose expression can vary widely based on specific molecular aberrations, which can cause some widely used software tools or methods to perform inaccurately. We therefore considered each step in the pipeline with regards to the clinical use cases. Many read counting algorithms have been developed and used for whole transcriptome RNAseq. As shown in Table S3, different read counting methods vary in speed performance. Our prior method, Samtools built-in program “view” calculates sequence reads in the targeted regions from alignment BAMs accurately when validated with the proper option settings [ 17 ]. But the long processing times can cause delays in clinical reporting. CuffDiff is another widely used software tool used for DGE analysis in whole transcriptome RNAseq [ 19 ]. It has a function to generate the number of sequence reads per gene in a sample. HTSeq-count is another widely used tool [ 18 ] BEDTools/coverageBed with the “-count” option is another option that allows multiple samples to be analyzed in parallel [ 20 ]. Corchete, et al [ 4 ] compared six methods of read counting for RNAseq data and reported HTSeq-count was superior. However, their study did not include coverageBed, Samtools-view or featureCounts which is optimized for efficient chromosome hashing and feature blocking techniques and parallel analysis. [ 17 ]. They also did not perform assessments of speed performance which is critical for clinical applications. Comparing speeding, precision accuracy and ease of use, we found that featureCounts and coverageBed had a significant advantage over other methods in output time, with featureCounts being the fastest. By statistical measures, featureCounts, coverageBed, HTSeq-count and Samtools-view have highly similar outputs. However, we found that CuffDiff underestimates sequence reads significantly which was traced to the limitations in the maximum number of reads in buffer of the software allowance which produces trimming of some reads from larger genes. Given that featureCounts and HTSeq-count require conversion of BED into GTF file, the ability of coverageBed to directly consume a BAM file is a great feature for smaller panels where the BAM files are not very large. When small gene panels are employed for DGE, minimizing skewing of the overall gene expression distribution (or centralization) is critical. To find the optimal normalization method of sequence reads for DGE from small RNAseq fusion panels, we compared multiple common employed strategies and a lab/assay-specific method accounting for the sample and extraction patterns in our laboratory. In the latter, we identified the five most consistently and stably expressed genes across the tumors typically analyzed by our RNAseq fusion panel, aka empirically determined mostly stable genes. We then compared normalization using these 5 genes to other normalization methods. This approach has impact on centralizing the distribution of log2FC values of DGE and improving or enhancing the separation of PCA clusters. In addition, this method could reduce intragroup variances significantly, as did the AllGenes and Top10S methods. The significantly expressed genes identified by DGE were also mostly similar with these three methods. Overall, for this assay/indication, normalization method with the empirically top 5 most stable genes was the superior normalization method. The traditional normalization methods using a single housekeeping gene such as ACTB and multiple highly expressed genes (e.g. Top10H) as references, respectively performed poorly in reducing intra-tumor group variation. Similarly, the default normalization method for DESeq2 (AllGenes) was not optimal for these targeted panels. These findings emphasize that assessment of different centralization/normalization methods, especially for more targeted RNAseq panel, should be included in the validation of each new assay with sample sets tuned for the specific clinical applications. Consideration of different data display methods is also critical given the limited time available for analysis for clinical assays. When mapping expression patterns of new/unknown tumors to existing data sets, we found that PCA is the best sample clustering method. By displaying typical tumor groups clearly (eg carcinoma vs sarcoma, high-grade vs low-grade), it can efficiently identify when an unknown sample maps well within a group as opposed to being an outlier where the study is non-informative. Although heatmaps are a useful visualization method when tumor expression patterns are highly similar (eg cultured tumor cells exposed to different drug doses), it is most difficult to visualize tumor-group differences in a single new sample. Using a carcinoma versus sarcoma tumor set, an optimized targeted panel with some lineage-specific genes included (panel 2) showed equivalent separation by PCA with parallel data from total RNAseq. In addition, many of the top tumor class discriminator were shared by the two assays. This similar result was further noted by IPA comparison of the top differentially regulated canonical signaling pathways (Fig. 5 ). However, disease states mapping by IPA among the carcinomas (red) and sarcomas (blue) was largely distinct between targeted panel and total RNAseq. This was likely due to the effects of genes not present in the targeted panel in providing subclass differentiation in small tumor sets. In summary, we identified featureCounts as the optimal read counting method, five assay-specific most stable genes as the best method for normalization prior to DESeq2 and PCA as the easiest to interpret sample clustering method for this targeted RNAseq panel clinical application. Methods NGS library design and sample sets The gene fusion NGS panels employed in this study included a custom 190-gene (panel 1) and 230-gene (panel 2) custom RNAseq panels used at The Ohio State University James Molecular Laboratory to detect diagnostically relevant gene fusions in human tumors. The 190-gene design included only a few typical housekeeping genes ( ACTB, MYH9 ) with no gene content explicitly designed for DGE. Over 500 clinical cases, the frequency of reportable gene fusions by panel 1 varied by diagnosis approaching 20% for soft tissue tumors versus ~ 5% for carcinomas. To improve utility for gene expression, panel 2 incorporated additional stably expressed/housekeeping and lineage-specific genes to aid in separation of poorly differentiated carcinomas and sarcoma. In 250 clinical cases using panel 2, 226 genes were routinely expressed at some level, with other genes only expressed when a particular fusion was present. The frequency of oncogenic fusion detection was similar to panel 1. For validation of the design of panel 2 using an optimized pipeline, full RNAseq (total transcriptome except for ribosomal genes) was performed for a set of diagnostically challenging spindled cell tumors. These represented 32 sarcomatoid poorly differentiated tumors that were diagnostically as either carcinoma or sarcoma following routine histopathology/immunohistochemistry workup. Final diagnosis was rendered by a soft tissue pathologist following comprehensive DNA mutation profiling and correlation with radiologic appearances and clinical presentation. All assays used and samples procured were obtained as part of routine clinical testing. The use of this data for development of the analytic pipeline was reviewed by the institutional review board with waiver of consent provided. All NGS protocols utilized total tumor RNA extracted from formalin-fixed paraffin-embedded tissue (FFPE) using PureLink FFPE RNA Isolation Kit (Invitrogen/ThermoFisher) and employed library preparation with DNA digestion and ribosomal RNA depletion (KAPA RNA HyperPrep Kit or Watchmaker Polaris). For the targeted panels, this was followed by hybridization with probes covering the full exonic regions of all target genes (xGen, IDT, Coralville, IA) designed through the GOAL consortium [ 14 ]. This size of the panel allowed 10–16 samples to be run on the Illumina NovaSeq SP flow cell, with adequate depth of coverage (~ 20–50,000 mean reads per sample per targeted area). To assess the suitability of methods across different gene sets, tumors sequenced with both panel 1 and panel 2 were included. For comparative analysis of performance of different read counting methods, 10 soft tissue tumor samples from panel 2 were used. For the comparison of normalization and clustering methods, 36 tumor samples from panel 2 were used, with 10 carcinomas and 11 soft tissue tumors as model comparators for DGE analysis ( Table S1 ). Initial clinical pipeline For clinical reporting, we employed a custom pipeline for fusion detection (“FindRNAFusion”): paired-end FASTQ files were downloaded to a high-performance HPE Linux server (256 CPU processors, 755 GB RAM memory, and RHEL 9.0 OS) from a mounted Illumina BaseSpace instance. The raw sequence reads of each sample were mapped to Hg19 (Human Genome version 19) to generate alignment BAM files using STAR [ 15 ]. Each BAM file was indexed by SAMTOOLS. The BAM files were then used to make fusion calls using Arriba [ 16 ]. After filtering out artifacts and low-level fusion calls (< 10–20 supporting reads), clinically reportable fusions were summarized in a text report file and a pdf file with graphical display. In this pipeline, the number of sequence reads/coverage in each targeted gene and in each targeted region/exon were also computed from the BAM files using each read counting method, as discussed below. These sequence reads were then normalized, using the methods described below. Subsequently, genes with very low reads were filtered out using a data-based threshold for maximum-based filters proposed by Rau et al [ 17 ] and the remaining genes with normalized read counts were used for DGE analysis using DESeq2 [ 18 ]. Read count method comparisons Five read counting methods, Samtools-view [ 19 ], HTSeq-count [ 20 ], featureCounts [ 21 ], coverageBed [ 22 ], and CuffLinks [ 23 ], were selected for comparative analysis. Samtools-view, which was our previously validated referent method, uses a BED file and a BAM file as input and outputs the number of reads in each targeted range (command line: “samtools view -F 0x04 -q 20 -c -@ 8”). This accuracy of this output was validated by comparing the count of mapped reads to manual inspection/calculation of read depth for a set of genes in the Integrative Genomics Viewer (IGV) [ 24 ] and proportion of fusion-positive tumor cells by in situ hybridization for selected common gene fusions. In contrast, featureCounts uses targeted genes in GTF file and BAM file as input (command line: “featureCounts -p -M -O -C -g gene_name --minOverlap 1 --maxMOp 30 -Q 20 –T 8”). HTSeq-count uses the same files as featureCounts as input (command line: “htseq-count -i gene_name --max-reads-in-buffer 50000000 -s no”). CuffDiff was run with option “--total-hits-norm TRUE” with the same files as featureCounts as input. coverageBed was run with option “-count” and “-a BED file” and “-b BAM file” as input. For this study, a Perl script was written to run program commands in parallel with a batch of NGS samples. Concordance of different methods was assessed among the 226 routinely expressed genes in panel 2 with the fully validated Samtools-view method as the referent. The precision of each method was assessed using the median rank from the coefficient of variation (CV) values of the 226 routinely expressed genes for each sample [ 25 ]. We also ranked CV values for each method. We also computed the median of the ranks of the 226 commonly expressed genes as the precision index of each method with lower values correlated with increased precision. The Pearson correlation coefficient (r) was computed to assess the association between the outputs of Samtools-view and the other methods. These correlation coefficients were calculated using the statistical functions in Microsoft Excel. Normalization comparisons DESeq2 is a widely used software tool for differential gene expression analysis and was chosen here as it performs well for larger RNAseq datasets (sample size > = 6) [ 16 ]. The standard normalization implemented in DESeq2 is relative log expression (RLE). The calculation of size factors is performed through the “estimateSizeFactor” function. The default for RLE normalization in DESeq2 is using all genes (“AllGenes”). We compared this result with those using a highly expressed “housekeeping” gene ACTB (HK), a panel of five traditionally used housekeeping genes, namely ACTB, MYH9, RANBP2, PRKACA and TFG (TraditionHK5) and to the top ten highly expressed genes as determined by the average number of reads in all the samples (Top10H) and the top 10 genes with smallest CV values in all the samples in each dataset (Top10S). Five genes ( CREBBP, BRAF, BRD4, ATF1 and CREB1 ) that were recurrently identified as among the top ten most stable in the range of test samples for panels 1 and 2 were also assessed (IdentifiedHK5, Table S2 ). These six normalization approaches were compared to no normalization for their effects on statistical and visualized centralization of the data. The CV% value representing the percentage of the standard deviation to the mean per gene across samples was computed using the formula CV% = STDEV.P * 100/AVERAGE. Total RNAseq data produced from the same RNA samples and library preparation as the targeted panels were analyzed using the core pipelines normalized using the 5 most stable genes from the targeted panel as normalization and subjected to DEQ2 as above. For pathway analysis, QIAGEN Ingenuity Pathway Analysis (IPA) was used ( https://www.qiagenbioinformatics.com ). Clustering and visualization The three clustering and data visualization methods employed are summarized in Table S4 . Principal component analysis (PCA), t-distributed stochastic neighbor embedding (t-SNE) and heatmap-clustering were selected for comparative analysis in sample clustering and gene expression. SRplot [ 26 ] was used to perform principal component analysis (PCA) and heatmap-clustering. t-SNE-Java ( https://github.com/lejon/T-SNE-Java ) was implemented to generate t-SNE graphics for clustering and visualization. Declarations Relevant Disclosures/Conflicts of Interest None. Funding statement: No relevant funding to declare. Author Contribution The Authors participated in the bioinformatic analysis (X.P., M.H., D.J.), data analysis (R.S., N.M., D.C., D.J.), performed the experiments (R.S., Y.C., D.C.) or provided molecular pathology review, data integration and case selection (D.J, Y.C., Y.H., C.M., W.Z., M.A.). D.J. and X.P. wrote the main manuscript. All authors reviewed the manuscript and provided feedback. Acknowledgement The Authors thank the medical technologists and the data scientists of the Polaris Molecular Laboratory for performing and analyzing the targeted RNAseq panels. The diagnostic work of soft tissue pathologists Hans Iwenofu and Swati Satturwar for the diagnostic challenging dataset is noted. Data Availability The data underlying this article will be shared on reasonable request to the corresponding author. References Curion, F. et al. Targeted RNA sequencing enhances gene expression profiling of ultra-low input samples. RNA Biol. 17 , 1741–1753 (2020). Heyer, E. E. et al. Diagnosis of fusion genes using targeted RNA sequencing. Nat. Commun. 10 , 1388 (2019). Capone, I. et al. Targeted RNA sequencing analysis for fusion transcript detection in tumour diagnostics: assessment of bioinformatic tools reliability in FFPE samples. Explor. Target. Antitumor Ther. 3 , 582–597 (2022). Corchete, L. A. et al. Systematic comparison and assessment of RNA-seq procedures for gene expression quantitative analysis. Sci. Rep. 10 , 19737 (2020). de Brito, M. W., de Carvalho, S. S., Mota, M. B. & Mesquita, R. D. RNA-seq validation: software for selection of reference and variable candidate genes for RT-qPCR. BMC Genom. 25 , 697 (2024). Fu, C. et al. Targeted RNA-seq assay incorporating unique molecular identifiers for improved quantification of gene expression signatures and transcribed mutation fraction in fixed tumour samples. BMC Cancer . 21 , 114 (2021). Li, J., Fu, C., Speed, T. P., Wang, W. & Symmans, W. F. Accurate RNA sequencing from formalin-fixed cancer tissue to represent high-quality transcriptome from frozen tissue. JCO Precision Oncol. 2 , 1–9 (2018). Välikangas, T., Suomi, T. & Elo, L. L. A systematic evaluation of normalization methods in quantitative label-free proteomics. Brief. Bioinform . 19 , 1–11 (2018). Zhao, Y. et al. TPM, FPKM or normalized counts? A comparative study of quantification measures for the analysis of RNA-seq data from the NCI Patient-Derived Models Repository. J. Transl Med. 19 , 269 (2021). Bushel, P. R. et al. Comparison of normalization methods for analysis of TempO-Seq targeted RNA sequencing data. Front. Genet. 11 , 594 (2020). Evans, C., Hardin, J. & Stoebel, D. M. Selecting between-sample RNA-seq normalization methods from the perspective of their assumptions. Brief. Bioinform . 19 , 776–792 (2018). Avenarius, M. R. et al. Integrated molecular profiling of rhabdomyosarcoma subtypes by targeted RNA-seq. medRxiv ( (2024). Amin, M. B. et al. AJCC Cancer Staging Manual, 8th ednSpringer, New York,. 13 (2017). Aisner, D. L. et al. The Genomics Organization for Academic Laboratories (GOAL): a vision for a genomics future for academic pathology. Acad. Pathol. 10 , 100090 (2023). Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29 , 15–21 (2013). Uhrig, S. et al. Accurate and efficient detection of gene fusions from RNA sequencing data. Genome Res. 31 , 448–460 (2021). Rau, A., Gallopin, M., Celeux, G. & Jaffrézic, F. Data-based filtering for replicated high-throughput transcriptome sequencing experiments. Bioinformatics 29 , 2146–2152 (2013). Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15 , 550 (2014). Danecek, P. et al. Twelve years of SAMtools and BCFtools. GigaScience. 10, giab008 (2021). Anders, S., Pyl, P. T. & Huber, W. HTSeq—a Python framework to work with high-throughput sequencing data. Bioinformatics 31 , 166–169 (2015). Liao, Y., Gordon, K., Smyth, G. K. & Shi, W. featureCounts: an efficient general-purpose program for assigning sequence reads to genomic features. Bioinformatics 30 , 923–930 (2014). Quinlan, A. R. & Hall, I. M. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26 , 841–842 (2010). Trapnell, C. et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat. Protoc. 7 , 562–578 (2012). Robinson, J. T. et al. Integrative Genomics Viewer. Nat. Biotechnol. 29 , 24–26 (2011). Li, D. et al. An evaluation of RNA-seq differential analysis methods. PLoS One . 17 , e0264246 (2022). Tang, D. et al. SRplot: a free online platform for data visualization and graphing. PLoS One . 18 , e0294236 (2023). Additional Declarations No competing interests reported. Supplementary Files PanetalSupplementaryTables012525.docx 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-8695099","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Article","associatedPublications":[],"authors":[{"id":581606712,"identity":"b2e5d10f-0b50-48ee-b159-f3b2e3b63185","order_by":0,"name":"Xiaokang Pan","email":"","orcid":"","institution":"The James Cancer Hospital","correspondingAuthor":false,"prefix":"","firstName":"Xiaokang","middleName":"","lastName":"Pan","suffix":""},{"id":581606713,"identity":"0b6a7870-3d87-4306-8b2d-07a4e4d304e2","order_by":1,"name":"Ashley Patton","email":"","orcid":"","institution":"The Ohio State University Wexner Medical Center","correspondingAuthor":false,"prefix":"","firstName":"Ashley","middleName":"","lastName":"Patton","suffix":""},{"id":581606714,"identity":"0b6a0e86-eb1c-4b66-adea-619290e83496","order_by":2,"name":"Yi Seok Chang","email":"","orcid":"","institution":"The James Cancer Hospital","correspondingAuthor":false,"prefix":"","firstName":"Yi","middleName":"Seok","lastName":"Chang","suffix":""},{"id":581606715,"identity":"1ab326cf-9aa3-4aab-8771-9c9e9653f327","order_by":3,"name":"Ryan Stevens","email":"","orcid":"","institution":"The James Cancer Hospital","correspondingAuthor":false,"prefix":"","firstName":"Ryan","middleName":"","lastName":"Stevens","suffix":""},{"id":581606716,"identity":"ffc39a7d-3016-458a-9216-c4b44b1ded28","order_by":4,"name":"Nehad Mohamed","email":"","orcid":"","institution":"The James Cancer Hospital","correspondingAuthor":false,"prefix":"","firstName":"Nehad","middleName":"","lastName":"Mohamed","suffix":""},{"id":581606717,"identity":"3179e43b-cf49-4fb3-9ee9-1a0466196bc3","order_by":5,"name":"Matthew Hunt","email":"","orcid":"","institution":"The James Cancer Hospital","correspondingAuthor":false,"prefix":"","firstName":"Matthew","middleName":"","lastName":"Hunt","suffix":""},{"id":581606718,"identity":"0537470c-14f2-44ce-ae51-87e5cc5a3a2c","order_by":6,"name":"Daniel Chappell","email":"","orcid":"","institution":"The James Cancer Hospital","correspondingAuthor":false,"prefix":"","firstName":"Daniel","middleName":"","lastName":"Chappell","suffix":""},{"id":581606719,"identity":"fc86adbb-2a84-4004-aeb3-7950b79cfac2","order_by":7,"name":"Yan Hu","email":"","orcid":"","institution":"The Ohio State University Wexner Medical Center","correspondingAuthor":false,"prefix":"","firstName":"Yan","middleName":"","lastName":"Hu","suffix":""},{"id":581606720,"identity":"6a9765c1-3a1a-4935-8a9f-f5fe818c47a8","order_by":8,"name":"Cecelia Miller","email":"","orcid":"","institution":"The Ohio State University Wexner Medical Center","correspondingAuthor":false,"prefix":"","firstName":"Cecelia","middleName":"","lastName":"Miller","suffix":""},{"id":581606721,"identity":"ef865feb-b0e4-489b-ae8d-682e423b7d41","order_by":9,"name":"Weiqiang Zhao","email":"","orcid":"","institution":"The Ohio State University Wexner Medical Center","correspondingAuthor":false,"prefix":"","firstName":"Weiqiang","middleName":"","lastName":"Zhao","suffix":""},{"id":581606722,"identity":"6fb78597-e3fb-4c45-aee7-dfa795ddc7b7","order_by":10,"name":"Matthew Avenarius","email":"","orcid":"","institution":"The Ohio State University Wexner Medical Center","correspondingAuthor":false,"prefix":"","firstName":"Matthew","middleName":"","lastName":"Avenarius","suffix":""},{"id":581606723,"identity":"e607ed25-4459-42bd-b1dc-c63881d15c63","order_by":11,"name":"Dan Jones","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA2UlEQVRIiWNgGAWjYBACxgYw9U/Ovr35AEToAHFaDhgb8BxLIE4LFBxINJDIMSBOC3N778GHPxjuJJjznPn48GcbgxzfjQQCDus5l2zMw/Asz7K9d7MxbxuDsSRBLTNyzKSBlhUznDm7TZqxjSFxAxFazH/+YGBObLiR80wS6LB6YrSYMfAwHAYansMmAXRYggFhv5wxluYxSDOW7DlmbMxzTsJw5pkH+LUYtvcYfvxRYSPHz9788OGPMht5vuMEbDFsAJEGcL4EfuUgIE9YySgYBaNgFIx4AADyqUd8mYANHQAAAABJRU5ErkJggg==","orcid":"","institution":"The Ohio State University Wexner Medical Center","correspondingAuthor":true,"prefix":"","firstName":"Dan","middleName":"","lastName":"Jones","suffix":""}],"badges":[],"createdAt":"2026-01-25 22:38:26","currentVersionCode":1,"declarations":{"humanSubjects":false,"vertebrateSubjects":false,"conflictsOfInterestStatement":false,"humanSubjectEthicalGuidelines":false,"humanSubjectConsent":false,"humanSubjectClinicalTrial":false,"humanSubjectCaseReport":false,"vertebrateSubjectEthicalGuidelines":false},"doi":"10.21203/rs.3.rs-8695099/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-8695099/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":101790568,"identity":"c3309a91-c669-434e-b904-038348d24322","added_by":"auto","created_at":"2026-02-03 16:06:20","extension":"jpg","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":63050,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003ePerformance of different sequence read counting methods for using a targeted RNAseq panel\u003c/strong\u003e. a. Ranks of CV values ie precision ranks for the 222 expressed genes in panel 2 are displayed for each read counting tool. b. The 5 other tools were compared to the validated reference method (Samtools-view) for concordance, with ranks of Pearson correlation coefficients for the 222 expressed genes for each method as compared to the referent counting method (ie accuracy ranks).\u003c/p\u003e","description":"","filename":"1.jpg","url":"https://assets-eu.researchsquare.com/files/rs-8695099/v1/1a63707aa9e143212e7ed5ab.jpg"},{"id":101790565,"identity":"688ff14a-2560-4aaf-ae7d-bb8d325c5f83","added_by":"auto","created_at":"2026-02-03 16:06:19","extension":"jpg","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":97435,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eEffect on centralization of different\u003c/strong\u003e \u003cstrong\u003enormalization methods for 21 tumor samples profiled using panel 2.\u003c/strong\u003e \u0026nbsp;The distribution of log2FC values generated by DESeq2 after read counting by featureCounts and normalization of read counts by 5 methods. a. Test set included 10 biopsied-confirmed carcinomas (red dots) and 11 soft tissue tumors (light blue dots) profiled by panel 2. b. Effects for a distinct set of low-grade and high-grade soft tissue tumors, graded according to the FNCLCC (French Federation of Cancer Centers Sarcoma Group) grading system [13] were compared with DESeq2 as in (a). Method abbreviations as in text.\u003c/p\u003e","description":"","filename":"2.jpg","url":"https://assets-eu.researchsquare.com/files/rs-8695099/v1/7f2ba54580370a89ec56617d.jpg"},{"id":101790561,"identity":"f1e01ac2-1d66-4536-a2b8-2994dc13a6e1","added_by":"auto","created_at":"2026-02-03 16:06:19","extension":"jpg","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":99554,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eEffect on PCA clustering of different\u003c/strong\u003e \u003cstrong\u003enormalization methods comparing soft tissue tumors and carcinomas. \u003c/strong\u003eDGE and PCA were performed as in Methods for 10 carcinomas (red dots) and 11 soft tissue tumors (light blue) profiled by panel 2, with read counting by featureCounts and normalization of reads by the indicated method. The same 21 samples were used as in Fig. 2a. Method abbreviations as in text.\u003c/p\u003e","description":"","filename":"3.jpg","url":"https://assets-eu.researchsquare.com/files/rs-8695099/v1/39e6234150260f168c7a8389.jpg"},{"id":101790567,"identity":"b6fb93a5-e68c-4197-b830-073d33c582e6","added_by":"auto","created_at":"2026-02-03 16:06:19","extension":"jpg","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":55810,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eComparison of\u003c/strong\u003e \u003cstrong\u003enormalization methods on outlier gene identification. \u003c/strong\u003eDGE was performed as in Methods, after read counting by featureCounts and normalization of reads by the indicated method, including identified 5 most stable genes across all samples type in the panel (IHK5), top most stable (Top10S) and using all genes. Venn diagrams illustrated the number of shared and significantly expressed genes in each category from these three normalization methods. Tumors analyzed for carcinoma versus soft tissue tumor comparison (a) and low-grade versus high-grade soft tissue tumors (b), as in Table S1.\u003c/p\u003e","description":"","filename":"4.jpg","url":"https://assets-eu.researchsquare.com/files/rs-8695099/v1/b56f1f9ca06e556d911c082d.jpg"},{"id":101880674,"identity":"c62041d3-af41-4dc9-b61c-a2802ccaaebc","added_by":"auto","created_at":"2026-02-04 15:05:06","extension":"jpg","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":98549,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eValidation of optimized the targeted RNAseq pipeline by comparison with total RNAseq data. \u003c/strong\u003eDESeq2 was performed using the targeted panel 2 (T) and total RNAseq was performed on 32 diagnostically challenging spindled cell tumors (19 carcinomas, 13 sarcomas) and compared with total RNAseq performed on the same RNA. a.\u003cstrong\u003e \u003c/strong\u003ePCA pattern shows a similar separation of carcinomas and sarcomas with targeted panel and total RNAseq. b. Ingenuity pathway analysis (IPA) shows similar canonical pathways that are upregulated (red) or downregulated (blue) in the carcinoma as compared to sarcomas. (c) Disease states similarity identified by IPA among the carcinomas (red) and sarcomas (blue) however are largely distinct between targeted panel and total RNAseq.\u003c/p\u003e","description":"","filename":"5.jpg","url":"https://assets-eu.researchsquare.com/files/rs-8695099/v1/672a318335ba9f9278a31a78.jpg"},{"id":102019069,"identity":"0df936a0-6c45-497b-8e50-b1e68b532d81","added_by":"auto","created_at":"2026-02-06 08:11:08","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":1546127,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-8695099/v1/ebc7e191-4c57-42f7-b7be-35b6509bde41.pdf"},{"id":101880908,"identity":"a9b53b0c-c8c0-4271-8136-18ec698a89f1","added_by":"auto","created_at":"2026-02-04 15:07:44","extension":"docx","order_by":0,"title":"","display":"","copyAsset":false,"role":"supplement","size":58105,"visible":true,"origin":"","legend":"","description":"","filename":"PanetalSupplementaryTables012525.docx","url":"https://assets-eu.researchsquare.com/files/rs-8695099/v1/6bdc505f01ce966bba87a708.docx"}],"financialInterests":"No competing interests reported.","formattedTitle":"Optimizing bioinformatic workflows to extract clinically usable gene expression data from targeted RNA sequencing panels: comparison with total RNAseq","fulltext":[{"header":"Introduction","content":"\u003cp\u003eRNA sequencing (RNAseq) by next-generation sequencing (NGS) is commonly used to detect diagnostic or therapy-related gene fusions in formalin-fixed paraffin-embedded (FFPE) tumor samples. To achieve rapid signout and maximal sensitivity through deep coverage, clinical laboratories typically perform fusion detection using limited/targeted gene sets of 50\u0026ndash;500 genes rather than assessing the entire transcriptome. Methods for targeted enrichment include bait-probe hybridization or anchored PCR amplification [\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e, \u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e, \u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e]. However, since relevant gene fusions are present in only a minority of tumors, many targeted RNAseq studies do not advance the diagnostic workup. Combining gene fusion detection with differential gene expression (DGE) analysis can increase the value of such panels. Clinically relevant goals of DGE in tumors include distinguishing site of origin (eg. bladder versus lung), distinguishing lower grade from higher grade soft tissue tumors and particularly the classification of poorly differentiated tumors where the differential includes carcinoma, melanoma and sarcoma.\u003c/p\u003e \u003cp\u003eThe first step in RNAseq-based DGE analysis is read counting where the number of sequence reads that align to each gene segment are quantified. With bait-probe capture and bridge amplification (Illumina) sequencing, NGS read counts have been shown to correlate with expression as determined by microarray or real-time PCR analysis [\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e, \u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e]. Although bioinformatics tools for extracting valid tumor gene expression from whole transcriptome data have been well-established [\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e], methods optimized for small/targeted RNAseq data sets are currently limited. Furthermore, the accuracy and reproducibility of any given read counting method needs to be established for a fully validated DGE assay to be used clinically [\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e].\u003c/p\u003e \u003cp\u003eGiven the wide variation in the quality of FFPE samples, fixation artifacts and tumor dilutional effects on differential gene expression can be significant, especially as compared to fresh/frozen tumor tissues [\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e]. Normalization or centralization methods for RNAseq DGE are usually applied to reduce bias due to variable tumor content and/or RNA degradation. Many normalization methods are available and have been evaluated for whole transcriptome RNAseq data [\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e, \u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e, \u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e, \u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e]. Whether these studies are valid or applicable to limited gene RNAseq panels has not been well-elucidated.\u003c/p\u003e \u003cp\u003eFor ease of clinical signout, data presentation using diagnostically meaningful outputs that accurately represent differences between tumors are critical for clinical utility. For example, heatmaps display comparisons of the entre sample and are most useful for classification discovery with static datasets but are difficult to apply dynamically. Principal component analysis (PCA), a linear dimensionality reduction technique that preserves large pairwise distances and variance in the data, may more easily highlight clustering of individual samples without explicitly specifying distinguishing gene sets. Similarly, t-distributed stochastic neighbor embedding (t-SNE) is another data reduction technique that can potentially easily highlight similarity of a new sample to the reference set of different tumors.\u003c/p\u003e \u003cp\u003eIn this study, we employed clinically validated\u0026thinsp;~\u0026thinsp;200-gene RNAseq panels optimized for fusion detection in solid tumors to systemically consider optimal methods of read counting, normalization and data visualization. We compared the outputs of the optimized pipelines for a targeted panel that included a limited number of genes included cell lineage and tumor grade assessment with parallel full RNAseq and showed highly comparable results in a set of poorly differentiated tumors.\u003c/p\u003e"},{"header":"Results","content":"\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e \u003ch2\u003eOptimizing read counting\u003c/h2\u003e \u003cp\u003eTo convert NGS output from the targeted RNAseq panels into gene expression data, raw sequence reads obtained from each gene segment were counted for each sample using 6 different tools: Samtools-view (referent clinical assay method), featureCounts, HTSeq-count, coverageBed, cuffdiff and countOverlaps. A sample output from a 16-sample RNAseq fusion assay is shown comparing usability, time-to-output and resource use. Using 8 CPU processors, featureCounts, took only 2 minutes to complete as compared to coverageBed at 9 minutes (\u003cb\u003eTable S3\u003c/b\u003e). In contrast, HTSeq-count, CuffDiff and Samtools-view produced outputs in 3, 4 and 6 hours, respectively.\u003c/p\u003e \u003cp\u003eThe total reads in the outputs of Samtools-view, featureCounts, HTSeq-count and coverageBed were similar for 10 samples (Table\u0026nbsp;\u003cspan refid=\"Tab1\" class=\"InternalRef\"\u003e1\u003c/span\u003e), whereas the read counts from cuffdiff outputs were much smaller due to larger genes producing high read counts being skipped in the calculation by the software. Excluding cuffdiff, the CV values for each of the 222 expressed genes were computed for each program. The median CV ranks of the 222 expressed genes from the four programs were also nearly identical (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e1\u003c/span\u003ea \u003cb\u003eand\u003c/b\u003e Table\u0026nbsp;\u003cspan refid=\"Tab2\" class=\"InternalRef\"\u003e2\u003c/span\u003e). The correlation of each expressed gene in the panel was compared with the validated Samtools-view value using Pearson correlation coefficients and showed high correlations (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e1\u003c/span\u003eb); the median rank of the correlation coefficients showed only small differences (Table\u0026nbsp;\u003cspan refid=\"Tab2\" class=\"InternalRef\"\u003e2\u003c/span\u003e). Considering the performance metrics of output time, precision and accuracy, featureCounts was rated as the best read counting method for this application.\u003c/p\u003e \u003cp\u003e \u003cdiv class=\"gridtable\"\u003e\u003ctable float=\"Yes\" id=\"Tab1\" border=\"1\"\u003e \u003ccaption language=\"En\"\u003e \u003cdiv class=\"CaptionNumber\"\u003eTable 1\u003c/div\u003e \u003cdiv class=\"CaptionContent\"\u003e \u003cp\u003eNumber of read counts in 10-sample tun using different read counting methods (assay 2).\u003c/p\u003e \u003c/div\u003e \u003c/caption\u003e \u003ccolgroup cols=\"11\"\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c1\" colnum=\"1\"\u003e\u003c/div\u003e \u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c2\" colnum=\"2\"\u003e\u003c/div\u003e \u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c3\" colnum=\"3\"\u003e\u003c/div\u003e \u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c4\" colnum=\"4\"\u003e\u003c/div\u003e \u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c5\" colnum=\"5\"\u003e\u003c/div\u003e \u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c6\" colnum=\"6\"\u003e\u003c/div\u003e \u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c7\" colnum=\"7\"\u003e\u003c/div\u003e \u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c8\" colnum=\"8\"\u003e\u003c/div\u003e \u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c9\" colnum=\"9\"\u003e\u003c/div\u003e \u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c10\" colnum=\"10\"\u003e\u003c/div\u003e \u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c11\" colnum=\"11\"\u003e\u003c/div\u003e \u003cthead\u003e \u003ctr\u003e \u003cth align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cem\u003eMethod\u003c/em\u003e\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c2\"\u003e \u003cp\u003eS-107\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c3\"\u003e \u003cp\u003eS-143\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c4\"\u003e \u003cp\u003eS-153\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c5\"\u003e \u003cp\u003eS-175\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c6\"\u003e \u003cp\u003eS-359\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c7\"\u003e \u003cp\u003eS-394\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c8\"\u003e \u003cp\u003eS-458\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c9\"\u003e \u003cp\u003eS-504\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c10\"\u003e \u003cp\u003eS-525\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c11\"\u003e \u003cp\u003eS-614\u003c/p\u003e \u003c/th\u003e \u003c/tr\u003e \u003c/thead\u003e \u003ctbody\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eSamtools-view\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c2\"\u003e \u003cp\u003e23863915\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e \u003cp\u003e84929701\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e38449531\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e103286313\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c6\"\u003e \u003cp\u003e61451551\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c7\"\u003e \u003cp\u003e52339994\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c8\"\u003e \u003cp\u003e173011118\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c9\"\u003e \u003cp\u003e32439398\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c10\"\u003e \u003cp\u003e60600048\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c11\"\u003e \u003cp\u003e72260163\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eFeatureCounts\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c2\"\u003e \u003cp\u003e23559590\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e \u003cp\u003e84701209\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e37895291\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e102612244\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c6\"\u003e \u003cp\u003e61105419\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c7\"\u003e \u003cp\u003e51932562\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c8\"\u003e \u003cp\u003e171937331\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c9\"\u003e \u003cp\u003e32037416\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c10\"\u003e \u003cp\u003e60076885\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c11\"\u003e \u003cp\u003e71446450\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eHTSeq-count\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c2\"\u003e \u003cp\u003e23505495\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e \u003cp\u003e84600839\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e37801803\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e102395919\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c6\"\u003e \u003cp\u003e60991722\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c7\"\u003e \u003cp\u003e51824717\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c8\"\u003e \u003cp\u003e171667472\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c9\"\u003e \u003cp\u003e31965039\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c10\"\u003e \u003cp\u003e59889358\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c11\"\u003e \u003cp\u003e71322334\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eCoverageBED\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c2\"\u003e \u003cp\u003e24306691\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e \u003cp\u003e85794224\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e39043406\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e104359380\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c6\"\u003e \u003cp\u003e62160975\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c7\"\u003e \u003cp\u003e53154627\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c8\"\u003e \u003cp\u003e174990470\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c9\"\u003e \u003cp\u003e32974310\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c10\"\u003e \u003cp\u003e61740756\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c11\"\u003e \u003cp\u003e73204253\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eCuffDiff\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c2\"\u003e \u003cp\u003e3863910\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e \u003cp\u003e3428790\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e62443716\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e7726710\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c6\"\u003e \u003cp\u003e4301250\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c7\"\u003e \u003cp\u003e5106000\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c8\"\u003e \u003cp\u003e9124866\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c9\"\u003e \u003cp\u003e4543470\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c10\"\u003e \u003cp\u003e9254736\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c11\"\u003e \u003cp\u003e10259952\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003c/tbody\u003e \u003c/colgroup\u003e \u003c/table\u003e\u003c/div\u003e \u003c/p\u003e \u003cp\u003e \u003cdiv class=\"gridtable\"\u003e\u003ctable float=\"Yes\" id=\"Tab2\" border=\"1\"\u003e \u003ccaption language=\"En\"\u003e \u003cdiv class=\"CaptionNumber\"\u003eTable 2\u003c/div\u003e \u003cdiv class=\"CaptionContent\"\u003e \u003cp\u003ePrecision and accuracy performance of different read counting methods.\u003c/p\u003e \u003c/div\u003e \u003c/caption\u003e \u003ccolgroup cols=\"5\"\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c1\" colnum=\"1\"\u003e\u003c/div\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c2\" colnum=\"2\"\u003e\u003c/div\u003e \u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c3\" colnum=\"3\"\u003e\u003c/div\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c4\" colnum=\"4\"\u003e\u003c/div\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c5\" colnum=\"5\"\u003e\u003c/div\u003e \u003cthead\u003e \u003ctr\u003e \u003cth align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cem\u003eRank\u003c/em\u003e\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c2\"\u003e \u003cp\u003e\u003cem\u003eCounting method\u003c/em\u003e\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c3\"\u003e \u003cp\u003e\u003cem\u003eMedian precision\u003c/em\u003e\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c4\"\u003e \u003cp\u003e\u003cem\u003eMedian accuracy\u003c/em\u003e\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c5\"\u003e \u003cp\u003e\u003cem\u003eOverall\u003c/em\u003e\u003c/p\u003e \u003c/th\u003e \u003c/tr\u003e \u003c/thead\u003e \u003ctbody\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e1\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eSamtools-view\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e \u003cp\u003e110\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003ena\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c5\"\u003e \u003cp\u003ena\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e2\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003ecoverageBed\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e \u003cp\u003e110\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e93\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c5\"\u003e \u003cp\u003e203\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e3\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003efeatureCounts\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e \u003cp\u003e110\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e111\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c5\"\u003e \u003cp\u003e221\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e4\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eHTSeq\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e \u003cp\u003e110\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e113\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c5\"\u003e \u003cp\u003e223\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003c/tbody\u003e \u003c/colgroup\u003e \u003ctfoot\u003e \u003ctr\u003e\u003ctd colspan=\"5\"\u003eNA: Not applicable as Samtools-view is the referent method.\u003c/td\u003e\u003c/tr\u003e \u003c/tfoot\u003e \u003c/table\u003e\u003c/div\u003e \u003c/p\u003e \u003c/div\u003e\n\u003ch3\u003eEffect of expression normalization strategies on tumor grouping and data centralization\u003c/h3\u003e\n\u003cp\u003eClinical biopsies/resections show a range of tumor cellularity and varying admixture of neoplastic and non-neoplastic cell types. Therefore, minimizing these effects on overall gene expression distribution through read count normalization is critical in highlighting diagnostically relevant sample differences. Especially in suboptimal FFPE samples, different normalization strategies are expected to have large effects. Using the raw read counts obtained with featureCounts as input (\u0026ldquo;NoNorm\u0026rdquo;), we compared the impact of normalization with \u003cem\u003eACTB\u003c/em\u003e, TraditionHK5, IdentifiedHK5, Top10H, Top10S and AllGenes (see Methods for definitions of these groups) followed by gene expression analysis in DESeq2.\u003c/p\u003e \u003cp\u003eOne readout of effective normalization assessed was a reduction in intragroup variation across tumors of similar types. This was assessed by calculating the average CV percentage values of read counts before or after normalization for groups of carcinomas, low-grade soft tissue tumors and high-grade sarcomas. The lowest CV values for each of the 3 tumor groups were achieved with the IdentifiedHK5, Top10S and AllGenes normalization methods (Table\u0026nbsp;\u003cspan refid=\"Tab3\" class=\"InternalRef\"\u003e3\u003c/span\u003e). The intragroup variation for these methods was significantly reduced compared to the NoNorm control condition (P\u0026thinsp;\u0026lt;\u0026thinsp;0.01, Wilcoxon signed rank test). Intra-tumor group variation was higher than NoNorm for the TraditionHK5 and \u003cem\u003eACTB\u003c/em\u003e normalization methods, with the latter having significantly higher average CV values for both the carcinoma and sarcoma groups.\u003c/p\u003e \u003cp\u003e \u003cdiv class=\"gridtable\"\u003e\u003ctable float=\"Yes\" id=\"Tab3\" border=\"1\"\u003e \u003ccaption language=\"En\"\u003e \u003cdiv class=\"CaptionNumber\"\u003eTable 3\u003c/div\u003e \u003cdiv class=\"CaptionContent\"\u003e \u003cp\u003eAverage CV percentage values after normalization by different methods using DESeq2.\u003c/p\u003e \u003c/div\u003e \u003c/caption\u003e \u003ccolgroup cols=\"8\"\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c1\" colnum=\"1\"\u003e\u003c/div\u003e \u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c2\" colnum=\"2\"\u003e\u003c/div\u003e \u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c3\" colnum=\"3\"\u003e\u003c/div\u003e \u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c4\" colnum=\"4\"\u003e\u003c/div\u003e \u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c5\" colnum=\"5\"\u003e\u003c/div\u003e \u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c6\" colnum=\"6\"\u003e\u003c/div\u003e \u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c7\" colnum=\"7\"\u003e\u003c/div\u003e \u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c8\" colnum=\"8\"\u003e\u003c/div\u003e \u003cthead\u003e \u003ctr\u003e \u003cth align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cem\u003eGroup\u003c/em\u003e\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c2\"\u003e \u003cp\u003e\u003cem\u003eNoNorm\u003c/em\u003e\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c3\"\u003e \u003cp\u003e\u003cem\u003eACTB\u003c/em\u003e\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c4\"\u003e \u003cp\u003e\u003cem\u003eTraditionHK5\u003c/em\u003e\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c5\"\u003e \u003cp\u003e\u003cem\u003eIdentifiedHK5\u003c/em\u003e\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c6\"\u003e \u003cp\u003e\u003cem\u003eTop10H\u003c/em\u003e\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c7\"\u003e \u003cp\u003e\u003cem\u003eTop10S\u003c/em\u003e\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c8\"\u003e \u003cp\u003e\u003cem\u003eAllGenes\u003c/em\u003e\u003c/p\u003e \u003c/th\u003e \u003c/tr\u003e \u003c/thead\u003e \u003ctbody\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eCarcinomas\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c2\"\u003e \u003cp\u003e93.17\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e \u003cp\u003e93.99\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e88.45*\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e84.99**\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c6\"\u003e \u003cp\u003e104.65**\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c7\"\u003e \u003cp\u003e86.67**\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c8\"\u003e \u003cp\u003e85.40**\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eSoft Tissue Tumors\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c2\"\u003e \u003cp\u003e82.93\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e \u003cp\u003e85.30*\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e84.20\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e83.35\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c6\"\u003e \u003cp\u003e89.09**\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c7\"\u003e \u003cp\u003e82.57\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c8\"\u003e \u003cp\u003e82.16\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eLow-grade Sarcoma\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c2\"\u003e \u003cp\u003e100.15\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e \u003cp\u003e96.67\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e84.36**\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e82.53**\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c6\"\u003e \u003cp\u003e101.60\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c7\"\u003e \u003cp\u003e83.33**\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c8\"\u003e \u003cp\u003e81.54**\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eHigh-grade Sarcoma\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c2\"\u003e \u003cp\u003e119.15\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e \u003cp\u003e127.86**\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e121.65\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e105.67**\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c6\"\u003e \u003cp\u003e141.57**\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c7\"\u003e \u003cp\u003e105.44**\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c8\"\u003e \u003cp\u003e105.76**\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003c/tbody\u003e \u003c/colgroup\u003e \u003ctfoot\u003e \u003ctr\u003e\u003ctd colspan=\"8\"\u003eNote: Comparisons \u003cem\u003evia\u003c/em\u003e \u0026ldquo;estimateSizeFactors\u0026rdquo; function in DESeq2 across different number of genes as control. ** and * represent 0.01 and 0.05 significantly different from the value without normalization NoNorm; red and black stars indicate significantly higher and lower, respectively.\u003c/td\u003e\u003c/tr\u003e \u003c/tfoot\u003e \u003c/table\u003e\u003c/div\u003e \u003c/p\u003e \u003cp\u003eUsing DESeq2 to compare gene expression in distinct tumor groups, the distributions of all the logFC values were calculated for downregulated and upregulated genes in carcinoma as compared to soft tissue tumors in assay 2 (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e2\u003c/span\u003ea). The distribution of the logFC values of upregulated genes was highly biased without normalization likely influenced by a small set of highly expressed set of genes related to extracellular matrix production such as collagen isoforms. Normalization using \u003cem\u003eACTB\u003c/em\u003e or multiple HK genes (TraditionHK5) exaggerated this effect. Normalization using Top10H better highlighted differences among down-regulated genes. Normalization using AllGenes and Top10S moved the distribution of the logFC values toward to 0 for upregulated genes, whereas the IdentifiedHK5 method centralized the distribution of the logFC values for both upregulated and downregulated genes.\u003c/p\u003e \u003cp\u003eTo assess the effects of normalization strategies on determining tumor grade by gene expression patterns, the distribution of the logFC values of the read counts was compared to a set of low-grade soft tissue tumors and high-grade sarcomas in assay 1 (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e2\u003c/span\u003eb). Without normalization, upregulated genes biased the distribution as did including all genes, whereas \u003cem\u003eACTB\u003c/em\u003e or the TraditionHK5 and Top10H methods improved the separation between the groups. IdentifiedHK5 and Top10S produced centralization of the distribution. Possibly due to effects on minimizing differences simply due to tumor cellularity, normalization using IdentifiedHK5 and Top10S produced excellent results for this assay indication.\u003c/p\u003e\n\u003ch3\u003eEffect of expression normalization strategies on tumor type clustering\u003c/h3\u003e\n\u003cp\u003eAn important rationale for normalization strategies in targeted tumor RNAseq assays is to enhance separation of tumor type clusters to improve mapping of samples of unknown lineage. Here, we evaluated several clustering and visualization methods, including PCA, t-SNE and Heatmap-clustering algorithms for their ease of interpretation by several molecular pathologists among the authors.\u003c/p\u003e \u003cp\u003eUsing a comparison of carcinoma and soft tissue tumor groups, the effects of read counts normalized by different methods on clustering by PCA was assessed. The IdentifiedHK5 and Top10S methods enhanced the separation of PCA clusters with IdentifiedHK5 showing best results (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e3\u003c/span\u003e). The other normalization methods reduced the separation of PCA clusters as compared to no normalization. Although a heatmap of the same data did highlight the tumor groups, the complicated display did not easily facilitate expression pattern of new samples during a quick review (not shown). With the IdentifiedHK5 normalization method, there was no obvious clustering with t-SNE for the carcinoma versus soft tissue diagnostic indication (not shown). Similar results are seen for the low-grade and high-grade tumor comparison (not shown).\u003c/p\u003e \u003cp\u003eUsing DESeq2, the most differentially expressed genes comparing carcinoma and soft tissue tumors, the IdentifiedHK5 method again showed the best results with a few more distinguishing genes (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e5\u003c/span\u003ea). A similar pattern was seen with the low-grade versus high-grade tumor sets as a separate distinct diagnostic indication (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e5\u003c/span\u003eb).\u003c/p\u003e\n\u003ch3\u003eComparison of targeted RNAseq outputs with total RNAseq\u003c/h3\u003e\n\u003cp\u003eTo validate the DGE output of the targeted panel 2, we compared the outputs with those obtained by total RNAseq on a diagnostically relevant set of 32 poorly differentiated tumors (19 carcinoma, 13 sarcoma). Although the number of expressed genes across all samples was vast differently (226 versus 20218; Table\u0026nbsp;\u003cspan refid=\"Tab4\" class=\"InternalRef\"\u003e4\u003c/span\u003e), the significantly differentially expressed genes between carcinoma and sarcoma cases (at adjusted p\u0026thinsp;\u0026lt;\u0026thinsp;0.05 level) were very similar for the shared genes (\u003cb\u003eSupplementary Table\u0026nbsp;6\u003c/b\u003e). The mapping of these cases into carcinoma and sarcoma clusters was also similar by PCA (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e5\u003c/span\u003ea). Similarly, the top differentially upregulated or downregulated canonical pathways by IPA analysis were also similar (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e5\u003c/span\u003eb). However, most of the disease states differentially regulated in sarcomatoid carcinoma and sarcoma, as identified by IPA, were distinct when exome data was compared to the targeted panel (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e5\u003c/span\u003ec). Although S100 family signaling was a shared upregulated pathway in sarcomas by both analyses, the targeted panel identified more specific signaling pathways in the carcinoma group.\u003c/p\u003e \u003cp\u003e \u003cdiv class=\"gridtable\"\u003e\u003ctable float=\"Yes\" id=\"Tab4\" border=\"1\"\u003e \u003ccaption language=\"En\"\u003e \u003cdiv class=\"CaptionNumber\"\u003eTable 4\u003c/div\u003e \u003cdiv class=\"CaptionContent\"\u003e \u003cp\u003eComparison of expressed genes in optimized outputs for targeted RNAseq (panel 2) and full RNAseq.\u003c/p\u003e \u003c/div\u003e \u003c/caption\u003e \u003ccolgroup cols=\"6\"\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c1\" colnum=\"1\"\u003e\u003c/div\u003e \u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c2\" colnum=\"2\"\u003e\u003c/div\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c3\" colnum=\"3\"\u003e\u003c/div\u003e \u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c4\" colnum=\"4\"\u003e\u003c/div\u003e \u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c5\" colnum=\"5\"\u003e\u003c/div\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c6\" colnum=\"6\"\u003e\u003c/div\u003e \u003cthead\u003e \u003ctr\u003e \u003cth align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cem\u003eDatasets\u003c/em\u003e\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c2\"\u003e \u003cp\u003e\u003cem\u003eTotal genes\u003c/em\u003e\u003c/p\u003e \u003cp\u003e\u003cem\u003estudied\u003c/em\u003e\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c3\"\u003e \u003cp\u003e\u003cem\u003eCorrelation coefficient for shared genes*\u003c/em\u003e\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c4\"\u003e \u003cp\u003e\u003cem\u003eExpressed\u003c/em\u003e\u003c/p\u003e \u003cp\u003e\u003cem\u003egenes\u003c/em\u003e\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c5\"\u003e \u003cp\u003e\u003cem\u003eExpressed genes (%)\u003c/em\u003e\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c6\"\u003e \u003cp\u003e\u003cem\u003eExpressed genes shared\u003c/em\u003e\u003c/p\u003e \u003c/th\u003e \u003c/tr\u003e \u003c/thead\u003e \u003ctbody\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eTotal RNAseq\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c2\"\u003e \u003cp\u003e20218\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c3\" morerows=\"1\" rowspan=\"2\"\u003e \u003cp\u003e0.95\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e2591\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e12.8\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c6\" morerows=\"1\" rowspan=\"2\"\u003e \u003cp\u003e32\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eTargeted RNAseq\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c2\"\u003e \u003cp\u003e226\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e47\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e14.2\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003c/tbody\u003e \u003c/colgroup\u003e \u003ctfoot\u003e \u003ctr\u003e\u003ctd colspan=\"6\"\u003eNote: Expressed genes selected by cutoffs of adjust p-value 0.05 and log2fc \u0026lt;= -0.7 \u0026amp; log2fc\u0026thinsp;\u0026gt;\u0026thinsp;=\u0026thinsp;0.7.\u003c/td\u003e\u003c/tr\u003e \u003c/tfoot\u003e \u003c/table\u003e\u003c/div\u003e \u003c/p\u003e"},{"header":"Discussion","content":"\u003cp\u003eUsing clinical-grade\u0026thinsp;~\u0026thinsp;200-gene RNAseq assays originally developed for gene fusion detection, we have systematically evaluated their suitability for tumor lineage assessment and grading. The best methods for read counting (for pipeline optimization) and data normalization using DGE and PCA for evaluation, were determined. The output of the optimized targeted panel pipeline was then validated for a difficult application (distinguishing sarcomatoid poorly differentiated tumors representing either carcinoma or sarcoma) using total RNAseq.\u0026nbsp;These comparisons provide a model for validating these methods for routine clinical use and highlight the importance of closely matching the methods to the design and goals of each assay.\u003c/p\u003e \u003cp\u003eAs highlighted by several use cases, the purposes for performing differential gene expression by RNAseq in tumors are varied. A common goal is to highlight one or more highly overexpressed lineage-associated genes in any given sample that may assist diagnosis. These uncommon expression patterns often signal cellular differentiation characteristics of a specific tumor type (eg high \u003cem\u003eMYOD1\u003c/em\u003e indicating skeletal muscle differentiation in rhabdomyosarcoma) [12{]. Another goal is to superimpose the overall gene expression pattern of any given tumor against a model set of different classes of tumors (e.g. carcinoma, melanoma and soft tissue tumors) to help determine site of origin when other biomarkers and microscopic studies are not informative. However, given the clinical impact of both indications, RNAseq analysis methods need to be carefully validated, and outputs need to be relatively easy to interpret. Targeted RNAseq often includes targets whose expression can vary widely based on specific molecular aberrations, which can cause some widely used software tools or methods to perform inaccurately. We therefore considered each step in the pipeline with regards to the clinical use cases.\u003c/p\u003e \u003cp\u003eMany read counting algorithms have been developed and used for whole transcriptome RNAseq.\u0026nbsp;As shown in Table S3, different read counting methods vary in speed performance. Our prior method, Samtools built-in program \u0026ldquo;view\u0026rdquo; calculates sequence reads in the targeted regions from alignment BAMs accurately when validated with the proper option settings [\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e]. But the long processing times can cause delays in clinical reporting. CuffDiff is another widely used software tool used for DGE analysis in whole transcriptome RNAseq [\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e]. It has a function to generate the number of sequence reads per gene in a sample. HTSeq-count is another widely used tool [\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e] BEDTools/coverageBed with the \u0026ldquo;-count\u0026rdquo; option is another option that allows multiple samples to be analyzed in parallel [\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e]. Corchete, et al [\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e] compared six methods of read counting for RNAseq data and reported HTSeq-count was superior. However, their study did not include coverageBed, Samtools-view or featureCounts which is optimized for efficient chromosome hashing and feature blocking techniques and parallel analysis. [\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e]. They also did not perform assessments of speed performance which is critical for clinical applications.\u003c/p\u003e \u003cp\u003eComparing speeding, precision accuracy and ease of use, we found that featureCounts and coverageBed had a significant advantage over other methods in output time, with featureCounts being the fastest. By statistical measures, featureCounts, coverageBed, HTSeq-count and Samtools-view have highly similar outputs. However, we found that CuffDiff underestimates sequence reads significantly which was traced to the limitations in the maximum number of reads in buffer of the software allowance which produces trimming of some reads from larger genes. Given that featureCounts and HTSeq-count require conversion of BED into GTF file, the ability of coverageBed to directly consume a BAM file is a great feature for smaller panels where the BAM files are not very large.\u003c/p\u003e \u003cp\u003eWhen small gene panels are employed for DGE, minimizing skewing of the overall gene expression distribution (or centralization) is critical. To find the optimal normalization method of sequence reads for DGE from small RNAseq fusion panels, we compared multiple common employed strategies and a lab/assay-specific method accounting for the sample and extraction patterns in our laboratory. In the latter, we identified the five most consistently and stably expressed genes across the tumors typically analyzed by our RNAseq fusion panel, aka empirically determined mostly stable genes. We then compared normalization using these 5 genes to other normalization methods. This approach has impact on centralizing the distribution of log2FC values of DGE and improving or enhancing the separation of PCA clusters. In addition, this method could reduce intragroup variances significantly, as did the AllGenes and Top10S methods. The significantly expressed genes identified by DGE were also mostly similar with these three methods. Overall, for this assay/indication, normalization method with the empirically top 5 most stable genes was the superior normalization method.\u003c/p\u003e \u003cp\u003eThe traditional normalization methods using a single housekeeping gene such as ACTB and multiple highly expressed genes (e.g. Top10H) as references, respectively performed poorly in reducing intra-tumor group variation. Similarly, the default normalization method for DESeq2 (AllGenes) was not optimal for these targeted panels. These findings emphasize that assessment of different centralization/normalization methods, especially for more targeted RNAseq panel, should be included in the validation of each new assay with sample sets tuned for the specific clinical applications.\u003c/p\u003e \u003cp\u003eConsideration of different data display methods is also critical given the limited time available for analysis for clinical assays. When mapping expression patterns of new/unknown tumors to existing data sets, we found that PCA is the best sample clustering method. By displaying typical tumor groups clearly (eg carcinoma vs sarcoma, high-grade vs low-grade), it can efficiently identify when an unknown sample maps well within a group as opposed to being an outlier where the study is non-informative. Although heatmaps are a useful visualization method when tumor expression patterns are highly similar (eg cultured tumor cells exposed to different drug doses), it is most difficult to visualize tumor-group differences in a single new sample.\u003c/p\u003e \u003cp\u003eUsing a carcinoma versus sarcoma tumor set, an optimized targeted panel with some lineage-specific genes included (panel 2) showed equivalent separation by PCA with parallel data from total RNAseq.\u0026nbsp;In addition, many of the top tumor class discriminator were shared by the two assays. This similar result was further noted by IPA comparison of the top differentially regulated canonical signaling pathways (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e5\u003c/span\u003e). However, disease states mapping by IPA among the carcinomas (red) and sarcomas (blue) was largely distinct between targeted panel and total RNAseq.\u0026nbsp;This was likely due to the effects of genes not present in the targeted panel in providing subclass differentiation in small tumor sets.\u003c/p\u003e \u003cp\u003eIn summary, we identified featureCounts as the optimal read counting method, five assay-specific most stable genes as the best method for normalization prior to DESeq2 and PCA as the easiest to interpret sample clustering method for this targeted RNAseq panel clinical application.\u003c/p\u003e"},{"header":"Methods","content":" \u003cdiv id=\"Sec8\" class=\"Section2\"\u003e \u003cdiv id=\"Sec9\" class=\"Section3\"\u003e \u003ch2\u003eNGS library design and sample sets\u003c/h2\u003e \u003cp\u003eThe gene fusion NGS panels employed in this study included a custom 190-gene (panel 1) and 230-gene (panel 2) custom RNAseq panels used at The Ohio State University James Molecular Laboratory to detect diagnostically relevant gene fusions in human tumors. The 190-gene design included only a few typical housekeeping genes (\u003cem\u003eACTB, MYH9\u003c/em\u003e) with no gene content explicitly designed for DGE. Over 500 clinical cases, the frequency of reportable gene fusions by panel 1 varied by diagnosis approaching 20% for soft tissue tumors versus ~\u0026thinsp;5% for carcinomas. To improve utility for gene expression, panel 2 incorporated additional stably expressed/housekeeping and lineage-specific genes to aid in separation of poorly differentiated carcinomas and sarcoma. In 250 clinical cases using panel 2, 226 genes were routinely expressed at some level, with other genes only expressed when a particular fusion was present. The frequency of oncogenic fusion detection was similar to panel 1.\u003c/p\u003e \u003cp\u003eFor validation of the design of panel 2 using an optimized pipeline, full RNAseq (total transcriptome except for ribosomal genes) was performed for a set of diagnostically challenging spindled cell tumors. These represented 32 sarcomatoid poorly differentiated tumors that were diagnostically as either carcinoma or sarcoma following routine histopathology/immunohistochemistry workup. Final diagnosis was rendered by a soft tissue pathologist following comprehensive DNA mutation profiling and correlation with radiologic appearances and clinical presentation. All assays used and samples procured were obtained as part of routine clinical testing. The use of this data for development of the analytic pipeline was reviewed by the institutional review board with waiver of consent provided.\u003c/p\u003e \u003cp\u003eAll NGS protocols utilized total tumor RNA extracted from formalin-fixed paraffin-embedded tissue (FFPE) using PureLink FFPE RNA Isolation Kit (Invitrogen/ThermoFisher) and employed library preparation with DNA digestion and ribosomal RNA depletion (KAPA RNA HyperPrep Kit or Watchmaker Polaris). For the targeted panels, this was followed by hybridization with probes covering the full exonic regions of all target genes (xGen, IDT, Coralville, IA) designed through the GOAL consortium [\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e]. This size of the panel allowed 10\u0026ndash;16 samples to be run on the Illumina NovaSeq SP flow cell, with adequate depth of coverage (~\u0026thinsp;20\u0026ndash;50,000 mean reads per sample per targeted area).\u003c/p\u003e \u003cp\u003eTo assess the suitability of methods across different gene sets, tumors sequenced with both panel 1 and panel 2 were included. For comparative analysis of performance of different read counting methods, 10 soft tissue tumor samples from panel 2 were used. For the comparison of normalization and clustering methods, 36 tumor samples from panel 2 were used, with 10 carcinomas and 11 soft tissue tumors as model comparators for DGE analysis (\u003cb\u003eTable \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e\u003c/b\u003e).\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e\n\u003ch3\u003eInitial clinical pipeline\u003c/h3\u003e\n\u003cp\u003eFor clinical reporting, we employed a custom pipeline for fusion detection (\u0026ldquo;FindRNAFusion\u0026rdquo;): paired-end FASTQ files were downloaded to a high-performance HPE Linux server (256 CPU processors, 755 GB RAM memory, and RHEL 9.0 OS) from a mounted Illumina BaseSpace instance. The raw sequence reads of each sample were mapped to Hg19 (Human Genome version 19) to generate alignment BAM files using STAR [\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e]. Each BAM file was indexed by SAMTOOLS. The BAM files were then used to make fusion calls using Arriba [\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e]. After filtering out artifacts and low-level fusion calls (\u0026lt;\u0026thinsp;10\u0026ndash;20 supporting reads), clinically reportable fusions were summarized in a text report file and a pdf file with graphical display. In this pipeline, the number of sequence reads/coverage in each targeted gene and in each targeted region/exon were also computed from the BAM files using each read counting method, as discussed below. These sequence reads were then normalized, using the methods described below. Subsequently, genes with very low reads were filtered out using a data-based threshold for maximum-based filters proposed by Rau et al [\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e] and the remaining genes with normalized read counts were used for DGE analysis using DESeq2 [\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e].\u003c/p\u003e \u003cdiv id=\"Sec11\" class=\"Section2\"\u003e \u003ch2\u003eRead count method comparisons\u003c/h2\u003e \u003cp\u003eFive read counting methods, Samtools-view [\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e], HTSeq-count [\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e], featureCounts [\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e], coverageBed [\u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e], and CuffLinks [\u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e23\u003c/span\u003e], were selected for comparative analysis. Samtools-view, which was our previously validated referent method, uses a BED file and a BAM file as input and outputs the number of reads in each targeted range (command line: \u0026ldquo;samtools view -F 0x04 -q 20 -c -@ 8\u0026rdquo;). This accuracy of this output was validated by comparing the count of mapped reads to manual inspection/calculation of read depth for a set of genes in the Integrative Genomics Viewer (IGV) [\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e] and proportion of fusion-positive tumor cells by in situ hybridization for selected common gene fusions. In contrast, featureCounts uses targeted genes in GTF file and BAM file as input (command line: \u0026ldquo;featureCounts -p -M -O -C -g gene_name --minOverlap 1 --maxMOp 30 -Q 20 \u0026ndash;T 8\u0026rdquo;). HTSeq-count uses the same files as featureCounts as input (command line: \u0026ldquo;htseq-count -i gene_name --max-reads-in-buffer 50000000 -s no\u0026rdquo;). CuffDiff was run with option \u0026ldquo;--total-hits-norm TRUE\u0026rdquo; with the same files as featureCounts as input. coverageBed was run with option \u0026ldquo;-count\u0026rdquo; and \u0026ldquo;-a BED file\u0026rdquo; and \u0026ldquo;-b BAM file\u0026rdquo; as input. For this study, a Perl script was written to run program commands in parallel with a batch of NGS samples.\u003c/p\u003e \u003cp\u003eConcordance of different methods was assessed among the 226 routinely expressed genes in panel 2 with the fully validated Samtools-view method as the referent. The precision of each method was assessed using the median rank from the coefficient of variation (CV) values of the 226 routinely expressed genes for each sample [\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e]. We also ranked CV values for each method. We also computed the median of the ranks of the 226 commonly expressed genes as the precision index of each method with lower values correlated with increased precision. The Pearson correlation coefficient (r) was computed to assess the association between the outputs of Samtools-view and the other methods. These correlation coefficients were calculated using the statistical functions in Microsoft Excel.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec12\" class=\"Section2\"\u003e \u003ch2\u003eNormalization comparisons\u003c/h2\u003e \u003cp\u003eDESeq2 is a widely used software tool for differential gene expression analysis and was chosen here as it performs well for larger RNAseq datasets (sample size\u0026thinsp;\u0026gt;\u0026thinsp;=\u0026thinsp;6) [\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e]. The standard normalization implemented in DESeq2 is relative log expression (RLE). The calculation of size factors is performed through the \u0026ldquo;estimateSizeFactor\u0026rdquo; function. The default for RLE normalization in DESeq2 is using all genes (\u0026ldquo;AllGenes\u0026rdquo;). We compared this result with those using a highly expressed \u0026ldquo;housekeeping\u0026rdquo; gene \u003cem\u003eACTB\u003c/em\u003e (HK), a panel of five traditionally used housekeeping genes, namely \u003cem\u003eACTB, MYH9, RANBP2, PRKACA\u003c/em\u003e and \u003cem\u003eTFG\u003c/em\u003e (TraditionHK5) and to the top ten highly expressed genes as determined by the average number of reads in all the samples (Top10H) and the top 10 genes with smallest CV values in all the samples in each dataset (Top10S). Five genes (\u003cem\u003eCREBBP, BRAF, BRD4, ATF1\u003c/em\u003e and \u003cem\u003eCREB1\u003c/em\u003e) that were recurrently identified as among the top ten most stable in the range of test samples for panels 1 and 2 were also assessed (IdentifiedHK5, \u003cb\u003eTable S2\u003c/b\u003e). These six normalization approaches were compared to no normalization for their effects on statistical and visualized centralization of the data. The CV% value representing the percentage of the standard deviation to the mean per gene across samples was computed using the formula CV% = STDEV.P * 100/AVERAGE.\u003c/p\u003e \u003cp\u003eTotal RNAseq data produced from the same RNA samples and library preparation as the targeted panels were analyzed using the core pipelines normalized using the 5 most stable genes from the targeted panel as normalization and subjected to DEQ2 as above. For pathway analysis, QIAGEN Ingenuity Pathway Analysis (IPA) was used (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://www.qiagenbioinformatics.com\u003c/span\u003e\u003cspan address=\"https://www.qiagenbioinformatics.com\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e).\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec13\" class=\"Section2\"\u003e \u003ch2\u003eClustering and visualization\u003c/h2\u003e \u003cp\u003eThe three clustering and data visualization methods employed are summarized in \u003cb\u003eTable S4\u003c/b\u003e. Principal component analysis (PCA), t-distributed stochastic neighbor embedding (t-SNE) and heatmap-clustering were selected for comparative analysis in sample clustering and gene expression. SRplot [\u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e26\u003c/span\u003e] was used to perform principal component analysis (PCA) and heatmap-clustering. t-SNE-Java (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://github.com/lejon/T-SNE-Java\u003c/span\u003e\u003cspan address=\"https://github.com/lejon/T-SNE-Java\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) was implemented to generate t-SNE graphics for clustering and visualization.\u003c/p\u003e \u003c/div\u003e"},{"header":"Declarations","content":"\u003ch2\u003eRelevant Disclosures/Conflicts of Interest\u003c/strong\u003e\u003c/h2\u003e\n\u003cp\u003eNone.\u003c/p\u003e\n\u003ch2\u003eFunding statement:\u003c/h2\u003e\n\u003cp\u003eNo relevant funding to declare.\u003c/p\u003e\n\u003ch2\u003eAuthor Contribution\u003c/h2\u003e\n\u003cp\u003eThe Authors participated in the bioinformatic analysis (X.P., M.H., D.J.), data analysis (R.S., N.M., D.C., D.J.), performed the experiments (R.S., Y.C., D.C.) or provided molecular pathology review, data integration and case selection (D.J, Y.C., Y.H., C.M., W.Z., M.A.). D.J. and X.P. wrote the main manuscript. All authors reviewed the manuscript and provided feedback.\u003c/p\u003e\n\u003ch2\u003eAcknowledgement\u003c/h2\u003e\n\u003cp\u003eThe Authors thank the medical technologists and the data scientists of the Polaris Molecular Laboratory for performing and analyzing the targeted RNAseq panels. The diagnostic work of soft tissue pathologists Hans Iwenofu and Swati Satturwar for the diagnostic challenging dataset is noted.\u003c/p\u003e\n\u003ch2\u003eData Availability\u003c/h2\u003e\n\u003cp\u003eThe data underlying this article will be shared on reasonable request to the corresponding author.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\u003cli\u003e\u003cspan\u003eCurion, F. et al. Targeted RNA sequencing enhances gene expression profiling of ultra-low input samples. \u003cem\u003eRNA Biol.\u003c/em\u003e \u003cb\u003e17\u003c/b\u003e, 1741\u0026ndash;1753 (2020).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHeyer, E. E. et al. Diagnosis of fusion genes using targeted RNA sequencing. \u003cem\u003eNat. Commun.\u003c/em\u003e \u003cb\u003e10\u003c/b\u003e, 1388 (2019).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eCapone, I. et al. Targeted RNA sequencing analysis for fusion transcript detection in tumour diagnostics: assessment of bioinformatic tools reliability in FFPE samples. \u003cem\u003eExplor. Target. Antitumor Ther.\u003c/em\u003e \u003cb\u003e3\u003c/b\u003e, 582\u0026ndash;597 (2022).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eCorchete, L. A. et al. Systematic comparison and assessment of RNA-seq procedures for gene expression quantitative analysis. \u003cem\u003eSci. Rep.\u003c/em\u003e \u003cb\u003e10\u003c/b\u003e, 19737 (2020).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ede Brito, M. W., de Carvalho, S. S., Mota, M. B. \u0026amp; Mesquita, R. D. RNA-seq validation: software for selection of reference and variable candidate genes for RT-qPCR. \u003cem\u003eBMC Genom.\u003c/em\u003e \u003cb\u003e25\u003c/b\u003e, 697 (2024).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eFu, C. et al. Targeted RNA-seq assay incorporating unique molecular identifiers for improved quantification of gene expression signatures and transcribed mutation fraction in fixed tumour samples. \u003cem\u003eBMC Cancer\u003c/em\u003e. \u003cb\u003e21\u003c/b\u003e, 114 (2021).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLi, J., Fu, C., Speed, T. P., Wang, W. \u0026amp; Symmans, W. F. Accurate RNA sequencing from formalin-fixed cancer tissue to represent high-quality transcriptome from frozen tissue. \u003cem\u003eJCO Precision Oncol.\u003c/em\u003e \u003cb\u003e2\u003c/b\u003e, 1\u0026ndash;9 (2018).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eV\u0026auml;likangas, T., Suomi, T. \u0026amp; Elo, L. L. A systematic evaluation of normalization methods in quantitative label-free proteomics. \u003cem\u003eBrief. Bioinform\u003c/em\u003e. \u003cb\u003e19\u003c/b\u003e, 1\u0026ndash;11 (2018).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eZhao, Y. et al. TPM, FPKM or normalized counts? A comparative study of quantification measures for the analysis of RNA-seq data from the NCI Patient-Derived Models Repository. \u003cem\u003eJ. Transl Med.\u003c/em\u003e \u003cb\u003e19\u003c/b\u003e, 269 (2021).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBushel, P. R. et al. Comparison of normalization methods for analysis of TempO-Seq targeted RNA sequencing data. \u003cem\u003eFront. Genet.\u003c/em\u003e \u003cb\u003e11\u003c/b\u003e, 594 (2020).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eEvans, C., Hardin, J. \u0026amp; Stoebel, D. M. Selecting between-sample RNA-seq normalization methods from the perspective of their assumptions. \u003cem\u003eBrief. Bioinform\u003c/em\u003e. \u003cb\u003e19\u003c/b\u003e, 776\u0026ndash;792 (2018).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eAvenarius, M. R. et al. Integrated molecular profiling of rhabdomyosarcoma subtypes by targeted RNA-seq. \u003cem\u003emedRxiv (\u003c/em\u003e (2024).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eAmin, M. B. et al. AJCC Cancer Staging Manual, 8th ednSpringer, New York,. 13 (2017).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eAisner, D. L. et al. The Genomics Organization for Academic Laboratories (GOAL): a vision for a genomics future for academic pathology. \u003cem\u003eAcad. Pathol.\u003c/em\u003e \u003cb\u003e10\u003c/b\u003e, 100090 (2023).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDobin, A. et al. STAR: ultrafast universal RNA-seq aligner. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cb\u003e29\u003c/b\u003e, 15\u0026ndash;21 (2013).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eUhrig, S. et al. Accurate and efficient detection of gene fusions from RNA sequencing data. \u003cem\u003eGenome Res.\u003c/em\u003e \u003cb\u003e31\u003c/b\u003e, 448\u0026ndash;460 (2021).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eRau, A., Gallopin, M., Celeux, G. \u0026amp; Jaffr\u0026eacute;zic, F. Data-based filtering for replicated high-throughput transcriptome sequencing experiments. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cb\u003e29\u003c/b\u003e, 2146\u0026ndash;2152 (2013).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLove, M. I., Huber, W. \u0026amp; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. \u003cem\u003eGenome Biol.\u003c/em\u003e \u003cb\u003e15\u003c/b\u003e, 550 (2014).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDanecek, P. et al. Twelve years of SAMtools and BCFtools. \u003cem\u003eGigaScience.\u003c/em\u003e 10, giab008 (2021).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eAnders, S., Pyl, P. T. \u0026amp; Huber, W. HTSeq\u0026mdash;a Python framework to work with high-throughput sequencing data. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cb\u003e31\u003c/b\u003e, 166\u0026ndash;169 (2015).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLiao, Y., Gordon, K., Smyth, G. K. \u0026amp; Shi, W. featureCounts: an efficient general-purpose program for assigning sequence reads to genomic features. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cb\u003e30\u003c/b\u003e, 923\u0026ndash;930 (2014).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eQuinlan, A. R. \u0026amp; Hall, I. M. BEDTools: a flexible suite of utilities for comparing genomic features. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cb\u003e26\u003c/b\u003e, 841\u0026ndash;842 (2010).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eTrapnell, C. et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. \u003cem\u003eNat. Protoc.\u003c/em\u003e \u003cb\u003e7\u003c/b\u003e, 562\u0026ndash;578 (2012).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eRobinson, J. T. et al. Integrative Genomics Viewer. \u003cem\u003eNat. Biotechnol.\u003c/em\u003e \u003cb\u003e29\u003c/b\u003e, 24\u0026ndash;26 (2011).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLi, D. et al. An evaluation of RNA-seq differential analysis methods. \u003cem\u003ePLoS One\u003c/em\u003e. \u003cb\u003e17\u003c/b\u003e, e0264246 (2022).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eTang, D. et al. SRplot: a free online platform for data visualization and graphing. \u003cem\u003ePLoS One\u003c/em\u003e. \u003cb\u003e18\u003c/b\u003e, e0294236 (2023).\u003c/span\u003e\u003c/li\u003e\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":true,"hideJournal":true,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":true,"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":"Read counting, normalization, clustering, tumor grading, cell lineage assessment","lastPublishedDoi":"10.21203/rs.3.rs-8695099/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-8695099/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eTargeted RNA sequencing (RNAseq) is widely used to detect gene fusions in tumors but clinical diagnostic use of expression data from these panels in fusion-negative cases has been limited. To facilitate this application, we evaluated methods for sequence read counting and gene normalization to optimize them for smaller gene sets. We present comparative methods to derive differential gene expression (DGE) data using\u0026thinsp;~\u0026thinsp;200-gene clinically validated RNAseq fusion panels and compared them to parallel full RNAseq.\u0026nbsp;We compared five methods for read counting, demonstrating that featureCounts is the most rapid and robust. For normalization prior to DGE with DESeq2, we compared five different normalization strategies and showed normalization using the 5 most stably expressed genes provided optimal centralization for these smaller gene sets. DGE output was assessed by principal component analysis (PCA), t-SNE and heatmap-clustering. The final pipeline was validated using PCA and pathway analysis by comparison with full RNAseq separately performed on a common set of challenging tumors with comparable results observed with the targeted gene panel. Overall, we show using an optimized bioinformatic pipeline that usable gene expression data can be obtained from smaller targeted RNAseq panels to maximize the clinical utility of these assays.\u003c/p\u003e","manuscriptTitle":"Optimizing bioinformatic workflows to extract clinically usable gene expression data from targeted RNA sequencing panels: comparison with total RNAseq","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2026-02-03 16:06:14","doi":"10.21203/rs.3.rs-8695099/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":"ab82aaff-8e39-40fc-8dea-bfcdc2e99db0","owner":[],"postedDate":"February 3rd, 2026","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[{"id":61872361,"name":"Biological sciences/Biological techniques"},{"id":61872362,"name":"Health sciences/Biomarkers"},{"id":61872363,"name":"Biological sciences/Cancer"},{"id":61872364,"name":"Biological sciences/Computational biology and bioinformatics"},{"id":61872365,"name":"Biological sciences/Genetics"},{"id":61872366,"name":"Health sciences/Oncology"}],"tags":[],"updatedAt":"2026-02-06T08:10:41+00:00","versionOfRecord":[],"versionCreatedAt":"2026-02-03 16:06:14","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-8695099","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-8695099","identity":"rs-8695099","version":["v1"]},"buildId":"XKTyCvWXoU3ODBz1xrDgd","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}