A Batch-Processing Pipeline for RNA-seq Gene Abundance Estimation

Background

Public RNA-seq repositories provide access to extensive datasets generated across species, tissues, developmental stages, and experimental conditions. The NCBI Sequence Read Archive (SRA) enables the reuse of raw sequencing data for transcriptomic analyses beyond the objectives of the original studies [1]. These resources facilitate comparative analyses, hypothesis generation, and data integration without requiring additional sequencing experiments [2]. However, datasets generated by independent studies often differ in sequencing platform, read length, library layout, sequencing depth, metadata quality, and analytical processing, hindering their standardized integration and comparative interpretation. Several computational strategies are available for RNA-seq processing. Alignment-based workflows commonly use splice-aware aligners such as STAR [3] or HISAT2 [4], whereas lightweight approaches rely on pseudoalignment or selective-alignment tools such as kallisto [5] and Salmon [6]. In this pipeline, STAR is a requirement of the architecture rather than an interchangeable preference, because RSEM quantifies directly from the transcriptome-coordinate BAM file that STAR produces via --quantMode TranscriptomeSAM; existing protocols guide users through quality assessment, read preprocessing, alignment, transcript abundance estimation, differential expression analysis, and downstream visualization [7,8]. More comprehensive solutions based on workflow managers, including Nextflow and nf-core/rnaseq, provide scalable and portable processing across computing environments [9,10]. Nevertheless, these approaches may require familiarity with workflow-management systems, substantial computational infrastructure, or advanced bioinformatics expertise. Furthermore, the selection and version of alignment and quantification tools can influence abundance estimates and downstream analyses, particularly for low-abundance genes and transcripts with complex isoform structures [11,12]. This protocol presents a modular batch-processing pipeline for RNA-seq abundance estimation from multiple sequencing runs. The workflow integrates data retrieval, quality assessment, read preprocessing, alignment, transcript abundance estimation, and matrix generation within a transparent and configurable framework. Gene- and isoform-level abundances are estimated using RSEM, which accounts for ambiguously mapped reads and has demonstrated competitive performance relative to alternative RNA-seq quantification methods [11,13]. The pipeline generates standardized TPM and FPKM gene abundance matrices together with individual RSEM abundance files, quality-control summaries, and aggregated MultiQC reports. Transcript-level estimates, reported as individual isoform-level (.isoforms.results) outputs, are deliberately preserved alongside the gene-level files to support downstream gene-level inference through tools such as tximport combined with DESeq2 or edgeR, following Soneson et al. [14], when differential expression analysis is the end goal. This pipeline is executed directly on Linux workstations or servers and does not require a dedicated workflow manager. It incorporates sample-level status tracking, automatic download retries, restartable execution, configurable computational resources, and controlled removal of intermediate files. These features enable reproducible and storage-efficient batch processing of multiple RNA-seq runs while avoiding unnecessary repetition of successfully completed analyses. The workflow can be applied to any organism for which compatible reference genomes, gene annotations, and RNA-seq metadata are available.