Comparison of Gene Expression Data Obtained from MPSS and Microarrays

We believe that there should be some biologically meaningful relationships among different samples from different technology platforms. For example, samples from same tissues with different condition/treatment are expected to have higher correlation than samples from different tissues, and the same samples measured by different technology platforms would be expected to have higher correlation than others. To test our hypothesis, we make comparisons of abundance expression level and comparisons of expression level ratios.

Comparison of Abundance Expression Levels

We wanted to compare the abundance of transcripts as measured by the expression data from MPSS, Affymetrix and Agilent technologies (Table 1). It is important to note that not all genes were measurable with all three technologies (Figure 1), therefore our analyses focused on the subset of genes that were represented on both microarray platforms and contained a Dpn II site detectable by MPSS (16,269 genes). The MPSS data showed a wider dynamic range than Affymetrix and Agilent data (Table 2). The standard deviations for MPSS data were greater than those for Affymetrix and Agilent data.

#	Tissue/Treament	MPSS	Affymetrix	Agilent
1	Callus	MCas	AFCaA/AFCaB	AGCaA/AGCaB
2	Inflorescence	MIns	AFInA/AFInB	AGInA/AGInB, AGSiInA/AGSiInB
3	Leaf	MLes	AFLeA/AFLeB	AGLeA/AGLeB, AGRoLeA/AGRoLeB
4	Root	MRos	AFRoA/AFRoB	AGRoA/AGRoB, AGRoLeA/AGRoLeB
5	Silique	MSis	AFSiA	AGSiInA/AGSiInB
6	Agamous Inflorescence.	MAgm	AFAgmA/AFAgmB
7	Ap1-10 Inflorescence.	MAp1	AFAp1A/AFAp1B
8	Ap3-6 Inflorescence.	MAp3	AFAp3A/AFAp3B
9	Sup/Ap1 Inflorescence.	MSap	AFSapA/AFSapB
10	Leaves SA 4 hr.	MS04	AFS04A/AFS04B
11	Leaves SA 52 hr.	MS52	AFS52A/AFS52B

The eleven RNA sources and corresponding platform or sample codes used in this paper.

Code names:

first two letters of the code: M-MPSS; AF-Affymetrix; AG-Agilent

middle letters of the code: Tissue/Treatment

last letter of the code: A/B-technical replicates

Code	Minimum	Maximum	1st Quartile	Median	3rd Quartile	Std
MCas	0	13615	0	39.74	138.13	460.66
MIns	0	13783	9.69	36.34	109.03	514.51
MLes	0	30874	7.38	41.84	120.59	699.07
MRos	0	17858	7.68	42.23	124.76	498.47
MSis	0	10614	0	44.54	152.37	387.94
AFCa	4.88	2452.5	43.69	82.95	166.42	214.09
AFIn	4.99	2422.8	51.88	88.06	162.84	206.16
AFLe	6.91	3056.1	43.3	78.11	155.42	249.48
AFRo	3.86	3238.8	37.9	69.95	147.11	271.89
AFSi	5.44	2421	45.92	81.65	158.98	225.03
AGCa	0.47	4229.7	36.48	74.02	150.55	304.25
AGIn	1.01	3345.6	44.8	83.74	162.09	242.15
AGLe	0.75	5092.1	39.17	76.78	160	285.71
AGRo	1.41	3375.4	39.38	77.51	166.92	240.88

Std represents the standard error of the mean.

To make the data directly comparable, we rescaled the normalized data into Normalized Expression Unit (NEU) values (see Data Analysis in Methods), which we used for abundance expression levels comparison. We calculated the Pearson's correlation and Spearman's rank correlation of logged NEU normalized expression data for 1,648 genes in MPSS, Affymetrix and Agilent (i.e., the genes that were selected to have a MPSS tag count greater than 4 TPM (Transcripts Per Million) and "present" calls determined by the Affymetrix software), and we performed hierarchical clustering using these data (Figure 2). The NEU normalization did not affect the correlation statistics for the logged expression measurements since it corresponds to a shift in the mean for logged data. The within-technology correlations across the tissues and mutants were greater than the inter-technology correlations of identical RNA samples (Table 3). The patterns of clustering and correlation among the samples were similar within technologies, although not identical, when based on either MPSS, Affymetrix or Agilent measurements. Different tissue samples in within-technology comparisons exhibited similar patterns: root was mostly correlated with callus; silique was primarily correlated with flower and secondarily correlated with leaf. The treated leaf samples were clustered together, as were the mutant flower samples.

	MCas	MIns	MLes	MRos	AFCa	AFIn	AFLe	AFRo	AGCa	AGIn	AGLe	AGRo
