Comparing the normalization methods for the differential analysis of Illumina high-throughput RNA-Seq data

Transcriptome deep sequencing (RNA-seq) data analysis is often affected by reference bias introduced by the use of a single linear reference (SLR) genome. Graph-based pangenomes can mitigate this bias by integrating the SLR genome with complex genetic variations within a species; however, their application remains limited owing to a lack of dedicated analytical tools. Here, we present PanGraphRNA, an integrated bioinformatics platform for RNA-seq data analysis using a graph pangenome as reference. Built on the Galaxy web-based framework, PanGraphRNA provides functional modules for constructing, evaluating, and applying graph pangenomes across different population scales, thus enabling accessibility, traceability, and reproducibility throughout the analysis. Applied to both real and simulated RNA-seq data sets from Arabidopsis (Arabidopsis thaliana), PanGraphRNA outperformed the SLR approach, achieving higher read alignment accuracy and more precise gene expression quantification. PanGraphRNA enabled the identification of drought stress-induced genes and flowering time-related quantitative trait loci that were previously missed with the conventional SLR approach. Furthermore, we successfully applied PanGraphRNA to process RNA-seq data sets from rice (Oryza sativa) and maize (Zea mays). By providing standardized, containerized workflows, PanGraphRNA will facilitate transcriptomic research in key plant species, including Arabidopsis, rice, and maize.

The sequencing of the full transcriptome (RNA-seq) has become the preferred choice for the measurement of genome-wide gene expression. Despite its widespread use, challenges remain in RNA-seq data analysis. One often-overlooked aspect is normalization. Despite the fact that a variety of factors or ‘batch effects’ can contribute unwanted variation to the data, commonly used RNA-seq normalization methods only correct for sequencing depth. The study of gene expression is particularly problematic when it is influenced simultaneously by a variety of biological factors in addition to the one of interest. Using examples from experimental neuroscience, we show that batch effects can dominate the signal of interest; and that the choice of normalization method affects the power and reproducibility of the results. While commonly used global normalization methods are not able to adequately normalize the data, more recently developed RNA-seq normalization can. We focus on one particular method, RUVSeq and show that it is able to increase power and biological insight of the results. Finally, we provide a tutorial outlining the implementation of RUVSeq normalization that is applicable to a broad range of studies as well as meta-analysis of publicly available data.

Recent advances in high-throughput RNA sequencing (RNA-seq) have enabled tremendous leaps forward in our understanding of bacterial transcriptomes. However, computational methods for analysis of bacterial transcriptome data have not kept pace with the large and growing data sets generated by RNA-seq technology. Here, we present new algorithms, specific to bacterial gene structures and transcriptomes, for analysis of RNA-seq data. The algorithms are implemented in an open source software system called Rockhopper that supports various stages of bacterial RNA-seq data analysis, including aligning sequencing reads to a genome, constructing transcriptome maps, quantifying transcript abundance, testing for differential gene expression, determining operon structures and visualizing results. We demonstrate the performance of Rockhopper using 2.1 billion sequenced reads from 75 RNA-seq experiments conducted with Escherichia coli, Neisseria gonorrhoeae, Salmonella enterica, Streptococcus pyogenes and Xenorhabdus nematophila. We find that the transcriptome maps generated by our algorithms are highly accurate when compared with focused experimental data from E. coli and N. gonorrhoeae, and we validate our system’s ability to identify novel small RNAs, operons and transcription start sites. Our results suggest that Rockhopper can be used for efficient and accurate analysis of bacterial RNA-seq data, and that it can aid with elucidation of bacterial transcriptomes.

