Nf-rnaSeqCount: A Nextflow pipeline for obtaining raw read counts from RNA-seq data.

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.

BackgroundRNA sequencing (RNA-seq) is a powerful tool for genome-wide expression profiling of biological samples with the advantage of high-throughput and high resolution. There are many existing algorithms nowadays for quantifying expression levels and detecting differential gene expression, but none of them takes the misaligned reads that are mapped to non-exonic regions into account. We developed a novel algorithm, XBSeq, where a statistical model was established based on the assumption that observed signals are the convolution of true expression signals and sequencing noises. The mapped reads in non-exonic regions are considered as sequencing noises, which follows a Poisson distribution. Given measureable observed and noise signals from RNA-seq data, true expression signals, assuming governed by the negative binomial distribution, can be delineated and thus the accurate detection of differential expressed genes.ResultsWe implemented our novel XBSeq algorithm and evaluated it by using a set of simulated expression datasets under different conditions, using a combination of negative binomial and Poisson distributions with parameters derived from real RNA-seq data. We compared the performance of our method with other commonly used differential expression analysis algorithms. We also evaluated the changes in true and false positive rates with variations in biological replicates, differential fold changes, and expression levels in non-exonic regions. We also tested the algorithm on a set of real RNA-seq data where the common and different detection results from different algorithms were reported.ConclusionsIn this paper, we proposed a novel XBSeq, a differential expression analysis algorithm for RNA-seq data that takes non-exonic mapped reads into consideration. When background noise is at baseline level, the performance of XBSeq and DESeq are mostly equivalent. However, our method surpasses DESeq and other algorithms with the increase of non-exonic mapped reads. Only in very low read count condition XBSeq had a slightly higher false discovery rate, which may be improved by adjusting the background noise effect in this situation. Taken together, by considering non-exonic mapped reads, XBSeq can provide accurate expression measurement and thus detect differential expressed genes even in noisy conditions.

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.

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.

In this chapter, we describe methods for analyzing RNA-Seq data, presented as a flow along a pipeline beginning with raw data from a sequencer and ending with an output of differentially expressed genes and their functional characterization. The first section covers de novo transcriptome assembly for organisms lacking reference genomes or for those interested in probing against the background of organism-specific transcriptomes assembled from RNA-Seq data. Section 2 covers both gene- and transcript-level quantifications, leading to the third and final section on differential expression analysis between two or more conditions. The pipeline starts with raw sequence reads, followed by quality assessment and preprocessing of the input data to ensure a robust estimate of the transcripts and their differential regulation. The preprocessed data can be inputted into the de novo transcriptome flow to assemble transcripts, functionally annotated using tools such as InterProScan or Blast2Go and then forwarded to differential expression analysis flow, or directly inputted into the differential expression analysis flow if a reference genome is available. An online repository containing sample data has also been made available, as well as custom Python scripts to modify the output of the programs within the pipeline for various downstream analyses.

BackgroundDifferential gene expression analysis is a prominent technique for the analysis of biomolecular data to identify genetic features associated with phenotypes. Limma—for microarray data –, and edgeR and DESeq2—for RNA-Seq data–, are the most widely used tools for differential gene expression analysis of bulk transcriptomic data.ResultsWe present the differential expression features of InMoose, a Python implementation of R tools limma, edgeR, and DESeq2. We experimentally show that InMoose stands as a drop-in replacement for those tools, with nearly identical results. This ensures reproducibility when interfacing both languages in bioinformatic pipelines. InMoose is an open source software released under the GPL3 license, available at www.github.com/epigenelabs/inmoose and https://inmoose.readthedocs.io.ConclusionsWe present a new Python implementation of state-of-the-art tools limma, edgeR, and DESeq2, to perform differential gene expression analysis of bulk transcriptomic data. This new implementation exhibits results nearly identical to the original tools, improving interoperability and reproducibility between Python and R bioinformatics pipelines.

In differential gene expression data analysis, one objective is to identify groups of co-expressed genes from a large dataset in order to detect the association between such a group of genes and an experimental condition. This is often done through a clustering approach, such as k-means or bipartition hierarchical clustering, based on particular similarity measures in the grouping process. In such a dataset, the gene differential expression itself is an innate attribute that can be used in the feature extraction process. For example, in a dataset consisting of multiple treatments versus their controls, the expression of a gene in each treatment would have three possible behaviors, upregulated, downregulated, or unchanged. We present in this chapter, a differential expression feature extraction (DEFE) method by using a string consisting of three numerical values at each character to denote such behavior, i.e., 1=up, 2=down, and 0=unchanged, which results in up to 3B differential expression patterns across all B comparisons. This approach has been successfully applied in many research projects, and among these, we demonstrate the strength of DEFE in a case study on RNA-sequencing (RNA-seq) data analysis of wheat challenged with the phytopathogenic fungus, Fusarium graminearum. Combinations of multiple schemes of DEFE patterns revealed groups of genes putatively associated with resistance or susceptibility to FHB.