MCas		0.46	0.41	0.64	0.66	0.35	0.29	0.50	0.54	0.26	0.19	0.36
MIns			0.62	0.56	0.44	0.71	0.58	0.50	0.34	0.56	0.42	0.39
MLes				0.50	0.40	0.55	0.70	0.43	0.32	0.44	0.53	0.35
MRos					0.63	0.51	0.44	0.76	0.50	0.39	0.31	0.55
AFCa						0.68	0.60	0.83	0.74	0.45	0.35	0.54
AFIn							0.85	0.77	0.44	0.69	0.53	0.51
AFLe								0.69	0.4	0.60	0.67	0.48
AFRo									0.59	0.52	0.41	0.65
AGCa										0.76	0.67	0.82
AGIn											0.85	0.82
AGLe												0.80
AGRo

The upper triangle represents the Pearson's correlation while the corresponding P-values are all smaller than 0.01. The correlation values greater than 0.8 are shown in bold; values between 0.6 to 0.8 are shown in italics. All genes had measurable abundance values greater than 4 TPM in all MPSS samples and were detected as present in Affymetrix.

When comparing the same RNA samples across technologies, the highest correlations were observed for MPSS versus Affymetrix microarrays (0.66 to 0.76) and for Affymetrix microarrays versus Agilent microarrays (0.65 to 0.74) while the lowest correlations were exhibited for Agilent microarrays versus MPSS (0.53 to 0.56) (Table 3). The correlations among different tissues within the MPSS data were lower than the correlations among different tissues measured by Affymetrix and Agilent microarrays. The correlation values among MPSS expression data ranged from 0.41 to 0.64 while the same measurements using Affymetrix or Agilent microarrays ranged from 0.60 to 0.85 and 0.67 to 0.82, respectively. Additional analyses using different Affymetrix and Agilent technical replicates also showed that each pair of technical replicates were closely correlated together, as expected (see Figure 1 in additional file 1, e.g, AFRoA and AFRoB).

We next assessed the NEU levels from MPSS and Affymetrix technologies in a gene-by-gene comparison of 6,501 genes, which were selected because their mean MPSS abundance across eleven samples was greater than 4 TPM and they were detected as present in the corresponding Affymetrix microarrays for the eleven RNA samples. There were numerous transcripts for which abundance was measurable with the Affymetrix microarrays, but for which no or very little expression was measureable in the MPSS analysis (Figure 3).

Comparison of Expression Level Ratios

In our comparison of the MPSS, Affymetrix and Agilent normalized abundance data, we observed that the RNA samples exhibited higher correlations within the technology platform used to measure RNA abundance rather than expected similarities due to the biological nature of the samples. Although MPSS and, in principle, Affymetrix data can be used to estimate abundance expression levels it is likely that probe-specific effects and other technological biases will have a large effect on the estimates. In practice, practitioners often focus on differences in expression between samples. Thus, we decided to use expression ratios of pairwise sample comparisons to reassess the data correlations. Nine samples were selected for this analysis: callus, inflorescence, leaf, root, silique, ap1, ap3, s04 and s52. We used the following eight pairs of samples for calculation of ratios: silique/inflorescence, silique/ap1, silique/ap3, ap1/ap3, callus/leaf, root/leaf, leaf/s04 and leaf/s52. For each technology platform, we computed the Pearson's correlations and Spearman's rank correlations between the ratios across pairs of samples. The hierarchical clustering of the correlations (Figure 4) showed that the samples were more closely clustered according to their biological nature than the technology platform. Therefore, the ratios were a more biologically meaningful measure to use for cluster analysis than abundance measurements of gene expression.

We next used a contour matrix plot [25] on expression ratios to assess the correlations among the different technology platforms (Figure 5). The Affymetrix and Agilent data were most similar to each other (e.g., AFRoLe and AGRoLe), indicating that the two microarray platforms were more consistent with each other in terms of detection of differential expression than with MPSS. The correlation values using the same tissue pairs for Affymetrix/Agilent and MPSS ranged from 0.38 to 0.71 while the values for Agilent and Affymetrix ranged from 0.72 to 0.91 (Figure 5). We also found that ratios comparing the inflorescence tissue within the different mutants (e.g., MSisAp1/AFSiAp1 and MSisAp3/AFSiAp3) or within leaf tissue from different treatments (e.g., MLesS04/AFLes04 and MLesS52/AFLes52) were not correlated as closely as those between mutants/tissues within the same technology platform (Figure 4). This indicates that comparisons of ratios are not always consistent enough across technologies to identify more subtle differences in a sample. However, there was sufficient consistency to distinguish ratios at a more coarse level (e.g, those involving completely different tissues). Also, the ratio data obtained from the Agilent microarrays which were non-self microarrays (each dye associated with a different sample) with two channels showed higher correlations than the ratio data from self-self microarrays (both dyes associated with the same sample) treated as single-channel microarrays. For example, the Pearson correlation between AFRoLe and AGRoLe was 0.91 while the correlation value between AFRoLe and AGRoLeZ was 0.83. Additional analyses using different Affymetrix and Agilent technical replicates also showed that each pair of ratio data from technical replicates were closely correlated as expected (see Figure 2 in additional file 1, e.g, AFLeS52A and AFLeS52B).