Introduction: Colon cancer is the fourth most common cancer in the United States and the third leading cause of cancer-related death. High throughput genomic sequencing has led to a number of significant advances in tumor biology and in identifying novel signaling molecules such as long non-coding RNA (lncRNA). The aim of this study was to identify differentially expressed lncRNA in colon cancer from a large RNA-sequencing (RNA-seq) data set. Methods: The raw RNA-seq files of 398 patients with colon cancer were downloaded from The Cancer Genome Atlas (TCGA). Sequencing files were aligned using STAR (Spliced Transcripts Alignment to a Reference) to the most recent genome annotation. Using a subset of patients with paired colon cancer and normal colon epithelium RNA-seq data (n=40), an exploratory binomial regression model was used to calculate differential RNA expression. The most differentially expressed lncRNA were identified from the exploratory analysis and verified by comparing the larger colon cancer RNA-seq data set (n=358) with the normal colon epithelium RNA-seq data set (n=40). Results: 33,514 genes were identified from the comparative analysis, and using differential expression cut off values of > +1.5 or < -1.5 log fold change and a false discovery rate of <0.05; 543 were upregulated and 1822 were downregulated. Within this dysregulated group, 60 lncRNAs were identified. Using the larger data set (n=358 vs n=40), 41/60 lncRNA remained differentially expressed, of which 15 were downregulated and 26 were upregulated. Twenty-four of these lncRNAs have not previously been described in colon cancer, 12 of which are upregulated and 12 of which are downregulated (Table 1). Conclusions: This analysis of RNA-seq data from TCGA has identified dysregulated lncRNAs which have not previously been described in human colon cancer. These lncRNAs may have significant roles in colon cancer tumor signaling. Table 1.Differential expression of lncRNA not described in colon cancerlncRNALog Fold ChangeP-valueDescription in other cancersKRT16P14.1041.36E-21Increased in lung squamous cell carcinomaBOK-AS13.7784.98E-09Increased in oral squamous cell carcinomaCELP3.1151.30E-09Not availableCLDN10-AS12.3944.61E-12Not availablePOU6F2-AS12.2678.25E-09Not availableSLC7A11-AS12.1931.60E-13Decreased in gastric adenocarcinomaCSAG42.1083.59E-10Not availableLUCAT12.0224.21E-11Increased in Ovarian cancer, renal clear cell carcinoma, head and neck squamous cell carcinoma, Non-small cell cancer, glioma, osteosarcoma, esophageal squamous cell carcinomaVPS9D1-AS11.9795.27E-21Decreased in gastric adenocarcinoma, Increased in non-small cell lung cancerLEF1-AS11.8464.48E-13Increased in glioblastomaCASC81.6485.64E-06Single nucleotide polymorphisms in colorectal adenocarcinomaLINC004911.5903.11E-04Not availableLINC00648-1.5021.02E-15Not availableLINC00163-1.6751.03E-22Not availableSATB2-AS1-1.8074.18E-09Increased in osteosarcomaLINC00702-2.1424.81E-18Not availableTRHDE-AS1-2.2623.44E-24Not availableLINC00908-2.3082.76E-103Not availableLINC00483-2.4012.87E-39Increased in gastric adenocarcinomaLINC00461-2.4401.42E-41Increased in glioma, multiple myelomaFENDRR-2.6872.12E-58Decreased in breast cancer, prostate cancer, osteosarcoma, gastric adenocarcinoma,ARHGEF26-AS1-2.7148.55E-32Not availableADAMTS9-AS2-2.9041.78E-74Increased in lung cancer salivary adenoid cystic carcinoma, Decreased in glioma, gastric adenocarcinomaMT1JP-4.4601.08E-114Decreased in retinoblastoma, Decreased in gastric adenocarcinoma Citation Format: Stephen J. O'Brien, Theodore Kalbfleisch, Sudhir Srivastava, Shesh Rai, Susan Galandiuk. Differential expression of long non-coding RNA in colon adenocarcinoma RNA-sequence data set [abstract]. In: Proceedings of the American Association for Cancer Research Annual Meeting 2019; 2019 Mar 29-Apr 3; Atlanta, GA. Philadelphia (PA): AACR; Cancer Res 2019;79(13 Suppl):Abstract nr 1817.

Normalization is crucial in RNA-seq data analyses. Due to the existence of excessive zeros and a large number of small measures, it is challenging to find reliable linear rescaling normalization parameters. We propose a Zipf plot based normalization method (ZN) assuming that all gene profiles have similar upper tail behaviors in their expression distributions. The new normalization method uses global information of all genes in the same profile without gene-level expression alteration. It doesn’t require the majority of genes to be not differentially expressed (DE), and can be applied to data where the majority of genes are weakly or not expressed. Two normalization schemes are implemented with ZN: a linear rescaling scheme and a non-linear transformation scheme. The linear rescaling scheme can be applied alone or together with the non-linear normalization scheme. The performance of ZN is benchmarked against five popular linear normalization methods for RNA-seq data. Results show that the linear rescaling normalization scheme by itself works well and is robust. The non-linear normalization scheme can further improve the normalization outcomes and is optional if the Zipf plots show parallel patterns.