Many scientific endeavors, such as molecular biology, have become dependent on big data and its analysis. For example, precision medicine depends on molecular measurements and data analysis per patient. Data analyses supporting medical decisions must be standardized and performed consistently across patients. While perhaps not life-threatening, data analyses in basic research have become increasingly complex. RNAseq data, for example, entails a multi-step analysis ranging from quality assessment of the measurements to statistical analyses. Workflow management systems (WFMS) enable the development of data analysis workflows (WF), their reproduction, and their application to datasets of the same type. However, far more than a hundred WFMS are available, and there is no way to convert data analysis WFs among WFMS. Therefore, the initial choice of a WFMS is important as it entails a lockin to the system. The reach in their particular field (number of citations) can be used as a proxy for selecting a WFMS, but of the about 25 WFMS we mention in this work, at least 5 have a large reach in scientific data analysis. Hence other criteria are needed to delineate among WFMS. By extracting such criteria from selected studies concerning WFMS and adding additional criteria, we arrived at five critical criteria: reproducibility, reusability, FAIRness, versioning support, and security. Another five criteria (providing a graphical user interface, WF flexibility, WF scalability, WF shareability, and computational transparency) we deemed important but not critical for the assessment of WFMS. We applied the criteria to the most cited WFMS in PubMed and found none that support all criteria. We hope that suggesting these criteria will spark a discussion on what features are important for WFMS in scientific data analysis and may lead to developing WFMS that fulfill such criteria.

RNA-Seq is increasingly being used for differential gene expression analysis which was dominated by the microarray technology in the past decade. However, inferring differential gene expression based on the observed difference of RNA-Seq read counts has unique challenges that were not present in microarray-based analysis. The differential expression estimation may be biased against low read count values such that the differential expression of genes with high read counts is more easily detected. The estimation bias may further propagate in downstream analyses at the systems biology level if it is not corrected. To obtain a better inference of differential gene expression, we propose a new efficient algorithm based on a markov random field (MRF) model, called MRFSeq, that uses additional gene coexpression data to enhance the prediction power. Our main technical contribution is the careful selection of the clique potential functions in the MRF so its maximum a posteriori (MAP) estimation can be reduced to the well-known maximum flow problem and thus solved in polynomial time. Our extensive experiments on simulated and real RNA-Seq datasets demonstrate that MRFSeq is more accurate and less biased against genes with low read counts than the existing methods based on RNA-Seq data alone. For example, on the well-studied MAQC dataset, MRFSeq improved the sensitivity from 11.6% to 38.8% for genes with low read counts. MRFSeq is implemented in C++ and available at http://www.cs.ucr.edu/~yyang027/mrfseq.htm.

BackgroundRecently, rapid improvements in technology and decrease in sequencing costs have made RNA-Seq a widely used technique to quantify gene expression levels. Various normalization approaches have been proposed, owing to the importance of normalization in the analysis of RNA-Seq data. A comparison of recently proposed normalization methods is required to generate suitable guidelines for the selection of the most appropriate approach for future experiments.ResultsIn this paper, we compared eight non-abundance (RC, UQ, Med, TMM, DESeq, Q, RPKM, and ERPKM) and two abundance estimation normalization methods (RSEM and Sailfish). The experiments were based on real Illumina high-throughput RNA-Seq of 35- and 76-nucleotide sequences produced in the MAQC project and simulation reads. Reads were mapped with human genome obtained from UCSC Genome Browser Database. For precise evaluation, we investigated Spearman correlation between the normalization results from RNA-Seq and MAQC qRT-PCR values for 996 genes. Based on this work, we showed that out of the eight non-abundance estimation normalization methods, RC, UQ, Med, TMM, DESeq, and Q gave similar normalization results for all data sets. For RNA-Seq of a 35-nucleotide sequence, RPKM showed the highest correlation results, but for RNA-Seq of a 76-nucleotide sequence, least correlation was observed than the other methods. ERPKM did not improve results than RPKM. Between two abundance estimation normalization methods, for RNA-Seq of a 35-nucleotide sequence, higher correlation was obtained with Sailfish than that with RSEM, which was better than without using abundance estimation methods. However, for RNA-Seq of a 76-nucleotide sequence, the results achieved by RSEM were similar to without applying abundance estimation methods, and were much better than with Sailfish. Furthermore, we found that adding a poly-A tail increased alignment numbers, but did not improve normalization results.ConclusionSpearman correlation analysis revealed that RC, UQ, Med, TMM, DESeq, and Q did not noticeably improve gene expression normalization, regardless of read length. Other normalization methods were more efficient when alignment accuracy was low; Sailfish with RPKM gave the best normalization results. When alignment accuracy was high, RC was sufficient for gene expression calculation. And we suggest ignoring poly-A tail during differential gene expression analysis.Electronic supplementary materialThe online version of this article (doi:10.1186/s12859-015-0778-7) contains supplementary material, which is available to authorized users.

