Comparison of RNA-Seq and microarray in transcriptome profiling of activated T cells

Abstract

To demonstrate the benefits of RNA-Seq over microarray in transcriptome profiling, both RNA-Seq and microarray analyses were performed on RNA samples from a human T cell activation experiment. In contrast to other reports, our analyses focused on the difference, rather than similarity, between RNA-Seq and microarray technologies in transcriptome profiling. A comparison of data sets derived from RNA-Seq and Affymetrix platforms using the same set of samples showed a high correlation between gene expression profiles generated by the two platforms. However, it also demonstrated that RNA-Seq was superior in detecting low abundance transcripts, differentiating biologically critical isoforms, and allowing the identification of genetic variants. RNA-Seq also demonstrated a broader dynamic range than microarray, which allowed for the detection of more differentially expressed genes with higher fold-change. Analysis of the two datasets also showed the benefit derived from avoidance of technical issues inherent to microarray probe performance such as cross-hybridization, non-specific hybridization and limited detection range of individual probes. Because RNA-Seq does not rely on a pre-designed complement sequence detection probe, it is devoid of issues associated with probe redundancy and annotation, which simplified interpretation of the data. Despite the superior benefits of RNA-Seq, microarrays are still the more common choice of researchers when conducting transcriptional profiling experiments. This is likely because RNA-Seq sequencing technology is new to most researchers, more expensive than microarray, data storage is more challenging and analysis is more complex. We expect that once these barriers are overcome, the RNA-Seq platform will become the predominant tool for transcriptome analysis.

Figure 1. Experimental design.

Human CCR6⁺ CD4 memory T cells were stimulated with anti-CD3/anti-CD28 coated beads under Th17 condition as described in Materials and Methods. RNA samples were prepared from cells collected at 0, 2, 4, 6, 24 and 72 hour post-stimulation. Gene expressions of these samples were studied with both Affymetrx microarray and RNA-Sequencing technologies.

Figure 2. Comparison of expression profiles of 18,306 common genes between both platforms.

Scatter plots show the averages (between biological duplicates) of log2 transformed expression values between two platforms, at each individual time point. The relationship between the expression profiles generated in both platforms is depicted as either a smoothing spline (black) or a linear regression line (red). The intercept (b) and the slope (m) of the linear regression line, and the correlation coefficient (r) are reported at the top-left corner in each plot corresponding to each time point. The plots show that the overall dynamic range of the 18,306 common genes generated by the two platforms is much broader in RNA-Seq (2.6×10⁵) than in microarray (3.6×10³). Similar dynamic ranges are displayed in both platforms for genes with relative expression level between 0.55 and 0.95. In each platform, the relative expression level of each gene was determined based on the average of log2 transformed expression values in all 12 samples.

Figure 3. Comparison of ANOVA results between both platforms.

Between-group (BG) and within-group (WG) variances are presented as the square root of the mean difference for each gene, and are plotted against the relative expression levels on both platforms, microarray (MA, panel A) and RNA-Seq (RS, panel B). The averages of variances (in square root, panel C) for between- and within-groups are plotted (panel C). The F-scores of 18,306 common genes are compared between two platforms (panel D), and the averages of F-scores (in square root, panel E) are also presented along with the relative expression levels by using smoothing spline for both platforms. The distributions of FDR-adjusted p-values, based on F-scores in both platforms are presented (panel F).

Figure 4. Comparing differential expression profiles of 18,306 common genes between microarray and RNA-Seq.

Scatter plots of log2 transformed ratios (vs. baseline at T = 0 h) between both platforms at selected time points (T = 2, 6, and 72 hour) show similar results are observed at T = 4 and 24 hour. Genes that are specifically differentially expressed in microarray or RNA-Seq are colored in red and green respectively, and genes that are differentially expressed in both platforms are colored in blue.

Figure 5. Detection of splicing variants with RNA-Seq approach.

Human RORγand RORγτ are the two isoforms of the RORC nuclear receptor generated from alternative splicing of the gene. RNA-Seq results showed that RORγ and RORγτ are encoded in the minus strand of the gene and their mRNA transcripts share most of the exons except for one or two exons at the 5′ end. RORγ mRNA utilizes specific exons 1 and 2 whereas RORγt mRNA has its specific exon 1, and the two transcripts are driven by their distinct promoters. In microarray studies, these two isoforms were indistinguishable since the two probe sets 228806_PM_at and 206419_PM_at hybridize with the exon regions that are common for these two isoforms. In contrast, RNA-Seq showed the specific expression of RORγt but not RORγ in CCR6⁺ CD4 T cells.

Figure 6. Comparison of RNA-Seq and Affymetrix microarray in detection of genes expressed at low levels.

Expression of MYCL1 was at low levels in CD4 T cells and the subtle change in MYCL1 expression in the process of T cell activation cannot be detected by microarray approach (center and right panels). The high sensitivity of RNA-Seq approach allows detection of a more than 32-fold decrease in MYCL1 mRNA expression upon T cell activation (left panel).

Figure 7. Comparison of RNA-Seq and Affymetrix microarray in detection of genes expressed at high levels.

β-actin (ACTB) is expressed at high levels in all conditions as measured by both RNA-Seq and microarray. Microarray detected similar levels of ACTB between resting and activated T cells whereas RNA-Seq showed a 2 to 4-fold increase in activated T cells. ACTB has been reported to be 5.3-fold up-regulated in in-vitro stimulated lymphocytes measured by quantitative polymerase chain reactions .

Figure 8. Detection of gene polymorphism with RNA-Seq approach.

A single nucleotide change was identified in the IL23R gene of this donor from sequence reads in RNA-Seq. The change results in a Gln to His mutation at the third amino acid of the N-terminal of IL23 receptor.

Figure 9. The controversy of redundant probe sets in microarray, and inconsistent results were obtained in Affymetrix microarray.

The bars in blue represent genes expressed at 0-axis indicates gene expression levels in log₂ scale.

Figure 10. Inaccurate annotation of probe set 224321_PM_at resulted in conflicting results between Affymetrix microarray and RNA-Seq approaches for TMEFF2 expression.

In Affymetrix microarray studies, probe set 224321_PM_at showed high expression of TMEFF2 at 0, 2, 4 and 6 hours. However, the high expression reported by probe set 224321_PM_at is supported neither by 233910_PM_at nor by 223557_PM_s_at. In contrast, there was no detectable expression of TMEFF2 in RNA-Seq studies. As a matter of fact, the probe set 224321_PM_at is mapped to a genomic region that is unrelated to TMEFF2, and thus this Affymetrix probe set is inaccurately annotated.

Figure 11. The gene expressions of EPB41 and TMEM200B are reported by microarray (top panel) and RNA-Seq (bottom panel).

Probe set 227386_PM_s_at targets the overlapping region of both gene EPB41 and TMEM200B (bottom panel). RNA-Seq clearly shows high expression of EPB41, but no expression of TMEM200B. The high expression of EPB41 is also shown by probe set 225051_PM_at (top panel). Because the Affymetrix library file associates probe set 227386_PM_s_at with only TMEM200B, it inaccurately reports high expression values for this gene.

Figure 12. The 6.3-fold decrease from 0 hour to 4 hour reported by probe set 205277_PM_at in microarray can only represent the change of two of the isoforms (NM_012231 and NM_001135610), but does not reflect the expression change for the entire PRDM2.

On the other hand, the reported 6.3-fold decrease was not supported by peer probe sets, either.