BackgroundThe basidiomycete Moniliophthora roreri is the causal agent of Frosty pod rot (FPR) disease of cacao (Theobroma cacao), the source of chocolate, and FPR is one of the most destructive diseases of this important perennial crop in the Americas. This hemibiotroph infects only cacao pods and has an extended biotrophic phase lasting up to sixty days, culminating in plant necrosis and sporulation of the fungus without the formation of a basidiocarp.ResultsWe sequenced and assembled 52.3 Mb into 3,298 contigs that represent the M. roreri genome. Of the 17,920 predicted open reading frames (OFRs), 13,760 were validated by RNA-Seq. Using read count data from RNA sequencing of cacao pods at 30 and 60 days post infection, differential gene expression was estimated for the biotrophic and necrotrophic phases of this plant-pathogen interaction. The sequencing data were used to develop a genome based secretome for the infected pods. Of the 1,535 genes encoding putative secreted proteins, 1,355 were expressed in the biotrophic and necrotrophic phases. Analysis of the data revealed secretome gene expression that correlated with infection and intercellular growth in the biotrophic phase and invasive growth and plant cellular death in the necrotrophic phase.ConclusionsGenome sequencing and RNA-Seq was used to determine and validate the Moniliophthora roreri genome and secretome. High sequence identity between Moniliophthora roreri genes and Moniliophthora perniciosa genes supports the taxonomic relationship with Moniliophthora perniciosa and the relatedness of this fungus to other basidiomycetes. Analysis of RNA-Seq data from infected plant tissues revealed differentially expressed genes in the biotrophic and necrotrophic phases. The secreted protein genes that were upregulated in the biotrophic phase are primarily associated with breakdown of the intercellular matrix and modification of the fungal mycelia, possibly to mask the fungus from plant defenses. Based on the transcriptome data, the upregulated secreted proteins in the necrotrophic phase are hypothesized to be actively attacking the plant cell walls and plant cellular components resulting in necrosis. These genes are being used to develop a new understanding of how this disease interaction progresses and to identify potential targets to reduce the impact of this devastating disease.

