RNA-seq Differential Expression Analysis Research Articles

Multiple model species selection for transcriptomics analysis of non-model organisms

BackgroundTranscriptomic sequencing (RNA-seq) related applications allow for rapid explorations due to their high-throughput and relatively fast experimental capabilities, providing unprecedented progress in gene functional annotation, gene regulation analysis, and environmental factor verification. However, with increasing amounts of sequenced reads and reference model species, the selection of appropriate reference species for gene annotation has become a new challenge.MethodsWe proposed a novel approach for finding the most effective reference model species through taxonomic associations and ultra-conserved orthologous (UCO) gene comparisons among species. An online system for multiple species selection (MSS) for RNA-seq differential expression analysis was developed, and comprehensive genomic annotations from 291 reference model eukaryotic species were retrieved from the RefSeq, KEGG, and UniProt databases.ResultsUsing the proposed MSS pipeline, gene ontology and biological pathway enrichment analysis can be efficiently achieved, especially in the case of transcriptomic analysis of non-model organisms. The results showed that the proposed method solved problems related to limitations in annotation information and provided a roughly twenty-fold reduction in computational time, resulting in more accurate results than those of traditional approaches of using a single model reference species or the large non-redundant reference database.ConclusionsSelection of appropriate reference model species helps to reduce missing annotation information, allowing for more comprehensive results than those obtained with a single model reference species. In addition, adequate model species selection reduces the computational time significantly while retaining the same order of accuracy. The proposed system indeed provides superior performance by selecting appropriate multiple species for transcriptomic analysis compared to traditional approaches.

BEL1-LIKE HOMEODOMAIN 11 regulates chloroplast development and chlorophyll synthesis in tomato fruit.

Chloroplast development and chlorophyll(Chl)metabolism in unripe tomato contribute to the growth and quality of the fruit, however these mechanisms are poorly understood. In this study, we initially investigated seven homeobox-containing transcription factors (TFs) with specific ripening-associated expression patterns using virus-induced gene silencing (VIGS) technology and found that inhibiting the expression of one of these TFs, BEL1-LIKE HOMEODOMAIN11 (SlBEL11), significantly increased Chl levels in unripe tomato fruit. This enhanced Chl accumulation was further validated by generating stable RNA interference (RNAi) transgenic lines. RNA sequencing (RNA-seq) of RNAi-SlBEL11 fruit at the mature green (MG) stage showed that 48 genes involved in Chl biosynthesis, photosynthesis and chloroplast development were significantly upregulated compared with the wild type (WT) fruit. Genomic global scanning for Homeobox TF binding sites combined with RNA-seq differential gene expression analysis showed that 22 of these 48 genes were potential target genes of SlBEL11 protein. These genes included Chl biosynthesis-related genes encoding for protochlorophyllide reductase (POR), magnesium chelatase H subunit (CHLH) and chlorophyllide a oxygenase (CAO), and chloroplast development-related genes encoding for chlorophyll a/b binding protein (CAB), homeobox protein knotted 2 (TKN2) and ARABIDOPSIS PSEUDO RESPONSE REGULATOR 2-LIKE (APRR2-like). Electrophoretic mobility shift assay (EMSA) and chromatin immunoprecipitation quantitative polymerase chain reaction (PCR) (ChIP-qPCR) assays were employed to verify that SlBEL11 protein could bind to the promoters for TKN2, CAB and POR. Taken together, our findings demonstrated that SlBEL11 plays an important role in chloroplast development and Chl synthesis in tomato fruit.

Optimization of an RNA-Seq Differential Gene Expression Analysis Depending on Biological Replicate Number and Library Size.