Background. A common research goal in transcriptome projects is to find genes that are differentially expressed in different phenotype classes. Biologists might wish to validate such gene candidates experimentally, or use them for downstream systems biology analysis. Producing a coherent differential gene expression analysis from RNA-seq count data requires an understanding of how numerous sources of variation such as the replicate size, the hypothesized biological effect size, and the specific method for making differential expression calls interact. We believe an explicit demonstration of such interactions in real RNA-seq data sets is of practical interest to biologists.Results. Using two large public RNA-seq data sets—one representing strong, and another mild, biological effect size—we simulated different replicate size scenarios, and tested the performance of several commonly-used methods for calling differentially expressed genes in each of them. We found that, when biological effect size was mild, RNA-seq experiments should focus on experimental validation of differentially expressed gene candidates. Importantly, at least triplicates must be used, and the differentially expressed genes should be called using methods with high positive predictive value (PPV), such as NOISeq or GFOLD. In contrast, when biological effect size was strong, differentially expressed genes mined from unreplicated experiments using NOISeq, ASC and GFOLD had between 30 to 50% mean PPV, an increase of more than 30-fold compared to the cases of mild biological effect size. Among methods with good PPV performance, having triplicates or more substantially improved mean PPV to over 90% for GFOLD, 60% for DESeq2, 50% for NOISeq, and 30% for edgeR. At a replicate size of six, we found DESeq2 and edgeR to be reasonable methods for calling differentially expressed genes at systems level analysis, as their PPV and sensitivity trade-off were superior to the other methods’.Conclusion. When biological effect size is weak, systems level investigation is not possible using RNAseq data, and no meaningful result can be obtained in unreplicated experiments. Nonetheless, NOISeq or GFOLD may yield limited numbers of gene candidates with good validation potential, when triplicates or more are available. When biological effect size is strong, NOISeq and GFOLD are effective tools for detecting differentially expressed genes in unreplicated RNA-seq experiments for qPCR validation. When triplicates or more are available, GFOLD is a sharp tool for identifying high confidence differentially expressed genes for targeted qPCR validation; for downstream systems level analysis, combined results from DESeq2 and edgeR are useful.

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.

The study of differential gene expression patterns through RNA-Seq comprises a routine task in the daily lives of molecular bioscientists, who produce vast amounts of data requiring proper management and analysis. Despite widespread use, there are still no widely accepted golden standards for the normalization and statistical analysis of RNA-Seq data, and critical biases, such as gene lengths and problems in the detection of certain types of molecules, remain largely unaddressed. Stimulated by these unmet needs and the lack of in-depth research into the potential of combinatorial methods to enhance the analysis of differential gene expression, we had previously introduced the PANDORA P-value combination algorithm while presenting evidence for PANDORA's superior performance in optimizing the tradeoff between precision and sensitivity. In this article, we present the next generation of the algorithm along with a more in-depth investigation of its capabilities to effectively analyze RNA-Seq data. In particular, we show that PANDORA-reported lists of differentially expressed genes are unaffected by biases introduced by different normalization methods, while, at the same time, they comprise a reliable input option for downstream pathway analysis. Additionally, PANDORA outperforms other methods in detecting differential expression patterns in certain transcript types, including long non-coding RNAs.

Identification of a Relapse Signature in t(8;21) Pediatric AML

Transcriptome sequencing, where RNA is isolated, converted to library of cDNA fragments, and sequenced using next-generation sequencing technology, has become the method of choice for the genome-wide characterization of mRNA levels. It offers a more accurate quantification of transcript levels than array-based methods, but also has the added benefit of allowing the discovery of novel gene/transcripts, alternative splice junctions, and novel RNAs. In addition, RNA sequencing may be used to investigate differential gene expression, allelic imbalance, eQTL mapping, RNA editing, RNA-protein interactions, and alternative splicing. A number of statistical methods and tools are available for differential expression analysis using RNA sequencing data and these are continually being developed and improved to handle more complex experimental designs. This chapter describes an example workflow for the quality control and analysis of raw RNA sequencing reads for the purposes of differential gene expression analysis, followed by pathway/enrichment analysis of significantly different genes. The methods and tools described are just one example of how this analysis can be conducted, but they can be applied to most standard RNA sequencing studies of differential gene expression. The methods covered are based on Illumina HiSeq single-end 50 bp reads. However, all programs used are capable of working with paired-end data, subsequent to minor adaptations.