The vast majority of human genes express more than one transcript isoform, and aberrant splicing patterns or changes in isoform expression ratios can lead to disease, notably cancer. Next generation sequencing platforms deliver the finest resolution yet of the transcriptome, and the wealth of sequencing data produced by TCGA provides an excellent basis for a broad study of alternative splicing in cancer. We have developed an interactive, exceptionally user-friendly tool, SpliceSeq, for investigation of alternative splicing in RNASeq data. The program is a follow-up to our SpliceCenter package (http://projects.insilico.us/SpliceCenter/SpliceOverview.jsp), which provides tools for use by the bench biologist in assessing experimental observations that may result from splice variation. Our results from SpliceSeq for TCGA RNASeq data from several cancer types will be presented. We identify pervasive changes in splicing patterns across different cancer types by comparing tumor vs. normal and clinical subgroups (e.g. tumor stage, hyper-mutator status, and MSI designation). Interesting observations include striking splice variations in APP, the amyloid beta precursor protein, and in PKM, which is centrally involved in the balance between oxidative and glycolytic metabolism in tumors. SpliceSeq is freely available at http://bioinformatics.mdanderson.org/main/SpliceSeq:Overview. 1. Ryan MC, Cleland J, Kim R, Wong WC, Weinstein JN. SpliceSeq: A Resource for Analysis and Visualization of RNA-Seq Data on Alternative Splicing and Its Functional Impacts. Bioinformatics, 10.1093, 2012. 2. Ryan MC, Zeeberg BR, Caplen NJ, Cleland JA, Kahn AB, Liu H, and Weinstein JN. SpliceCenter: a suite of web-based bioinformatic applications for evaluating the impact of alternative splicing on RT-PCR, RNAi, microarray, and peptide-based studies. BMC Bioinformatics 9:313, 2008. 3. Kahn AB, Ryan MC, Liu H, Zeeberg BR, Jamison DC, Weinstein JN. SpliceMiner: a high-throughput database implementation of the NCBI Evidence Viewer for microarray splice variant analysis. BMC Bioinformatics 8:752, 2006. Supported by NCI Grant No. U24CA143883. The content is solely the responsibility of the authors and does not necessarily represent official views of NCI or NIH. Also supported in part by a gift from the Mary K. Chapman Foundation and a grant from the Michael and Susan Dell Foundation honoring Lorraine Dell. Citation Format: Michael C. Ryan, John N. Weinstein. Analysis of TCGA RNASeq data using SpliceSeq provides a survey of alternative splicing in cancer. [abstract]. In: Proceedings of the 104th Annual Meeting of the American Association for Cancer Research; 2013 Apr 6-10; Washington, DC. Philadelphia (PA): AACR; Cancer Res 2013;73(8 Suppl):Abstract nr 1796. doi:10.1158/1538-7445.AM2013-1796

The rate of raw sequence production through Next-Generation Sequencing (NGS) has been growing exponentially due to improved technology and reduced costs. This has enabled researchers to answer many biological questions through “multi-omics” data analyses. Even though such data promises new insights into how biological systems function and understanding disease mechanisms, computational analyses performed on such large datasets comes with its challenges and potential pitfalls. The aim of this study was to develop a robust portable and reproducible bioinformatic pipeline for the automation of RNA sequencing (RNA-seq) data analyses. Using Nextflow as a workflow management system and Singularity for application containerisation, the nf-rnaSeqCount pipeline was developed for mapping raw RNA-seq reads to a reference genome and quantifying abundance of identified genomic features for differential gene expression analyses. The pipeline provides a quick and efficient way to obtain a matrix of read counts that can be used with tools such as DESeq2 and edgeR for differential expression analysis. Robust and flexible bioinformatic and computational pipelines for RNA-seq data analysis, from QC to sequence alignment and comparative analyses, will reduce analysis time, and increase accuracy and reproducibility of findings to promote transcriptome research.

This manuscript describes the development of a resource module that is part of a learning platform named “NIGMS Sandbox for Cloud-based Learning” https://github.com/NIGMS/NIGMS-Sandbox. The overall genesis of the Sandbox is described in the editorial NIGMS Sandbox at the beginning of this Supplement. This module delivers learning materials on RNA sequencing (RNAseq) data analysis in an interactive format that uses appropriate cloud resources for data access and analyses. Biomedical research is increasingly data-driven, and dependent upon data management and analysis methods that facilitate rigorous, robust, and reproducible research. Cloud-based computing resources provide opportunities to broaden the application of bioinformatics and data science in research. Two obstacles for researchers, particularly those at small institutions, are: (i) access to bioinformatics analysis environments tailored to their research; and (ii) training in how to use Cloud-based computing resources. We developed five reusable tutorials for bulk RNAseq data analysis to address these obstacles. Using Jupyter notebooks run on the Google Cloud Platform, the tutorials guide the user through a workflow featuring an RNAseq dataset from a study of prophage altered drug resistance in Mycobacterium chelonae. The first tutorial uses a subset of the data so users can learn analysis steps rapidly, and the second uses the entire dataset. Next, a tutorial demonstrates how to analyze the read count data to generate lists of differentially expressed genes using R/DESeq2. Additional tutorials generate read counts using the Snakemake workflow manager and Nextflow with Google Batch. All tutorials are open-source and can be used as templates for other analysis.

Measuring differential gene expression is a common task in the analysis of RNA-Seq data. To identify differentially expressed genes between two samples, it is crucial to normalize the datasets. While multiple normalization methods are available, all of them are based on certain assumptions that may or may not be suitable for the type of data they are applied on. Researchers therefore need to select an adequate normalization strategy for each RNA-Seq experiment. This selection includes exploration of different normalization methods as well as their comparison. Methods that agree with each other most likely represent realistic assumptions under the particular experimental conditions. We developed the NVT package, which provides a fast and simple way to analyze and evaluate multiple normalization methods via visualization and representation of correlation values, based on a user-defined set of uniformly expressed genes. The R package is freely available under https://github.com/Edert/NVT CONTACT: thomas.rattei@univie.ac.atSupplementary information: Supplementary data are available at Bioinformatics online.

IntroductionOur study undertakes a detailed exploration of gene expression dynamics within human lung organ tissue equivalents (OTEs) in response to Influenza A virus (IAV), Human metapneumovirus (MPV), and Parainfluenza virus type 3 (PIV3) infections. Through the analysis of RNA-Seq data from 19,671 genes, we aim to identify differentially expressed genes under various infection conditions, elucidating the complexities of virus-host interactions.MethodsWe employ Generalized Linear Models (GLMs) with Quasi-Likelihood (QL) F-tests (GLMQL) and introduce the novel Magnitude-Altitude Score (MAS) and Relaxed Magnitude-Altitude Score (RMAS) algorithms to navigate the intricate landscape of RNA-Seq data. This approach facilitates the precise identification of potential biomarkers, highlighting the host’s reliance on innate immune mechanisms. Our comprehensive methodological framework includes RNA extraction, library preparation, sequencing, and Gene Ontology (GO) enrichment analysis to interpret the biological significance of our findings.ResultsThe differential expression analysis unveils significant changes in gene expression triggered by IAV, MPV, and PIV3 infections. The MAS and RMAS algorithms enable focused identification of biomarkers, revealing a consistent activation of interferon-stimulated genes (e.g., IFIT1, IFIT2, IFIT3, OAS1) across all viruses. Our GO analysis provides deep insights into the host’s defense mechanisms and viral strategies exploiting host cellular functions. Notably, changes in cellular structures, such as cilium assembly and mitochondrial ribosome assembly, indicate a strategic shift in cellular priorities. The precision of our methodology is validated by a 92% mean accuracy in classifying respiratory virus infections using multinomial logistic regression, demonstrating the superior efficacy of our approach over traditional methods.DiscussionThis study highlights the intricate interplay between viral infections and host gene expression, underscoring the need for targeted therapeutic interventions. The stability and reliability of the MAS/RMAS ranking method, even under stringent statistical corrections, and the critical importance of adequate sample size for biomarker reliability are significant findings. Our comprehensive analysis not only advances our understanding of the host’s response to viral infections but also sets a new benchmark for the identification of biomarkers, paving the way for the development of effective diagnostic and therapeutic strategies.

Long noncoding RNAs (lncRNAs) are emerging as key regulators of diverse cellular processes, but their roles in breast cancer biology are just beginning to be elucidated. In this study, integration of powerful techniques, RNA-seq data from subcellular fractionated RNA with GRO-seq data has yielded a comprehensive catalog of estrogen-regulated lncRNAs in MCF-7 cells. Analysis of RNA-seq data from samples representing molecular subtypes of breast cancer and normal tissue types, revealed that many lncRNAs (such as lincRNA86) show distinct expression patterns. LincRNA86 shows highest normal expression in pancreas followed by testis in normal human tissues. The hypothesis is lincRNA86 regulate estrogen-dependent signaling in breast cancer. In functional assays, knockdown of lncRNA86 inhibits the growth of ER-positive breast cancer cells. Amplified expression of lncRNA86 in breast cancer correlates with clinical outcome. LncRNA86 have now been fully annotated (transcription start and stop site, 5’ cap, polyA tail, and exon/intron structure), and cloned. We are now performing detailed molecular analyses to better understand the underlying mechanisms of action of the lncRNA. We are also currently have experiments underway to view cancer phenotypes: estrogen-dependent tumor growth. Collectively, our preliminary results suggest that lincRNA86 plays a critical role in ERα-dependent pathways.

BackgroundThe use of RNA-sequencing (RNA-seq) in molecular biology research and clinical settings has increased significantly over the past decade. Despite its widespread adoption, there is a lack of simple and interactive tools to analyze and explore RNA-seq data. Many established tools require programming or Unix/Bash knowledge to analyze and visualize results. This requirement presents a significant barrier for many researchers to efficiently analyze and present RNA-seq data.ResultsHere we present BEAVR, a Browser-based tool for the Exploration And Visualization of RNA-seq data. BEAVR is an easy-to-use tool that facilitates interactive analysis and exploration of RNA-seq data. BEAVR is developed in R and uses DESeq2 as its engine for differential gene expression (DGE) analysis, but assumes users have no prior knowledge of R or DESeq2. BEAVR allows researchers to easily obtain a table of differentially-expressed genes with statistical testing and then visualize the results in a series of graphs, plots and heatmaps. Users are able to customize many parameters for statistical testing, dealing with variance, clustering methods and pathway analysis to generate high quality figures.ConclusionBEAVR simplifies analysis for novice users but also streamlines the RNA-seq analysis process for experts by automating several steps. BEAVR and its documentation can be found on GitHub at https://github.com/developerpiru/BEAVR. BEAVR is available as a Docker container at https://hub.docker.com/r/pirunthan/beavr.

p53 is inactivated in approximately half of human cancers. p53 execute various functions by modulating transcriptional regulation. To identify direct p53 transcriptional targets including non-coding RNAs, we performed gene expression analysis with next-generation sequencing (RNA-seq) in osteosarcoma U2OS cells treated with Nutlin-3a which activates endogenous p53. We identified 261 genes increased more than 4-fold and 373 long non-coding RNAs (lncRNAs) increased more than 2-fold. Next, we analyzed the data of chromatin immunoprecipitation with next-generation sequencing (ChIP-seq) and then searched for p53 binding sites in the whole human genome. p53 binding sites were located in the neighborhood of 126/261 (48.3%) genes including 28 known target genes and 122/373 (32.7%) lncRNAs including 7 lncRNAs we already identified in previous study. Recent reports have revealed that lncRNAs play an important role in various biological and pathological processes such as development, differentiation, stemness and carcinogenesis. Although all functions of lncRNAs have not been elucidated, many lncRNAs are involved in transcriptional regulation like transcriptional factors. Therefore, we performed a gene network analysis of cancer transcriptome including lncRNAs using public RNA-seq data of cancer tissues. As a result, we identified several hub lncRNAs which regulate the transcription of many downstream genes. Our results suggest that p53 and lncRNAs comprise a complex transcriptional network for various biological functions and tumor suppression. Citation Format: Masashi Idogawa, Tomoko Ohashi, Yasushi Sasaki, Takashi Tokino. The identification of p53 target genes and noncoding RNAs through the combined analysis of RNA-seq and ChIP-seq data. [abstract]. In: Proceedings of the 107th Annual Meeting of the American Association for Cancer Research; 2016 Apr 16-20; New Orleans, LA. Philadelphia (PA): AACR; Cancer Res 2016;76(14 Suppl):Abstract nr 1938.

RNA sequencing (RNA-seq) is a highly sensitive method for transcriptome analysis that allows simultaneous assessment of expression of thousands of genes and identification of expression patterns under various conditions. The existing variety of RNA-seq data formats, normalization methods, and approaches to statistical processing of results complicates comparison of data from different studies and reduces reproducibility of the analysis. This study presents an automated pipeline PipeSeq that combines standard steps of RNA-seq data processing: loading (SRA Toolkit), read alignment to the reference genome (HISAT2), transcript assembly (StringTie), transcript counting (FeatureCounts) and statistical analysis of differential gene expression under various experimental conditions (DESeq2). PipeSeq has a simple visual interface, supports multithreading, and generates ready-to-analyze gene expression heat maps, tables and graphs. The functionality of the pipeline is demonstrated on three sets of raw RNA-seq data from the green alga Chlamydomonas reinhardtii cells available in the NCBI SRA database. The data from these experiments were used to analyze the differential expression of C. reinhardtii genes encoding the GATA family transcription factors under different light cultivation conditions. The data obtained by in silico methods were verified by real-time reverse transcription polymerase chain reaction (RT-qPCR) for 12 GATA genes, which allowed us to hypothesize their functions and evaluate the correlation between the bulk (RNA-seq) and targeted (RT-qPCR) approaches. Our results showed that RNA-seq and RT-qPCR methods reveal similar directions of gene expression changes, but demonstrate differences in the effect size and sensitivity, which emphasizes the need for a combined use of the two approaches. Thus, the PipeSeq program is a tool for conducting a full cycle of bioinformatic analysis of RNA-seq data, additionally providing the opportunity to process RT-qPCR data and perform a comparative statistical analysis of the results obtained.