RNA-Seq is a widely used technology that allows an efficient genome-wide quantification of gene expressions for, for example, differential expression (DE) analysis. After a brief review of the main issues, methods and tools related to the DE analysis of RNA-Seq data, this article focuses on the impact of both the replicate number and library size in such analyses. While the main drawback of previous relevant studies is the lack of generality, we conducted both an analysis of a two-condition experiment (with eight biological replicates per condition) to compare the results with previous benchmark studies, and a meta-analysis of 17 experiments with up to 18 biological conditions, eight biological replicates and 100 million (M) reads per sample. As a global trend, we concluded that the replicate number has a larger impact than the library size on the power of the DE analysis, except for low-expressed genes, for which both parameters seem to have the same impact. Our study also provides new insights for practitioners aiming to enhance their experimental designs. For instance, by analyzing both the sensitivity and specificity of the DE analysis, we showed that the optimal threshold to control the false discovery rate (FDR) is approximately 2−r, where r is the replicate number. Furthermore, we showed that the false positive rate (FPR) is rather well controlled by all three studied R packages: DESeq, DESeq2, and edgeR. We also analyzed the impact of both the replicate number and library size on gene ontology (GO) enrichment analysis. Interestingly, we concluded that increases in the replicate number and library size tend to enhance the sensitivity and specificity, respectively, of the GO analysis. Finally, we recommend to RNA-Seq practitioners the production of a pilot data set to strictly analyze the power of their experimental design, or the use of a public data set, which should be similar to the data set they will obtain. For individuals working on tomato research, on the basis of the meta-analysis, we recommend at least four biological replicates per condition and 20 M reads per sample to be almost sure of obtaining about 1000 DE genes if they exist.

Analysis of RNA-Seq Data Using TEtranscripts.

Transposable elements (TE) are mobile genetic elements that can readily change their genomic position. When not properly silenced, TEs can contribute a substantial portion to the cell's transcriptome, but are typically ignored in most RNA-seq data analyses. One reason for leaving TE-derived reads out of RNA-seq analyses is the complexities involved in properly aligning short sequencing reads to these highly repetitive regions. Here we describe a method for including TE-derived reads in RNA-seq differential expression analysis using an open source software package called TEtranscripts. TEtranscripts is designed to assign both uniquely and ambiguously mapped reads to all possible gene and TE-derived transcripts in order to statistically infer the correct gene/TE abundances. Here, we provide a detailed tutorial of TEtranscripts using a published qPCR validated dataset.Barbara McClintock laid the foundation for TE research with her discoveries in maize of mobile genetic elements capable of inserting into novel locations in the genome, altering the expression of nearby genes [1]. Since then, our appreciation of the contribution of repetitive TE-derived sequences to eukaryotic genomes has vastly increased. With the publication of the first human genome draft by the Human Genome Project, it was determined that nearly half of the human genome is derived from TE sequences [2, 3], with varying levels of repetitive DNA present in most plant and animal species. More recent studies looking at distantly related TE-like sequences have estimated that up to two thirds of the human genome might be repeat-derived [4], with the vast majority of these sequences attributed to retrotransposons that require transcription as part of the mobilization process, as discussed below.

RNA-Seq differential expression analysis: An extended review and a software tool.

The correct identification of differentially expressed genes (DEGs) between specific conditions is a key in the understanding phenotypic variation. High-throughput transcriptome sequencing (RNA-Seq) has become the main option for these studies. Thus, the number of methods and softwares for differential expression analysis from RNA-Seq data also increased rapidly. However, there is no consensus about the most appropriate pipeline or protocol for identifying differentially expressed genes from RNA-Seq data. This work presents an extended review on the topic that includes the evaluation of six methods of mapping reads, including pseudo-alignment and quasi-mapping and nine methods of differential expression analysis from RNA-Seq data. The adopted methods were evaluated based on real RNA-Seq data, using qRT-PCR data as reference (gold-standard). As part of the results, we developed a software that performs all the analysis presented in this work, which is freely available at https://github.com/costasilvati/consexpression. The results indicated that mapping methods have minimal impact on the final DEGs analysis, considering that adopted data have an annotated reference genome. Regarding the adopted experimental model, the DEGs identification methods that have more consistent results were the limma+voom, NOIseq and DESeq2. Additionally, the consensus among five DEGs identification methods guarantees a list of DEGs with great accuracy, indicating that the combination of different methods can produce more suitable results. The consensus option is also included for use in the available software.