Probe-position Effects in Affymetrix and Tag-position Effects in MPSS

The position of the microarray probe or MPSS signature within a given transcript may impact the measurement of transcript abundance due to the use of different polyadenylation sites or technical biases [26]. Therefore, we compared the effects of the probe position on the Affymetrix microarray data and the effects of the tag position on the MPSS data. This type of analysis may be useful to distinguish differences in the measured transcript abundances that arise from biological causes in the mRNA sample from those that come from technological biases.

First, we investigated the relationship between the primary tag position and the measured RNA abundance in MPSS. Since MPSS tags are captured via the poly(A) tail and anchored by a four base restriction enzyme recognition site such as GATC, the primary tag was defined as the 3'-most tag in an annotated gene, which was determined based on an analysis of Dpn II sites. We analyzed the MPSS data considering the distance of the primary tag from the 5' and 3' ends. We found that primary tags closer to the 3' end tended to have lower measured transcript abundances (Figure 6C and 6D). With increasing distance of the primary tag from the 5' end of the gene, there was a decrease in the MPSS-measured transcript abundance (Figure 6A and 6B).

Next, we investigated the relationship between the probe position and perfect match (PM) probe values in the Affymetrix microarrays. The probe positions were determined relative to the annotated 5' and 3' ends of each gene. We found that the measured transcript abundances correlated with the position from the 3' end such that probes close to the 3' end tended to have lower PM values (Figure 7), but the variation was smaller than for MPSS tags. There was not as strong an association between the measured abundance and distance from the 5' end for the Affymetrix probes compared to MPSS.

Based on the analysis of probe position in the Affymetrix microarrays and the tag position in MPSS samples, we found that both probe position and tag position, as measured from the 3' end, demonstrated a consistent trend. With the increase in probe position or tag position from the 3' end, the PM value and tag-abundance tended to increase then plateau. One reason for this tendency could be that there are many alternative transcripts for some genes which are not correctly annotated in the current Arabidopsis genome which could be alternative splice products or transcripts utilizing alternative polyadenylation sites. In fact, some probes and tags that we used for this analysis were not located within the most recently annotated genome [27, 28]. To check whether this tag- or probe-position bias is not merely a consequence of alternative splicing, we were able to obtain similar results when only single exon-containing genes were included (results not shown).

Data Preparation and Pre-processing

The microarray data and MPSS data for the manuscript has been submitted to the GEO Omnibus [29]. The Series record number is GSE8994 and details relating to sample and data preparation are available there.

Eleven different distinct RNA samples were used in this study (Table 1). The RNA samples were obtained from different Arabidopsis thaliana tissues (samples 1 to 5 included callus, inflorescence, leaves, roots and silique), homeotic floral mutants (samples 6 to 9 included agamous inflorescence, apetala1-10 inflorescence, apetala3-6 inflorescence and a superman/apetala1-10 inflorescence double mutant) and two leaf treatments (samples 10 and 11 include leaves sampled 4 and 52 hours after salicylic acid treatment, respectively) [24]. Only one technical replication for each sample was available for the MPSS data. Two technical replicates for each sample were available for most of the Affymetrix and Agilent data (codes with last letter "A" or "B"). For the Agilent platform, some microarrays were self-self microarrays (e.g, AFCaA indicates both red and green channels are callus RNA samples) while others are were non-self microarrays (e.g, AGCaLeA indicates that red channel is callus and green channel is leaf RNA sample).

