Volume 14 Supplement 9
Parallel comparison of Illumina RNA-Seq and Affymetrix microarray platforms on transcriptomic profiles generated from 5-aza-deoxy-cytidine treated HT-29 colon cancer cells and simulated datasets
- Xiao Xu†1,
- Yuanhao Zhang†3,
- Jennie Williams1,
- Eric Antoniou2,
- W Richard McCombie2,
- Song Wu3,
- Wei Zhu3,
- Nicholas O Davidson4,
- Paula Denoya1 and
- Ellen Li1Email author
© Xu et al.; licensee BioMed Central Ltd. 2013
Published: 28 June 2013
High throughput parallel sequencing, RNA-Seq, has recently emerged as an appealing alternative to microarray in identifying differentially expressed genes (DEG) between biological groups. However, there still exists considerable discrepancy on gene expression measurements and DEG results between the two platforms. The objective of this study was to compare parallel paired-end RNA-Seq and microarray data generated on 5-azadeoxy-cytidine (5-Aza) treated HT-29 colon cancer cells with an additional simulation study.
We first performed general correlation analysis comparing gene expression profiles on both platforms. An Errors-In-Variables (EIV) regression model was subsequently applied to assess proportional and fixed biases between the two technologies. Then several existing algorithms, designed for DEG identification in RNA-Seq and microarray data, were applied to compare the cross-platform overlaps with respect to DEG lists, which were further validated using qRT-PCR assays on selected genes. Functional analyses were subsequently conducted using Ingenuity Pathway Analysis (IPA).
Pearson and Spearman correlation coefficients between the RNA-Seq and microarray data each exceeded 0.80, with 66%~68% overlap of genes on both platforms. The EIV regression model indicated the existence of both fixed and proportional biases between the two platforms. The DESeq and baySeq algorithms (RNA-Seq) and the SAM and eBayes algorithms (microarray) achieved the highest cross-platform overlap rate in DEG results from both experimental and simulated datasets. DESeq method exhibited a better control on the false discovery rate than baySeq on the simulated dataset although it performed slightly inferior to baySeq in the sensitivity test. RNA-Seq and qRT-PCR, but not microarray data, confirmed the expected reversal of SPARC gene suppression after treating HT-29 cells with 5-Aza. Thirty-three IPA canonical pathways were identified by both microarray and RNA-Seq data, 152 pathways by RNA-Seq data only, and none by microarray data only.
These results suggest that RNA-Seq has advantages over microarray in identification of DEGs with the most consistent results generated from DESeq and SAM methods. The EIV regression model reveals both fixed and proportional biases between RNA-Seq and microarray. This may explain in part the lower cross-platform overlap in DEG lists compared to those in detectable genes.
In recent years, RNA-Seq emerged as an appealing alternative to classical microarrays in measuring global genomic expressions [1, 2]. The RNA-Seq technology has been applied to many human pathological studies such as prostate cancer , neurodegenerative disease , retina defection , and colorectal cancer . Gene detection in RNA-Seq, unlike microarray, is not dependent on probe design; rather it relies on short nucleotide reads mapping which can attain exceedingly high resolution. Furthermore, the RNA-Seq gene counts cover a larger dynamic range than microarray probe hybridization based design. On the other hand, microarray technology is still widely used because of lower costs and wider availability . Previous studies comparing parallel RNA-Seq with microarray data have reported good correlation between the two platforms [1, 8–13]. While classical correlation approaches can evaluate the strength of the association between the two platforms, they have been insufficient in gauging proportional and fixed biases between the two platforms. Given the uncertainties in measuring gene expressions for both platforms, we have therefore applied the Errors-In-Variables (EIV) regression model . The EIV model is a more suitable regression method for this type of platform comparison because (1) it reflects measurement errors from both platforms, (2) its goodness-of-fit measure reflects the Pearson correlation, yet with the added advantages of (3) providing a measure for fixed bias and, a measure for proportional bias .
A major rationale for conducting global transcriptomic studies is to identify genes that are differentially expressed between two or more biological conditions. In previous comparisons of the differentially expressed gene (DEG) lists generated using parallel RNA-Seq and microarray data, the biological groups that were studied were often very different (e.g. liver vs. kidney, or malignant breast cell line vs. normal breast cell line) [1, 8]. In the current study, parallel sets of RNA-Seq and Affymetrix microarray data were generated on a single HT-29 colon cancer cell line that was treated with and without 5-aza-deoxy-cytidine (5-Aza), a DNA methylation enzyme inhibitor. The concentrations of 5-Aza used in the present study (0 μM, 5 μM and 10 μM), approximated or exceeded the concentration previously reported to reverse hypermethylation of the SPARC (EMBL: ENSG00000113140) gene promoter and reverse suppression of SPARC mRNA expression in HT-29 cells . In this study, paired ends 100bp RNA-Seq data was generated as opposed to single end RNA-Seq data described in similar reports [1, 8, 10, 11, 13]. Moreover, most of the previous studies comparing the two platforms were usually based on one or two DEG detection methods, which were relatively outdated or not inclusive [7, 15, 17]. Our study surveyed an array of currently used algorithms to identify DEGs in parallel for both microarray and RNA-Seq data. We sought to determine which pair of microarray and RNA-Seq algorithms would yield the largest overlap in the DEG lists under the same statistical significance level. A simulation study was further conducted using published parallel RNA-Seq and microarray datasets , to assess the consistency of different DEG methods across platforms and their ability in identifying true positives. Quantitative reverse transcriptase polymerase chain reaction assays (qRT-PCR) was used to assay expression of the SPARC gene and other DEGs selected by using 1) both datasets, 2) RNASeq data only and 3) microarray data only. Finally we determined which Ingenuity Pathways Analysis (IPA) canonical pathways were identified by 1) both datasets, 2) RNASeq data only and 3) microarray data only.
5-Aza treatment of HT-29 cells
The HT-29 (ATCC) colon cancer cell line was maintained in DMEM supplemented with 10% fetal bovine serum, 1% kanamycin, streptomycin-penicillin, and incubated at 37° C and 5% CO2. Three replicative 150 mm cultures were treated with: 1) dimethyl sulfoxide (vehicle alone, 0 μM 5-Aza); 2) 5 μM 5-Aza and 3) 10 μM 5-Aza; for five days. These 5-Aza concentrations are similar and greater than the 5-Aza concentration previously reported to increase apoptosis, alter genome methylation as well as mRNA gene expression in HT-29 cells . The HT-29 cells were washed with phosphate buffered saline on the plate prior to scraping and centrifuging the cells. Total RNA was extracted separately from each of these nine cultures using TRI Reagent according to the manufacturer's recommendations. The RNA quality was assessed using an Agilent 2100 Bioanalyzer (Agilent Tech., Palo Alto, CA) to have a RNA Integrity Number score ≥ 7. Each of the nine RNA samples was used to generate parallel RNA-Seq, microarray and qRT-PCR data.
Illumina RNA-Seq data
Aliquots (1 μg) of nine RNA samples (triplicate samples for each of three experimental conditions), were subjected to paired-ends 100 bp Illumina sequencing. The RNA-Seq libraries were prepared and sequenced at Cold Spring Harbor Laboratories using the TruSeq RNA Sample Preparation Kit (Illumina Inc., San Diego, CA). In brief, mRNA was purified and fragmented, followed by cDNA synthesis with random hexamers. This product then underwent end repair, adapter ligation, and size selection using AMPure XP beads (Beckman Coulter Inc., Brea, CA) to isolate DNA templates of 320nt fragments and to remove excess adapters. The cDNA was PCR amplified. Each library was sequenced using Illumina 2000 sequencer (Illumina Inc., San Diego, CA) on 2 lanes of the flow cell. Between 41 and 88 million reads were generated for each of the RNA samples. The sequences were filtered using FASTX-Toolkit  to remove sequences with low Phred scores (~ first 3 nucleotides). The short reads fastq files were processed using Tophat (v2.0.1)  and mapped to the reference Ensembl human genome 19 using default settings for paired reads. Cufflink program (version 1.3.0)  or HTSeq-count (v 0.5.3p7)  were subsequently employed to convert aligned short reads (BAM format) into Fragments Per Kilobase of exon model per Million mapped fragments (FPKM) or raw gene counts. In the following step, a filter procedure was applied to remove gene entries with max alignment number of < 10 in all three replicates of the experimental groups (control or 0 μM, 5 μM and 10 μM 5-Aza). The RNA-Seq data were deposited in NCBI's Gene Expression Omnibus database with accession number GSE41588.
Affymetrix microarray data
Aliquots (150 ng) of the same nine RNA samples were each labeled (single color), hybridized to Affymetrix hgu133plus2.0 (Affymetrix Inc., Santa Clara, CA) arrays, and the arrays were scanned in the Stony Brook University DNA Microarray Core Facility, according to the manufacturer's protocol. Note each RNA sample was hybridized to a separate microarray chip. Microarray data were preprocessed using Bioconductor's affy package followed by a custom filter procedure to retain the probe entries that were present in all three biological triplicates of one experimental group (control or 0 μM, 5 μM and 10 μM 5-Aza). RMA normalization was applied to scale the replicates to a comparable range. If multiple probes on the array corresponded to a single gene, the probe with the highest intensities was used to represent the gene intensity. The microarray data have been deposited in NCBI's Gene Expression Omnibus and are accessible through GEO Series accession number GSE41364.
Platform comparison based on gene expression levels and correlations
General between-platform association analysis was applied to compare RNA-Seq with microarray data profiles. This includes a detectable gene determination for each group after the filter procedure, in which detectable genes were identified and compared respectively between the two platforms. In addition, the general gene expression profiles from RNA-Seq or microarray were examined in a scatter plot with Pearson and Spearman correlation coefficients calculated for all the genes (including those removed by the filtering procedure). Detectable genes which are RNA-Seq exclusive were compared to the overlapped ones using expression intensity histograms. This analysis was performed to verify the sensitivity of RNA-Seq technology in detecting genes expressed at low levels.
Errors-In-Variables (EIV) regression model
Both normalized microarray data and RNA-Seq FPKM values were transformed into log2 scale and subsequently converted to unit-free ratios by dividing a pre-selected housekeeping gene, ZNF311(EMBL: ENSG00000197935). This gene was selected based on its moderate intensity and consistent expression levels (rank of expression intensity) across all samples on both platforms. We did not use GAPDH (EMBL: ENSG00000111640) as housekeeping gene because it is highly expressed. In our experiment a moderate expressed gene is more suitable as the reference for all measured genes.
Here Yij denotes the normalized value of RNA-Seq expression for gene i and sample j and Xij represents the normalized microarray expression intensity. Moreover, is the expected value of Y; and are independent platform measurement errors with mean zero and variances and . A prerequisite of this EIV model is the homoscedasticity assumption and in practice we removed the top 1% of genes with the largest variation and examined the remaining genes using Levene's test  to ensure equal error variances on both platforms. The ratio of error variances λ is estimable when we have multiple observations from the same sample, which we fortunately do in this study with 3 replicates per sample. When the errors are normally distributed we can obtain the point estimators of the model parameters via the maximum likelihood method . The confidence intervals for the regression slope and intercept can be obtained via the bootstrap resampling method.
In our study, an EIV regression model was constructed for each of the three experimental HT-29 cell groups (control or 0 μM, 5 μM, and 10 μM 5-Aza) and the R rootSolve package (v. 1.6.3) was used to compute the point estimators for each regression model. The bootstrap resampling method with 1000 times resampling were performed to derive the corresponding 95% confidence interval for the regression intercept α and the regression slope ß as an estimate of the fixed and the proportional bias respectively. Statistically, the confidence interval of α covering 0 indicates an absence of fixed bias; whereas the confidence interval of ß encompassing 1 implies the absence of proportional bias.
DEG algorithms for microarray and RNA-Seq data
The T-test with Benjamini-Hochberg correction , SAM  and eBayes  algorithms were applied to the filtered Affymetrix microarray data to generate DEG lists ( > 2-fold, FDR ≤ 0.05) for the following two pair-wise comparisons: 1) 5µM vs. 0µM 5-Aza groups and 2) 10µM vs. 5µM 5-Aza groups, respectively. The Cuffdiff , SAMSeq , DESeq , baySeq  algorithms were applied to the filtered RNA-Seq data to generate DEG lists based on the same cutoff (> 2-fold or <0.5, FDR ≤ 0.05). NOISeq  was applied to the RNA-Seq data and the DEG list was subsequently filtered for a threshold of (> 2-fold or <0.5). The popular edgeR algorithm  was not included since it closely resembled the DESeq algorithm .
Comparing DEG algorithms using simulated data
In our simulation study, we designed a simulation method which generated consistent RNA-seq and microarray data in comparing DEG algorithms of the two platforms. The RNA-Seq and microarray simulations were built upon parallel RNA-Seq and microarray datasets reported previously by Marioni et.al. (GSE11045) . For this simulated analysis we could not apply the Cuffdiff algorithm because the previously published RNA-Seq data was reported only as raw gene counts without exon level information.
By averaging across the , we can approximately eliminate the effect of η and ε; and the transformed data could then be used to build an empirical distribution of u. The true expression levels of simulated genes were sampled from this empirical distribution in such a way that: a histogram of true data was generated using 500 bins at first step; a simulated gene was then assinged to a bin based on the frequency with a small turbulance added to its value. Uniform distribution (ranging the length of the bin) was assumed to the turbulance term to differentiate genes in the same group. As a result, the transformed expression level of a gene at a certain x% quantile of a given sample is equal to the same x% quantile of g(y) in terms of distribution + a small turbulance. The turbulance added same effect of variation to every gene because of the variance stablizing transformation.
The RNA-Seq data were fitted in a negative binomial model as described by Kvam . The mean expression level λ was sampled from a gamma distribution whose parameters were determined by fitting the true data with maximum likelihood method; similarly, the over-dispersion parameter φ was also generated from a gamma distribution described before . In practice, the distribution of λ was slightly rescaled to the range of the real data. Subsequently, we sampled both microarray g(y) and RNA-Seq λ of 10,000 genes from each corresponding distribution with a strict rule on quantile consistency (any gene of the α percentile of one distribution shall have the same quantile in the other dataset). In reality, sample percentile converge to distribution percentile when sample size is large, therefore a Spearman correlation of 1 was approximated in our simulated data sets across platforms.
Pre-defined significant DEGs were randomly sampled so that the log fold changes of these preset DEGs were generated from a mixed normal distribution where the probabilities of being up and down regulated were both equal. Moreover, in our analysis, the absolute expectation of log fold changes and standard errors for up and down regulated genes were set to be the same.
In practice, we generated 10 sets of simulated data in practice with 5 replicates included in both treatment and control group and 1000 random selected genes were preset to be differently expressed in different levels using a method of 95% minimum fold change such that the 1000 preset bona fida DEGs were generated with their fold changes following a log-normal distribution with 95% of the 1000 genes having their fold changes above the given level. We implemented a total of seven algorithms in our study, namely T test, SAM and eBayes on microarray data and baySeq, DESeq, SAMseq, NOISeq for RNA-Seq data.
In this work, the sensitivity and false discovery rates (FDR) were firstly evaluated for each DEG method under the 95% minimum fold change of 2 for preset DEGs and FDR cutoff of 0.05. For NOISeq method, a q = 0.8 (recommended by NOISeq author) criterion was used due to the absence of FDR control in this method. We further evaluated the sensitivity and false positive rate of each DEG algorithm by varying the differential significance levels of the preset 1000 genes using 95% minimum fold change method. Specifically, a range of values from 0.5 to 4 by an increment of 0.5 were used to generate the simulated DEGs.
qRT-PCR analysis of control HT-29 RNA samples vs. 5 μM 5-Aza treated HT-29 RNA samples
Reverse transcription was performed on 5 μg aliquots of each of the six RNA samples using the SuperScript III First-Strand Synthesis System for RT-PCR (Invitrogen Life Technologies, Grand Island, NY) according to the manufacturer's protocol. QRT-PCR was initially performed on serial dilutions of the cDNA for each Taqman assay kit using Taqman assay kits (Invitrogen Inc., Carlsbad, CA) in order to confirm that each of the assays were conducted in the linear range and the slopes of the threshold cycle Ct when plotted against the dilution were the same for all of the assays. Thirteen genes were selected by the majority vote of platform specific DEG detection methods and are categorized into 3 groups, which are: 1) commonly identified on RNA-Seq and microarray datasets; 2) RNA-Seq data only and 3) microarray data only. This list of genes, with the additional SPARC and GAPDH (EMBL: ENSG00000111640), and the corresponding commercial Taqman assays are listed in Additional file 1. The qRT-PCR assays were conducted in triplicate for each RNA sample. The ΔCt values (Ct for GAPDH - Ct for the test gene) were calculated for each RNA sample. The Student t-test was used to analyze whether there was a significant difference between the mean ΔCt for the control vs. the 5 μM 5-Aza treated HT-29 groups, with a threshold significance level of 0.05. The fold change in gene expression was calculated as 2-ΔΔCt (ΔΔCt = ΔCt of 5-Aza group - ΔCt of control group).
Ingenuity Pathway Analysis of microarray and RNA-Seq data
Based on the results of the simulation, we performed IPA analysis (Ingenuity® Inc, Redwood city, CA) on up-regulated DEGs (5 µM vs. 0 µM 5-Aza) and down-regulated DEGs respectively. 5 DEG lists were generated by the SAM, eBayes, Cuffdiff, DESeq and baySeq algorithms. Significantly enriched canonical pathways were selected based on the p value cutoff of 0.05 and included gene number > 3 .
General association between the two platforms
Applying EIV model for platform comparison
Comparison of DEG algorithms applied to experimental microarray and RNA-Seq HT-29 data
Three microarray DEG algorithms (T-test, SAM, eBayes) and five RNA-Seq algorithms (Cuffdiff, SAMSeq, DESeq, baySeq and NOISeq) were applied to the experimental HT-29 microarray and RNA data, respectively (See Additional file 2). The threshold was set at fold-change > 2 or less than 0.5 and a false discovery rate ≤ 0.05 for all of the eight algorithms except NOISeq. Since setting a fold change was not an option for NOISeq, we set a threshold of q = 0.8 and then subsequently filtered the selected genes with a threshold of fold-change > 2 or less than 0.5. Treatment of HT-29 cells with 5 μM 5-Aza (compared to the control HT-29 cells) resulted in up-regulation (↑) and down-regulation (↓) of genes. The T-test identified 392↑ 148↓, SAM identified 794↑ 256↓ and eBayes identified 782↑ 259↓ using the same microarray data (~13,000 detectable genes). Cuffdiff found 1149↑ 558↓, SAMSeq found 2262 ↑ 282↓, DESeq found 1840↑ 300↓, baySeq found 2013↑ 293↓, and NOISeq identified 673↑ 151↓ using the same RNA-Seq data (~17,000 detectable genes). All of the algorithms demonstrated an overall upregulation of gene expression after treatment of 5 μM 5-Aza. This is consistent with the concept that 5-Aza treatment reverses hypermethylation of gene promoters in HT-29 colon cancer cells and thus activates corresponding genes. However, activation of SPARC gene expression, which was previously reported after treatment of HT-29 cells with 4 μM 5-Aza , was observed in the RNA-Seq data only, and not in the microarray data.
The effect of increasing the concentration of 5-Aza from 5 μM to 10 μM 5-Aza was also analyzed using the eight algorithms and the same threshold parameters. The T-test identified 0↑ 2↓, SAM identified 13↑ 285↓ and eBayes identified 41↑ 278↓ using the same microarray data (~13,000 detectable genes). Cuffdiff detected 15↑ 485↓, SAMSeq detected 0↑ 626↓, DESeq detected 43↑ 389↓, baySeq detected 58↑ 424↓, and NOISeq detected 95↑ 123↓ using the same RNA-Seq data (~17,000 detectable genes). With the exception of the T-test and NOISeq, the remaining six algorithms detected an overall down-regulation in gene expression when the concentration of 5-Aza was increased from 5 μM to 10 μM. This could reflect toxic effects of 5-Aza at the higher 10 μM concentration.
Cross-platform overlap in DEG lists using RNA-Seq and microarray HT-29 data.
5 µM vs 0 µM
10 µM vs 5 µM
Comparison of DEG algorithms applied to simulated microarray and RNA-Seq data
Intra- and cross-platform comparison of DEG lists generated from simulated microarray and RNA-Seq data.
Confirmation of DEGs selected using both RNA-Seq and microarray data, RNA-Seq only and microarray only by qRT-PCR.
Biological function analysis of DEG lists generated by microarray and RNA-Seq data
Pathways commonly detected by SAM, eBayes, DESeq and baySeq
Human Embryonic Stem Cell Pluripotency
Bladder Cancer Signaling
LPS/IL-1 Mediated Inhibition of RXR Function
Activation of IRF by Cytosolic Pattern Recognition Receptors
Antigen Presentation Pathway
Factors Promoting Cardiogenesis in Vertebrates
Role of BRCA1 in DNA Damage Response
Role of CHK Proteins in Cell Cycle Checkpoint Control
Hepatic Fibrosis/Hepatic Stellate Cell Activation
Type I Diabetes Mellitus Signaling
Estrogen-mediated S-phase Entry
Hereditary Breast Cancer Signaling
Neuroprotective Role of THOP1 in Alzheimer's Disease
Caveolar-mediated Endocytosis Signaling
Graft-versus-Host Disease Signaling
Protein Ubiquitination Pathway
Oncostatin M Signaling
Autoimmune Thyroid Disease Signaling
Role of IL-17A in Arthritis
Cell Cycle Regulation by BTG Family Proteins
Aryl Hydrocarbon Receptor Signaling
Role of Osteoblasts, Osteoclasts and Chondrocytes in Rheumatoid Arthritis
In order to evaluate the performance of paired-end RNA-Seq data with a widely used commercial microarray platform, we chose to generate parallel datasets in a well-characterized experimental system, treatment of HT-29 colon cancer cells with 5-Aza, a DNA methyltransferase inhibitor [16, 37]. The 5-Aza concentrations were chosen to approximate and exceed the concentration previously reported to increase apoptosis and alter genome methylation as well as mRNA gene expression in HT-29 cells . Specifically reversal of hypermethylation of the SPARC promoter and reversal of suppression of SPARC gene expression were reported . The RNA-Seq technology is rapidly advancing, hence paired-end rather than single end RNA-Seq data were generated for this study.
We first examined the detection sensitivity for both platforms. RNA-Seq detected more genes than microarray, particularly among genes expressed at low levels. This observation is consistent with previous studies [11, 38]. The higher sensitivity of RNA-Seq can be attributed to its detection mechanism based on single-read/nucleotide resolution . The microarray gene quantification results largely depend on the accuracy of probe fluorescence scanning; background signal and other confounding factors (e.g. stains on array surface) may conceal the real genetic signal for a probe having a low abundance. In this perspective, the difference in detection mechanism confers a natural advantage to RNA-Seq comparing to microarray. The genomic ranges covered by both platforms also differ significantly. In addition, RNA-Seq detects all sequences that are expressed and basically surveys all the known genes provided by hg19 reference genome (N = 23,368), whereas microarray only examines genes based on the predesigned probe sets included on the array (N = 18,209). The correlation analysis confirmed strong general concordance on the gene expression measurements across platforms. Both Pearson and the Spearman correlation coefficients between the two technologies were found well above 0.8 with P values << 0.001 indicating the data were in comparable quality to previously reported parallel microarray and RNA-Seq datasets [1, 11, 40]. Furthermore, the EIV regression model was applied since the classical correlation based analysis is insufficient in gauging the quantitative concordance of the two platforms and the existence of random errors in both measurements rendered the traditional ordinary least regression method unsuitable in the current case. As per our study, the EIV regression revealed the existence of both fixed and proportional biases between the microarray and RNA-Seq platforms. We found that the fixed bias plays a minor part while the proportional bias is the major source of discrepancy between the two platforms (Figure 2). Basically, an estimated fixed bias at -0.24 on the log2 scale reflected a trivial baseline difference, whereas an estimated ~1.45 proportional bias meant that a unit change on microarray gene intensity on the log2 scale corresponded to about 1.45 units change for RNA-Seq on the log2 scale. This regression model is consistent with the observation that RNA-Seq was more sensitive and exhibited a larger dynamic range than its microarray counterparts in measuring the expression level of the same transcript.
Since the major goal of conducting global transcriptomic studies is to identify genes that are differentially expressed between two or more biological groups, this study applied several DEG algorithms designed for either microarray or RNA-Seq data. Two of the most widely used microarray DEG algorithms in recent years, SAM and eBayes, are included in this study. The classical T-test, which is known to perform relatively poorly in microarray analysis was also evaluated as a "control" method . While microarray data produces a continuous intensity, which commonly follows a log-normal distribution , the RNA-Seq gene expression level is discrete or digital in nature. The microarray DEG algorithms are based on continuous distribution of random variables (after log transformation of the probe hybridization intensities). On the other hand, RNA-Seq DEG algorithms are rapidly evolving. The earlier studies mostly relied on algorithms assuming a Poisson distribution on the gene counts [1, 13, 39] while the more recent methods utilized a negative binomial model which was considered better than Poisson assumption in explaining biological variability of the RNA-Seq data [28, 43]. This study considers several of the currently used, popular RNA-Seq DEG algorithms: Cuffdiff, baySeq and DESeq which are roughly based on the negative binomial modeling of RNA-Seq data and the nonparametric SAMSeq and NOISeq methods, which are relatively model-free. Each of the methods has its own virtue and relevance: the Cuffdiff method is built to incorporate biological variability information (e.g. isoforms and fragment assignment uncertainty) from the initial short reads input. In baySeq algorithm, the estimate of significance is based on an empirical Bayes approach, which ranks the DEGs by posterior probabilities of the treatment group. DESeq assumes a locally linear relationship between variance and mean expression level. The SAMSeq algorithm, on the other hand, differs from the aforementioned algorithms by identifying DEGs using a Wilcoxon rank based nonparametric approach, which is relatively free from model biases. Lastly, the NOISeq algorithm evaluates the log-ratio of normalized counts (M value) versus their absolute difference (D value) and determined their differential significance by comparing to the noise distribution, and is designed to overcome the sequencing depth dependency commonly seen in other DEG methods.
Our simulation experiment using preset, true-positive genes at a minimal fold change of 2, demonstrated maximal cross-platform overlaps in the DEG lists generated by two of the RNA-Seq algorithms, baySeq and DESeq, and by two microarray methods, eBayes and SAM (Table 2). These observations are consistent with our results obtained using the HT-29 experimental data. Note however, that we were not able to evaluate the Cuffdiff algorithm using the simulated dataset. When the sensitivity of all the DEG methods were also examined in our study, the results showed that baySeq performed best among all RNA-Seq algorithms evaluated, in identifying true positive genes at each 95% minimal fold-change level. This observation is consistent with a previous study in which baySeq was found superior in ranking genes by significance to be declared . DESeq tails immediately after baySeq in sensitivity curves and performed comparably well at lower fold change levels (e.g. log2 fold change ~ 1.5). The microarray DEG algorithms, SAM and eBayes, were generally found less sensitive than RNA-Seq programs.
With respect to FDR evaluation, however, baySeq resulted in more false positive calls than most of the other RNA-Seq algorithms except for NOISeq, especially when the 95% minimum fold changes of true positive genes are higher (Figure 3B, right section). DESeq constantly results in the lowest FDR among all the RNA-Seq algorithms evaluated in the simulation experiments, indicating its superior reliability. The NOISeq showed a very poor performance on FDR evaluation curve particularly with lower 95% minimal fold change thresholds (Figure 3B, left section), reflecting the fact that NOISeq's DEG discerning power by comparing noise distribution against a true signal was seriously compromised when the 'true difference' is less remarkable. In practice, it is of theoretical importance to weigh more on preventing false positives than false negatives; we thus favor DESeq over Bayseq in RNA-Seq analysis as the former method controls FDR better than the latter in higher differential significance level (Figure 3B, right section).
Of the two microarray DEG algorithms, SAM slightly outperforms Ebayes in both sensitivity and FDR evaluation. The traditional T-test with BH correction, not surprisingly, showed a very poor performance in identifying true positives, probably due to its inappropriate independence assumption. When we view our results from the perspective of platform comparison, it is generally expected that DESeq and SAM can lead to consistent and reasonable DEG results -- an observation which is exactly reflected in our HT-29 experiment (Table 1).
Finally, to begin to address the biological significance of these studies, we undertook to validate that treatment of HT-29 colon cancer cells with 5 μM 5-Aza would relieve suppression of SPARC gene expression. While this anticipated outcome was confirmed using both the RNA-Seq data and qRT-PCR data, it was not observed in the microarray data. In addition a higher percentage of other DEGs identified using both platforms or RNA-Seq only was confirmed by qRT-PCR than the DEGs identified using microarray alone.
A strong correlation of genomic expression profiles was observed between the microarray and RNA-Seq platforms with the latter technology detecting more genes across the genome. Remarkable differences between the two platforms in terms of (1) the existence of both fixed and proportional biases detected by the errors-in-variable (EIV) regression model, and (2) discrepancies in DEG identification have been discovered in our study. We further confirmed that the DEG discrepancies are mostly related to the different algorithms used for both platforms. Among all the DEG algorithms surveyed in this study, the largest cross-platform overlaps were observed between the DEG lists generated by two RNA-Seq algorithms, baySeq and DESeq, and the DEG lists generated by two microarray algorithms, SAM and eBayes, from the HT-29 experimental dataset. The simulation studies, which did not include evaluation of Cuffdiff, indicate that the the DESeq algorithm outperformed the other RNA-Seq algorithms, based upon the combined considerations of sensitivity and false discovery rate. DESeq also demonstrated the highest overlap rate with the DEG list generated by SAM from the microarray data. Overall, the nonparametric based DEG methods such as SAMSeq or NOISeq exhibited suboptimal performance compared to their parametric counterparts, partly due to the limited number of replicates. QRT-PCR validated a higher percentage of the DEGs identified by both platforms and RNA-Seq only, than the DEGs identified by microarray only. Finally, while there were common IPA canonical pathways identified by both microarray and RNA-Seq data, a large number of additional canonical pathways were identified by RNA-Seq data alone. No additional canonical pathways were identified by microarray data alone.
List of abbreviations used
differentially expressed genes
false discovery rate
quantitative reverse transcriptase polymerase chain reaction.
The authors would like to thank Molly Hammel from Cold Spring Harbor Laboratory for her valuable suggestions in the RNA-Seq data analysis process. The authors would also like to express their sincere gratitude to Maxine McGredy for her valuable assistance in designing and performing the qRT-PCR experiments. Finally, the authors are thankful to all lab technicians from both Stony Brook University Health Science Center and Cold Spring Harbor Laboratory involved in experiments relevant to this study.
Publication of this article is funded by the Stony Brook Institute for Clinical and Translational Sciences Fusion Award (P.D.), the Simons Foundation (E.L.), and NIH grants R01CA140487 (J.W.), NIH RO1DK-56260 (N.D.), RO1HL-38180 (N.D.) and P30DK-52574 (N.D.). RNA-Sequencing was performed at the Cold Spring Harbor Laboratory Genome Sequencing Core, which is partially supported by NCI center grant CA045508.
This article has been published as part of BMC Bioinformatics Volume 14 Supplement 9, 2013: Selected articles from the 8th International Symposium on Bioinformatics Research and Applications (ISBRA'12). The full contents of the supplement are available online at http://www.biomedcentral.com/bmcbioinformatics/supplements/14/S9.
- Marioni JC, Mason CE, Mane SM, Stephens M, Gilad Y: RNA-seq: an assessment of technical reproducibility and comparison with gene expression arrays. Genome research. 2008, 18 (9): 1509-1517. 10.1101/gr.079558.108.PubMed CentralView ArticlePubMedGoogle Scholar
- Oshlack A, Robinson MD, Young MD: From RNA-seq reads to differential expression results. Genome biology. 2010, 11 (12): 220-10.1186/gb-2010-11-12-220.PubMed CentralView ArticlePubMedGoogle Scholar
- Ren S, Peng Z, Mao JH, Yu Y, Yin C, Gao X, Cui Z, Zhang J, Yi K, Xu W: RNA-seq analysis of prostate cancer in the Chinese population identifies recurrent gene fusions, cancer-associated long noncoding RNAs and aberrant alternative splicings. Cell research. 2012, 22 (5): 806-821. 10.1038/cr.2012.30.PubMed CentralView ArticlePubMedGoogle Scholar
- Courtney E, Kornfeld S, Janitz K, Janitz M: Transcriptome profiling in neurodegenerative disease. Journal of neuroscience methods. 2010, 193 (2): 189-202. 10.1016/j.jneumeth.2010.08.018.View ArticlePubMedGoogle Scholar
- Farkas MH, Grant GR, Pierce EA: Transcriptome analyses to investigate the pathogenesis of RNA splicing factor retinitis pigmentosa. Advances in experimental medicine and biology. 2012, 723: 519-525. 10.1007/978-1-4614-0631-0_65.View ArticlePubMedGoogle Scholar
- Castellarin M, Warren RL, Freeman JD, Dreolini L, Krzywinski M, Strauss J, Barnes R, Watson P, Allen-Vercoe E, Moore RA: Fusobacterium nucleatum infection is prevalent in human colorectal carcinoma. Genome research. 2012, 22 (2): 299-306. 10.1101/gr.126516.111.PubMed CentralView ArticlePubMedGoogle Scholar
- Fu X, Fu N, Guo S, Yan Z, Xu Y, Hu H, Menzel C, Chen W, Li Y, Zeng R: Estimating accuracy of RNA-Seq and microarrays with proteomics. BMC genomics. 2009, 10: 161-10.1186/1471-2164-10-161.PubMed CentralView ArticlePubMedGoogle Scholar
- Bradford JR, Hey Y, Yates T, Li Y, Pepper SD, Miller CJ: A comparison of massively parallel nucleotide sequencing with oligonucleotide microarrays for global transcription profiling. BMC genomics. 2010, 11: 282-10.1186/1471-2164-11-282.PubMed CentralView ArticlePubMedGoogle Scholar
- Mokry M, Hatzis P, Schuijers J, Lansu N, Ruzius FP, Clevers H, Cuppen E: Integrated genome-wide analysis of transcription factor occupancy, RNA polymerase II binding and steady-state RNA levels identify differentially regulated functional gene classes. Nucleic acids research. 2012, 40 (1): 148-158. 10.1093/nar/gkr720.PubMed CentralView ArticlePubMedGoogle Scholar
- Bottomly D, Walter NA, Hunter JE, Darakjian P, Kawane S, Buck KJ, Searles RP, Mooney M, McWeeney SK, Hitzemann R: Evaluating gene expression in C57BL/6J and DBA/2J mouse striatum using RNA-Seq and microarrays. PloS one. 2011, 6 (3): e17820-10.1371/journal.pone.0017820.PubMed CentralView ArticlePubMedGoogle Scholar
- Su Z, Li Z, Chen T, Li QZ, Fang H, Ding D, Ge W, Ning B, Hong H, Perkins RG: Comparing next-generation sequencing and microarray technologies in a toxicological study of the effects of aristolochic acid on rat kidneys. Chemical research in toxicology. 2011, 24 (9): 1486-1493. 10.1021/tx200103b.View ArticlePubMedGoogle Scholar
- Lahiry P, Lee LJ, Frey BJ, Rupar CA, Siu VM, Blencowe BJ, Hegele RA: Transcriptional profiling of endocrine cerebro-osteodysplasia using microarray and next-generation sequencing. PloS one. 2011, 6 (9): e25400-10.1371/journal.pone.0025400.PubMed CentralView ArticlePubMedGoogle Scholar
- Liu S, Lin L, Jiang P, Wang D, Xing Y: A comparison of RNA-Seq and high-density exon array for detecting differential gene expression between closely related species. Nucleic acids research. 2011, 39 (2): 578-588. 10.1093/nar/gkq817.PubMed CentralView ArticlePubMedGoogle Scholar
- Lancaste T: A Note on an Errors in Variables Model. J Am Stat Assoc. 1966, 61 (313): 128-&.Google Scholar
- Linnet K: Evaluation of Regression Procedures for Methods Comparison Studies. Clin Chem. 1993, 39 (3): 424-432.PubMedGoogle Scholar
- Cheetham S, Tang MJ, Mesak F, Kennecke H, Owen D, Tai IT: SPARC promoter hypermethylation in colorectal cancers can be reversed by 5-Aza-2'deoxycytidine to increase SPARC expression and improve therapy response. British journal of cancer. 2008, 98 (11): 1810-1819. 10.1038/sj.bjc.6604377.PubMed CentralView ArticlePubMedGoogle Scholar
- Bolstad BM, Irizarry RA, Astrand M, Speed TP: A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics. 2003, 19 (2): 185-193. 10.1093/bioinformatics/19.2.185.View ArticlePubMedGoogle Scholar
- FASTX-Toolkit. [http://hannonlab.cshl.edu/fastx_toolkit/]
- Trapnell C, Pachter L, Salzberg SL: TopHat: discovering splice junctions with RNA-Seq. Bioinformatics. 2009, 25 (9): 1105-1111. 10.1093/bioinformatics/btp120.PubMed CentralView ArticlePubMedGoogle Scholar
- Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, Salzberg SL, Wold BJ, Pachter L: Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nature biotechnology. 2010, 28 (5): 511-515. 10.1038/nbt.1621.PubMed CentralView ArticlePubMedGoogle Scholar
- HTSeq-count. [http://www-huber.embl.de/users/anders/HTSeq/doc/count.html]
- Levene H, Olkin II, Hotelling H: Robust tests for equality of variances. 1960, Stanford University PressGoogle Scholar
- Barnett VD: Fitting Straight Lines-The Linear Functional Relationship with Replicated Observations. Journal of the Royal Statistical Society Series C (Applied Statistics). 1970, 135-144.Google Scholar
- Benjamini Y, Hochberg Y: Controlling the False Discovery Rate - a Practical and Powerful Approach to Multiple Testing. J Roy Stat Soc B Met. 1995, 57 (1): 289-300.Google Scholar
- Tusher VG, Tibshirani R, Chu G: Significance analysis of microarrays applied to the ionizing radiation response (vol 98, pg 5116, 2001). Proceedings of the National Academy of Sciences of the United States of America. 2001, 98 (18): 10515-10515.View ArticleGoogle Scholar
- Smyth GK: Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Statistical applications in genetics and molecular biology. 2004, 3: Article3-View ArticlePubMedGoogle Scholar
- Li J, Tibshirani R: Finding consistent patterns: A nonparametric approach for identifying differential expression in RNA-Seq data. Statistical methods in medical research. 2011Google Scholar
- Anders S, Huber W: Differential expression analysis for sequence count data. Genome biology. 2010, 11 (10): R106-10.1186/gb-2010-11-10-r106.PubMed CentralView ArticlePubMedGoogle Scholar
- Hardcastle TJ, Kelly KA: baySeq: empirical Bayesian methods for identifying differential expression in sequence count data. BMC bioinformatics. 2010, 11: 422-10.1186/1471-2105-11-422.PubMed CentralView ArticlePubMedGoogle Scholar
- Tarazona S, Garcia-Alcalde F, Dopazo J, Ferrer A, Conesa A: Differential expression in RNA-seq: a matter of depth. Genome research. 2011, 21 (12): 2213-2223. 10.1101/gr.124321.111.PubMed CentralView ArticlePubMedGoogle Scholar
- Robinson MD, McCarthy DJ, Smyth GK: edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010, 26 (1): 139-140. 10.1093/bioinformatics/btp616.PubMed CentralView ArticlePubMedGoogle Scholar
- Rocke DM, Durbin B: A model for measurement error for gene expression arrays. Journal of computational biology: a journal of computational molecular cell biology. 2001, 8 (6): 557-569. 10.1089/106652701753307485.View ArticleGoogle Scholar
- Durbin BP, Hardin JS, Hawkins DM, Rocke DM: A variance-stabilizing transformation for gene-expression microarray data. Bioinformatics. 2002, 18 (Suppl 1): S105-110. 10.1093/bioinformatics/18.suppl_1.S105.View ArticlePubMedGoogle Scholar
- RMAExpress. [http://rmaexpress.bmbolstad.com/]
- Kvam VM, Liu P, Si Y: A comparison of statistical methods for detecting differentially expressed genes from RNA-seq data. American journal of botany. 2012, 99 (2): 248-256. 10.3732/ajb.1100340.View ArticlePubMedGoogle Scholar
- Zhang T, DeSimone RA, Jiao X, Rohlf FJ, Zhu W, Gong QQ, Hunt SR, Dassopoulos T, Newberry RD, Sodergren E: Host genes related to paneth cells and xenobiotic metabolism are associated with shifts in human ileum-associated microbial composition. PloS one. 2012, 7 (6): e30044-10.1371/journal.pone.0030044.PubMed CentralView ArticlePubMedGoogle Scholar
- Karpf AR, Peterson PW, Rawlins JT, Dalley BK, Yang Q, Albertsen H, Jones DA: Inhibition of DNA methyltransferase stimulates the expression of signal transducer and activator of transcription 1, 2, and 3 genes in colon tumor cells. Proceedings of the National Academy of Sciences of the United States of America. 1999, 96 (24): 14007-14012. 10.1073/pnas.96.24.14007.PubMed CentralView ArticlePubMedGoogle Scholar
- Sultan M, Schulz MH, Richard H, Magen A, Klingenhoff A, Scherf M, Seifert M, Borodina T, Soldatov A, Parkhomchuk D: A global view of gene activity and alternative splicing by deep sequencing of the human transcriptome. Science. 2008, 321 (5891): 956-960. 10.1126/science.1160342.View ArticlePubMedGoogle Scholar
- Wilhelm BT, Marguerat S, Watt S, Schubert F, Wood V, Goodhead I, Penkett CJ, Rogers J, Bahler J: Dynamic repertoire of a eukaryotic transcriptome surveyed at single-nucleotide resolution. Nature. 2008, 453 (7199): 1239-U1239. 10.1038/nature07002.View ArticlePubMedGoogle Scholar
- Davidson RM, Hansey CN, Gowda M, Childs KL, Lin H, Vaillancourt B, Sekhon RS, Natalia de Leon, Kaeppler SM, Jiang N: Utility of RNA Sequencing for Analysis of Maize Reproductive Transcriptomes. The Plant Genome. 2011, 4: 191-203. 10.3835/plantgenome2011.05.0015.View ArticleGoogle Scholar
- Jeanmougin M, de Reynies A, Marisa L, Paccard C, Nuel G, Guedj M: Should We Abandon the t-Test in the Analysis of Gene Expression Microarray Data: A Comparison of Variance Modeling Strategies. PloS one. 2010, 5 (9):
- Hoyle DC, Rattray M, Jupp R, Brass A: Making sense of microarray data distributions. Bioinformatics. 2002, 18 (4): 576-584. 10.1093/bioinformatics/18.4.576.View ArticlePubMedGoogle Scholar
- Robinson MD, Smyth GK: Moderated statistical tests for assessing differences in tag abundance. Bioinformatics. 2007, 23 (21): 2881-2887. 10.1093/bioinformatics/btm453.View ArticlePubMedGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.