Multiple functions of CREB-binding protein during postembryonic development: identification of target genes

BackgroundJuvenile hormones (JH) and ecdysteroids control postembryonic development in insects. They serve as valuable targets for pest management. Hence, understanding the molecular mechanisms of their action is of crucial importance. CREB-binding protein (CBP) is a universal transcriptional co-regulator. It controls the expression of several genes including those from hormone signaling pathways through co-activation of many transcription factors. However, the role of CBP during postembryonic development in insects is not well understood. Therefore, we have studied the role of CBP in postembryonic development in Tribolium, a model coleopteran insect.ResultsCBP is ubiquitously expressed in the red flour beetle, Tribolium castaneum. RNA interference (RNAi) mediated knockdown of CBP resulted in a decrease in JH induction of Kr-h1 gene expression in Tribolium larvae and led to a block in their development. Moreover, the injection of CBP double-stranded RNA (dsRNA) showed lethal phenotypes within 8 days of injection. RNA-seq and subsequent differential gene expression analysis identified CBP target genes in Tribolium. Knockdown of CBP caused a decrease in the expression of 1306 genes coding for transcription factors and other proteins associated with growth and development. Depletion of CBP impaired the expression of several JH response genes (e.g., Kr-h1, Hairy, early trypsin) and ecdysone response genes (EcR, E74, E75, and broad complex). Further, GO enrichment analyses of the downregulated genes showed enrichment in different functions including developmental processes, pigmentation, anatomical structure development, regulation of biological and cellular processes, etc.ConclusionThese data suggest diverse but crucial roles for CBP during postembryonic development in the coleopteran model insect, Tribolium. It can serve as a target for RNAi mediated pest management of this stored product pest.

A novel feature selection for RNA-seq analysis

RNA-seq data are challenging existing omics data analytics for its volume and complexity. Although quite a few computational models were proposed from different standing points to conduct differential expression (D.E.) analysis, almost all these methods do not provide a rigorous feature selection for high-dimensional RNA-seq count data. Instead, most or even all genes are invited into differential calls no matter they have real contributions to data variations or not. Thus, it would inevitably affect the robustness of D.E. analysis and lead to the increase of false positive ratios.In this study, we presented a novel feature selection method: nonnegative singular value approximation (NSVA) to enhance RNA-seq differential expression analysis by taking advantage of RNA-seq count data's non-negativity. As a variance-based feature selection method, it selects genes according to its contribution to the first singular value direction of input data in a data-driven approach. It demonstrates robustness to depth bias and gene length bias in feature selection in comparison with its five peer methods. Combining with state-of-the-art RNA-seq differential expression analysis, it contributes to enhancing differential expression analysis by lowering false discovery rates caused by the biases. Furthermore, we demonstrated the effectiveness of the proposed feature selection by proposing a data-driven differential expression analysis: NSVA-seq, besides conducting network marker discovery.

Network module-based model in the differential expression analysis for RNA-seq.

RNA-seq has emerged as a powerful technology for the detection of differential gene expression in the transcriptome. The commonly used statistical methods for RNA-seq differential expression analysis were designed for individual genes, which may detect too many irrelevant significantly genes or too few genes to interpret the phenotypic changes. Recently network module-based methods have been proposed as a powerful approach to analyze and interpret expression data in microarray and shotgun proteomics. But the module-based statistical model has not been adequately addressed for RNA-seq data. we proposed a network module-based generalized linear model for differential expression analysis of the count-based sequencing data from RNA-seq. The simulation studies demonstrated the effectiveness of the proposed model and the improvement of the statistical power for identifying the differentially expressed modules in comparison to the existing methods. We also applied our method to tissue datasets and identified 207 significantly differentially expressed kidney-active or liver-active modules. For liver cancer datasets, significantly differentially expressed modules, including Wnt signaling pathway and VEGF pathway, were found to be tightly associated with liver cancer. Besides, in comparison with the single gene-level analysis, our method could identify more significantly biological modules, which related to the liver cancer. The R package SeqMADE is available at https://cran.r-project.org/web/packages/SeqMADE/ . jing.li@sjtu.edu.cn. Supplementary data are available at Bioinformatics online.