The Arabidopsis MPSS data has been published previously and is available on our website [30]. In the present study, we utilized the "Signature" MPSS samples, rather than "Classic" MPSS samples, because there exists a significant tag-position bias in Classic MPSS samples [26, 31]. We also focused on our analyses on the 17-base signatures, rather than the 20-base signatures, because the shorter tags represent a less-biased data set due to a lower frequency of palindromic sequences or other "bad words" [32] that are poorly sequenced by MPSS technology. The raw abundance of each signature in each sample is normalized to a metric of TPM (Transcripts Per Million), which facilitates comparisons across samples and technologies. To obtain the expression level for each gene, we summed the abundance of signatures that match the sense-strand of the gene and are found only once in the Arabidopsis genome and transcriptome.

To generate data that was comparable across technologies, the same Arabidopsis RNA samples used for the MPSS analyses were hybridized to two different high-density Arabidopsis microarray platforms: the Affymetrix ATH1 GeneChip [33] (22,800 features, ~11 probes per gene and 25-nucleotides per probe) and Agilent 22 K Arabidopsis microarray [34] (21,500 features, 60 nucleotide probes and one probe per gene). MAS5.0 software [33] was used to obtain the expression level for each gene on the Affymetrix microarray; Agilent data was preprocessed by Agilent G2566AA Feature Extraction Software [34]. Two technical replicates were obtained for almost all Affymetrix microarrays used to analyze each RNA sample from mutants/treatments/tissues (e.g, AFRoA and AFRoB) and for all Agilent microarrays (e.g, AGRoA and AGRoB which are self-self microarrays and AGRoLeA and AGRoLeB which are non-self microarrays). For the Affymetrix microarray data, for each gene the transcript abundance (expression) measurements from each pair of technical replicates was used to obtain the average (e.g., data was averaged for AFRoA and AFRoB to obtain AFRo). For self-self Agilent microarrays (both dyes associated with the same sample), because each microarray has two-dye channels, for each gene we first averaged the data from the red channel and green channels, then averaged these data from each pair of technical replicates to obtain the final average (e.g., data was averaged for AGRoA and AGRoB to obtain AGRo). Likewise, for Agilent microarrays which were non-self microarrays (each dye associated with a different sample), we averaged the pairs of technical replicates to get the average log ratio (e.g, data averaged for AGRoLeA and AGRoLeB to get AGRoLe).

We downloaded the most recent Arabidopsis genome release data (TIGR version 6, November 2005) from TAIR [35] and focused on the 26,751 protein-coding genes among the total of 31,407 annotated genes. Because each of these microarray platforms and the MPSS platform contained a subset of the total annotated Arabidopsis genes, we compared only the genes represented or detected by all three platforms. A total of 16,269 genes were potentially detectable by all three platforms (Figure 1). For the MPSS data, signatures which were duplicated (i.e., matching multiple genomic locations) were excluded from the analysis; genes matching these signatures were also excluded to avoid a bias that could result if only a subset of the signatures from each gene are used. For the comparison of MPSS and Affymetrix platforms, 6,501 genes were selected since these were detected as present in all eleven Affymetrix microarray data sets and had mean MPSS tag counts across all eleven samples greater than 4 TPM (i.e., the lower identification limit of MPSS signatures [32]). Moreover, for the comparison of the three platforms, 1,648 genes were selected with all MPSS tag counts greater than 4 TPM, detected as present in Affymetrix and detected as present in Agilent. The Agilent ratio data was produced in two different ways: the Agilent ratio data named without "Z" as the last letter was from a non-self microarray, with the log ratio derived from the red channel intensities divided by the green channel intensities; the Agilent ratio data with "Z" as the last letter was from self-self microarrays, with the log ratio obtained from two different microarrays. For example, AGRoLe data is from the root/leaf microarray, while AGRoLeZ data is from root/root and leaf/leaf.

Data Analysis

where x represented the normalized expression level of gene calculated by each platform. For MPSS platform, NEU has the same meaning as Transcript Per Million (TPM).

Contour-plots were used to represent the density of data. Expression data were log₂-transformed. When investigating the correlation relationship between different samples, hierarchical clustering method was used with both Pearson's correlation and Spearman's rank correlation as the distance metrics and using a complete linkage algorithm [36, 37]. The distance for clustering was calculated by 1 - r, where r represents correlation coefficient value. Correlation coefficient and total least-square values [38] were calculated for measuring the relationship among data from different platforms. A contour-plot matrix [25] was used to show correlation relationships. Because the relationship between two variables was often non-linear, we used a regression method called LOWESS (local weighted scatter-plot smoother) [39] to obtain smooth curves which showed the mean tendency of dependence pattern between two variables. When using LOWESS, we used a Gaussian kernel to calculate the weight contributing to each point from its neighboring points. We also used the bootstrap method [40] to evaluate the confidence intervals (5–95%) of the LOWESS method.