Assessment of histone tail modifications and transcriptional profiling during colon cancer progression reveals a global decrease in H3K4me3 activity

During colon cancer, epigenetic alterations contribute to the dysregulation of major cellular functions and signaling pathways. Modifications in chromatin signatures such as H3K4me3 and H3K9ac, which are associated with transcriptionally active genes, can lead to genomic instability and perturb the expression of gene sets associated with oncogenic processes. In order to further elucidate early pre-tumorigenic epigenetic molecular events driving CRC, we integrated diverse, genome-wide, epigenetic inputs (by high throughput sequencing of RNA, H3K4me3, and H3K9ac) and compared differentially expressed transcripts (DE) and enriched regions (DER) in an in-vivo rat colon cancer progression model. Carcinogen (AOM) effects were detected genome-wide at the RNA (116 DE genes), K9ac (49 DERs including 24 genes) and K4me3 (7678 DERs including 3792 genes) level. RNA-seq differential expression and pathway analysis indicated that interferon-associated innate immune responses were impacted by AOM exposure. Despite extensive associations between K4me3 DERs and colon tumorigenesis (1210 genes were linked to colorectal carcinoma) including FOXO3, GNAI2, H2AFX, MSH2, NR3C1, PDCD4 and VEGFA, these changes were not reflected at the RNA gene expression level during early cancer progression. Collectively, our results indicate that carcinogen-induced changes in gene K4me3 DERs are harbingers of future transcriptional events, which drive malignant transformation of the colon.

Powerful differential expression analysis incorporating network topology for next-generation sequencing data.

RNA-seq has become the technology of choice for interrogating the transcriptome. However, most methods for RNA-seq differential expression (DE) analysis do not utilize prior knowledge of biological networks to detect DE genes. With the increased availability and quality of biological network databases, methods that can utilize this prior knowledge are needed and will offer biologists with a viable, more powerful alternative when analyzing RNA-seq data. We propose a three-state Markov Random Field (MRF) method that utilizes known biological pathways and interaction to improve sensitivity and specificity and therefore reducing false discovery rates (FDRs) when detecting differentially expressed genes from RNA-seq data. The method requires normalized count data (e.g. in Fragments or Reads Per Kilobase of transcript per Million mapped reads (FPKM/RPKM) format) as its input and it is implemented in an R package pathDESeq available from Github. Simulation studies demonstrate that our method outperforms the two-state MRF model for various sample sizes. Furthermore, for a comparable FDR, it has better sensitivity than DESeq, EBSeq, edgeR and NOISeq. The proposed method also picks more top Gene Ontology terms and KEGG pathways terms when applied to real dataset from colorectal cancer and hepatocellular carcinoma studies, respectively. Overall, these findings clearly highlight the power of our method relative to the existing methods that do not utilize prior knowledge of biological network. As an R package at https://github.com/MalathiSIDona/pathDESeq. install_github("MalathiSIDona/pathDESeq",build_vignettes = TRUE). After installation, type vignette("pathDESeq") to access the vignette. a.salim@latrobe.edu.au. Supplementary data are available at Bioinformatics online.

Sample size calculation while controlling false discovery rate for differential expression analysis with RNA-sequencing experiments.

BackgroundRNA-Sequencing (RNA-seq) experiments have been popularly applied to transcriptome studies in recent years. Such experiments are still relatively costly. As a result, RNA-seq experiments often employ a small number of replicates. Power analysis and sample size calculation are challenging in the context of differential expression analysis with RNA-seq data. One challenge is that there are no closed-form formulae to calculate power for the popularly applied tests for differential expression analysis. In addition, false discovery rate (FDR), instead of family-wise type I error rate, is controlled for the multiple testing error in RNA-seq data analysis. So far, there are very few proposals on sample size calculation for RNA-seq experiments.ResultsIn this paper, we propose a procedure for sample size calculation while controlling FDR for RNA-seq experimental design. Our procedure is based on the weighted linear model analysis facilitated by the voom method which has been shown to have competitive performance in terms of power and FDR control for RNA-seq differential expression analysis. We derive a method that approximates the average power across the differentially expressed genes, and then calculate the sample size to achieve a desired average power while controlling FDR. Simulation results demonstrate that the actual power of several popularly applied tests for differential expression is achieved and is close to the desired power for RNA-seq data with sample size calculated based on our method.ConclusionsOur proposed method provides an efficient algorithm to calculate sample size while controlling FDR for RNA-seq experimental design. We also provide an R package ssizeRNA that implements our proposed method and can be downloaded from the Comprehensive R Archive Network (http://cran.r-project.org).Electronic supplementary materialThe online version of this article (doi:10.1186/s12859-016-0994-9) contains supplementary material, which is available to authorized users.

The mRNA expression of soluble urokinase plasminogen activator surface receptor in human adipose tissue is positively correlated with body mass index.

suPLAUR is the transcript variant that encodes the soluble form of the urokinase plasminogen activator surface receptor (suPLAUR). This soluble protein has been shown to enhance leukocyte migration and adhesion, and its circulatory level is increased in inflammatory states. In this pilot study, we used RNA-Seq to examine the splicing pattern of PLAUR in omental adipose tissues from obese and lean individuals. Of the three transcript variants of the PLAUR gene, only the proportion of suPLAUR (transcript variant 2) increases in obesity. After removing the effects of gender and age, the expression of suPLAUR is positively correlated with body mass index. This observation was validated using RT-qPCR with an independent cohort of samples. Additionally, in our RNA-Seq differential expression analysis, we also observed, in obese adipose tissues, an up-regulation of genes encoding other proteins involved in the process of chemotaxis and leukocyte adhesion; of particular interest is the integrin beta 2 (ITGB2) that is known to interact with suPLAUR in leukocyte adhesion. These findings suggest an important role for suPLAUR in the recruitment of immune cells to obese adipose tissue, in the pathogenesis of obesity.

How are Bayesian and Non-Parametric Methods Doing a Great Job in RNA-Seq Differential Expression Analysis? : A Review

In a short history, RNA-seq data have established a revolutionary tool to directly decode various scenarios occurring on whole genome-wide expression profiles in regards with differential expression at gene, transcript, isoform, and exon specific quantification, genetic and genomic mutations, and etc. RNA-seq technique has been rapidly replacing arrays with seq-based platform experimental settings by revealing a couple of advantages such as identification of alternative splicing and allelic specific expression. The remarkable characteristics of high-throughput large-scale expression profile in RNA-seq are lied on expression levels of read counts, structure of correlated samples and genes, larger number of genes compared to sample size, different sampling rates, inevitable systematic RNA-seq biases, and etc. In this study, we will comprehensively review how robust Bayesian and non-parametric methods have a better performance than classical statistical approaches by explicitly incorporating such intrinsic RNA-seq specific features with flexible and more appropriate assumptions and distributions in practice.

Abstract 974: Targeting the tumor microenvironment with antimiRs that exploit oncomiR addiction in lymphoma

Abstract The purpose of this study was to develop a tumor-targeted therapeutic to antagonize aberrantly expressed microRNAs in lymphoma. MicroRNA dysfunction has been linked to the onset cancer. For example, in lymphoma miR-155 can affect carcinogenesis by acting as a hyperexpressed microRNA oncogene, or oncomiR. As ubiquitous regulators of gene expression, microRNAs are promising therapeutic targets. As such, numerous recent preclinical and clinical studies have utilized therapeutic oligonucleotides, known as antimiRs, to bind to and inhibit disease-associated microRNAs. However, the efficacy of systemically administered antimiRs is hindered by non-specific tissue targeting, which typically results in substantial liver uptake. We have developed an antimiR-peptide conjugate which preferentially targets the acidic tumor microenvironment. Intravital imaging and confocal tissue analysis show that this conjugate accumulates in tumors, while avoiding the liver, and has no observable systemic or cellular toxicity. The antimiR component of this conjugate comprises charge-neutral peptide nucleic acids (PNA). Due to their stability and high binding affinity, PNAs are effective antimiRs; however, they do not efficiently enter cells without a delivery vector. In addition to its tumor targeting benefits, we have previously shown that this peptide can deliver cargo into cells by translocating tethered molecules directly across lipid membranes. Via this non-endocytic mechanism, this peptide also promotes the effective delivery of PNA antimiRs into cultured B lymphocytes. Therefore, this conjugate both targets tumors and delivers antimiRs into cells via non-canonical pathways. Here we will detail the design and characterization of this conjugate, with specific emphasis on its ability to deliver antimiRs into cells via a non-endocytic route. In cultured cells we will demonstrate the efficacy of miR-155 inhibition using reporter and cell viability assays. For in vivo studies, we developed a metastatic mouse model of inducible diffuse large B cell lymphoma that is addicted to miR-155, such that lymphomagenesis is dependent on overexpression of the oncomiR, and withdrawal of miR-155 leads to cancer regression. Using this spontaneous miR-155-addicted lymphoma tumor model, we will discuss the biodistribution of the conjugate, as well as the effects of antimiR-155 therapy on lymphoid tumor growth and metastasis. In addition, we will outline how antagonizing miR-155 affects novel gene targets and biological pathways involved in oncomiR addiction, which were uncovered via RNA-seq differential gene expression analysis. OncomiRs have recently established a new paradigm for anti-cancer gene therapy. This work introduces a versatile and non-toxic tumor microenvironment-targeted system that exploits these lynchpin molecules as therapeutic tools. Citation Format: Christopher J. Cheng, Don M. Engelman, Mark Saltzman, Frank J. Slack. Targeting the tumor microenvironment with antimiRs that exploit oncomiR addiction in lymphoma. [abstract]. In: Proceedings of the 105th Annual Meeting of the American Association for Cancer Research; 2014 Apr 5-9; San Diego, CA. Philadelphia (PA): AACR; Cancer Res 2014;74(19 Suppl):Abstract nr 974. doi:10.1158/1538-7445.AM2014-974

RNA-Skim: a rapid method for RNA-Seq quantification at transcript level.

Motivation: RNA-Seq technique has been demonstrated as a revolutionary means for exploring transcriptome because it provides deep coverage and base pair-level resolution. RNA-Seq quantification is proven to be an efficient alternative to Microarray technique in gene expression study, and it is a critical component in RNA-Seq differential expression analysis. Most existing RNA-Seq quantification tools require the alignments of fragments to either a genome or a transcriptome, entailing a time-consuming and intricate alignment step. To improve the performance of RNA-Seq quantification, an alignment-free method, Sailfish, has been recently proposed to quantify transcript abundances using all k-mers in the transcriptome, demonstrating the feasibility of designing an efficient alignment-free method for transcriptome quantification. Even though Sailfish is substantially faster than alternative alignment-dependent methods such as Cufflinks, using all k-mers in the transcriptome quantification impedes the scalability of the method.Results: We propose a novel RNA-Seq quantification method, RNA-Skim, which partitions the transcriptome into disjoint transcript clusters based on sequence similarity, and introduces the notion of sig-mers, which are a special type of k-mers uniquely associated with each cluster. We demonstrate that the sig-mer counts within a cluster are sufficient for estimating transcript abundances with accuracy comparable with any state-of-the-art method. This enables RNA-Skim to perform transcript quantification on each cluster independently, reducing a complex optimization problem into smaller optimization tasks that can be run in parallel. As a result, RNA-Skim uses <4% of the k-mers and <10% of the CPU time required by Sailfish. It is able to finish transcriptome quantification in <10 min per sample by using just a single thread on a commodity computer, which represents >100 speedup over the state-of-the-art alignment-based methods, while delivering comparable or higher accuracy.Availability and implementation: The software is available at http://www.csbio.unc.edu/rs.Contact: weiwang@cs.ucla.eduSupplementary information: Supplementary data are available at Bioinformatics online.

MultiRankSeq: multiperspective approach for RNAseq differential expression analysis and quality control.

Background. After a decade of microarray technology dominating the field of high-throughput gene expression profiling, the introduction of RNAseq has revolutionized gene expression research. While RNAseq provides more abundant information than microarray, its analysis has proved considerably more complicated. To date, no consensus has been reached on the best approach for RNAseq-based differential expression analysis. Not surprisingly, different studies have drawn different conclusions as to the best approach to identify differentially expressed genes based upon their own criteria and scenarios considered. Furthermore, the lack of effective quality control may lead to misleading results interpretation and erroneous conclusions. To solve these aforementioned problems, we propose a simple yet safe and practical rank-sum approach for RNAseq-based differential gene expression analysis named MultiRankSeq. MultiRankSeq first performs quality control assessment. For data meeting the quality control criteria, MultiRankSeq compares the study groups using several of the most commonly applied analytical methods and combines their results to generate a new rank-sum interpretation. MultiRankSeq provides a unique analysis approach to RNAseq differential expression analysis. MultiRankSeq is written in R, and it is easily applicable. Detailed graphical and tabular analysis reports can be generated with a single command line.

Count-based differential expression analysis of RNA sequencing data using R and Bioconductor

RNA sequencing (RNA-seq) has been rapidly adopted for the profiling of transcriptomes in many areas of biology, including studies into gene regulation, development and disease. Of particular interest is the discovery of differentially expressed genes across different conditions (e.g., tissues, perturbations) while optionally adjusting for other systematic factors that affect the data-collection process. There are a number of subtle yet crucial aspects of these analyses, such as read counting, appropriate treatment of biological variability, quality control checks and appropriate setup of statistical modeling. Several variations have been presented in the literature, and there is a need for guidance on current best practices. This protocol presents a state-of-the-art computational and statistical RNA-seq differential expression analysis workflow largely based on the free open-source R language and Bioconductor software and, in particular, on two widely used tools, DESeq and edgeR. Hands-on time for typical small experiments (e.g., 4-10 samples) can be <1 h, with computation time <1 d using a standard desktop PC.

A Comparison of Methods for RNA-Seq Differential Expression Analysis and a New Empirical Bayes Approach.

Transcriptome-based biosensors are expected to have a large impact on the future of biotechnology. However, a central aspect of transcriptomics is differential expression analysis, where, currently, deep RNA sequencing (RNA-seq) has the potential to replace the microarray as the standard assay for RNA quantification. Our contributions here to RNA-seq differential expression analysis are two-fold. First, given the high cost of an RNA-seq run, biological replicates are rare, and therefore, information sharing across genes to obtain variance estimates is crucial. To handle such information sharing in a rigorous manner, we propose an hierarchical, empirical Bayes approach (R-EBSeq) that combines the Cufflinks model for generating relative transcript abundance measurements, known as FPKM (fragments per kilobase of transcript length per million mapped reads) with the EBArrays framework, which was previously developed for empirical Bayes analysis of microarray data. A desirable feature of R-EBSeq is easy-to-implement analysis of more than pairwise comparisons, as we illustrate with experimental data. Secondly, we develop the standard RNA-seq test data set, on the level of reads, where 79 transcripts are artificially differentially expressed and, therefore, explicitly known. This test data set allows us to compare the performance, in terms of the true discovery rate, of R-EBSeq to three other widely used RNAseq data analysis packages: Cuffdiff, DEseq and BaySeq. Our analysis indicates that DESeq identifies the first half of the differentially expressed transcripts well, but then is outperformed by Cuffdiff and R-EBSeq. Cuffdiff and R-EBSeq are the two top performers. Thus, R-EBSeq offers good performance, while allowing flexible and rigorous comparison of multiple biological conditions.