ReMOAT, re-annotation tables and Bioconductor

ReMOAT (Re-annotation and Mapping for Oligonucleotide Array Technologies) is a Web-based interface to the re-annotation data generated by the described pipeline. Data generated were processed using Perl and stored in a relational database. ReMOAT uses an Apache Web server, in a Linux environment, and a collection of PHP and Perl CGI scripts provide the user interface coupled to a MySQL RDBMS (http://www.mysql.com). The Web interface provides access to a tool capable of converting Illumina probe IDs, Entrez gene IDs or Ensembl gene IDs to several other formats such as Illumina probe IDs, Ensembl gene IDs, Entrez gene IDs, HGNC gene symbols, Unigene IDs or Lumi IDs. The resulting Web pages provide access to further information via links to Wikigenes (36), iHOP (37), HGNC nomenclature (38), Entrez Gene (39) and Ensembl (25) databases for each gene. A second search page enables extraction of the probe sequence, probe type, repeat masking results, GC contents, associated CpG islands, miRNAs and SNPs and quality assessment, generated by the pipeline for submitted Illumina IDs. Further search pages provide information on probe location, transcript and non-specific genomic hits (Supplementary Figure S1), and these are documented online.

Full re-annotation tables and their detailed description, as well as a link to ReMOAT, can be found online on http://www.compbio.group.cam.ac.uk/Resources/Annotation/.

A Bioconductor annotation package for each re-annotated platform has been built using the AnnotationDbi package. We chose reliable probes (i.e. with ‘Good’ or ‘Perfect’ annotation) and retrieved the RefSeq IDs as determined by our re-annotation. These IDs were then used by AnnotationDbi to create a database schema with predefined mappings. The packages can be installed by the same means as any other Bioconductor package and are named according to the platform annotated. For instance, the Human WG-6 V3 chip is available as the illuminaHumanv3.db (for data summarized using Illumina’s BeadStudio software) or illuminaHumanv3BeadID.db (for those starting from raw Illumina data) packages.

Annotation tables and ReMOAT are updated at least every 6 months, concomitantly with Bioconductor, and consistent version labelling is used for tracking old annotations.

Data sets

In order to evaluate the impact of probe annotation on the analysis and interpretation of Illumina expression data, we have looked at 27 Human WG-6 version 1 data sets (926 arrays) derived from a wide variety of conditions and tissue types (including breast, blood, artery, stem cell, sperm). These were obtained from GEO (40) on 18 November 2008, searching for platform GPL2507 and using the GEOquery Bioconductor package (41) to read the data into R (42).

The MAQC project (43,44) selected two human RNA samples for the comparison of differential expression on microarray platforms. These were the Universal Human Reference RNA (UHRR) from Stratagene (http://www.stratagene.com/manuals/740000.pdf) (a pool of 10 cancer cell lines from multiple tissues) and Human Brain Total RNA from Ambion (http://www.ambion.com/catalog/CatNum.php?AM6050). The original MAQC publication analyzed a dilution series of the two samples hybridized at three different locations. The Illumina portion of the data was generated using the Human WG-6 V1 platform and available in GEO (GSE5350) and ArrayExpress (E-TABM-132). For this study, we only consider the two extreme (pure) dilution levels. The same samples have since been run on Human WG-6 V2 BeadArrays by Asuragen Inc. and made publicly available on Illumina’s Website (http://www.switchtoi.com/datasets.ilmn). Finally, we have run those same samples in-house on Human WG-6 V3 BeadArrays.

Another illustrative data set consisted of Human WG-6 V2 arrays from the study in (45), which compared gene expression in hepatocytes between trisomic mice carrying human chromosome 21 and their wild-type litter-mates.

To illustrate possible effects of a SNP on the performance of a probe, we have used the Japanese HapMap (46) population. These 45 individuals have been examined using Illumina expression BeadArrays, and additionally genotyped, making them ideal for this purpose. Expression information comes from the Human WG-6 V1 platform and was obtained from the GENEVAR (1) project, and the corresponding genotypes were obtained from BioMart (24).

Two Illumina training samples were run on the DASL platform (47) using the standard cancer panel of 1506 probes. These consisted of RNA from the Illeum, and corresponding DNA (run in duplicate and triplicate, respectively).

Finally, we have looked at the preprocessed data from Miranda and colleagues, GEO series GSE13733. Samples consisted of MCF7 cells and were run on 10 Human WG-6 V3 arrays. There were two different drug treatments (DZNep and 5-Aza-CdR), with four technical replicates each, and two replicates of an untreated control sample.

Throughout the rest of the manuscript, these data sets are referred to as GEO, MAQC (V1, V2 and V3), Trisomy, HapMap, DASL and Miranda, respectively.

Filtering and differential expression analysis

Figure 2A shows how, for the GEO data set, the average expression ranks of genes vary between different annotation categories of probes. As could be expected, probes annotated as reliable are observed to have a high mean rank (~27 000 for ‘Perfect’ and ‘Good’), whereas other probes (e.g. matching to an intronic or intergenic region, or with no matches) have a much lower mean rank of below 20 000. The figure also shows a decreased mean rank of probes mapping closer to the 5′ end of transcripts, consistent with the instability of RNA molecules and their degradation typically starting from the 5′ end. This effect has been addressed in the probe-level models for the analysis of Affymetrix (20) and motivates most of the manufacturers, including Illumina, to design their gene expression platforms with probes targeting the 3′ end of transcripts.

Open in a separate window

Open in a separate window

Figure 2.

Impact of annotation on expression analysis. (A) Box plots of the average expression ranks of probes, calculated across the GEO arrays, for each annotation category (box widths proportional to the number of probes in the respective category); (B) proportions of probes retained by each filtering method applied to the MAQC V1 data set (y-axis) as a function of the chosen cut-off (x-axis); (C) proportions of probes of each category (y-axis) found in a gene list of certain length arising from a differential expression analysis of the MAQC V1 data (x-axis). For panels (B) and (C), the colour code is the same as used in (A); dashed lines are associated with quality grade categories and solid lines with coding zones.

Among the 100 bead types with the highest average rank (>46 891), 39 are found to target ribosomal proteins that function in protein biosynthesis and some of which have been previously found to be housekeeping or reference genes (55,56). Thus, these genes are often considered useful for normalization. Indeed, 4 genes from this list are also found in the table of top 15 candidate housekeeping genes presented in (55), which includes 13 ribosomal proteins. Other interesting features of these 100 probes are that 53 have multiple matches to the genome, 24 map to repeat regions and 41 map to SNPs (Supplementary Table S1).

For the MAQC V1 data set, Figure 2B shows that, for filtering based both on expression level (‘Detection’ and ‘Expression’) and variability (‘IQR’), probes annotated as ‘Perfect’ or ‘Good’ are more readily retained than those classified as ‘Bad’ or not matching.

Figure 2C shows, in terms of our annotation categories, the composition of the ranked gene lists from the differential expression analysis of the MAQC V1 data. As expected, bead types annotated as reliable occur towards the top of the list. On the other hand, genes mapping to intronic and intergenic regions, or having no match at all, are ranked lower in the list.

SNPs and mismatches

Having presented the general good performance of Illumina probe design, we now discuss some specific cases where the performance requires careful interpretation.

The retention, by the filtering methods, of probes with mismatches (named ‘Good’) and their generally high expression ranks suggest that hybridization is still possible when probe and target sequences do not match perfectly. Furthermore, for Human V2 over 100 of the detected mismatches can be explained by SNPs and 5791 probes map to annotated SNPs. Therefore, allelic ambiguity needs to be taken into account when interpreting gene expression data. Figure 3 and Supplementary Figure S3 illustrate an association between genotype and measured expression in the Japanese HapMap population that is consistent with a mismatch inducing a significant decrease in registered intensity. Similar effects have been previously reported for Affymetrix (57).

Open in a separate window

Figure 3.

Association between genotype and measured expression. (A) The sequence we examine is in the BCKDHB gene on chromosome 6, which contains a SNP (rs7740958—T/C) at position 7. The Illumina probe targeting this sequence (GI_34101271-I) contains a C at this location. There is also a probe (GI_34101266-A) targeting a constitutive splice junction of BCKDHB and matching no known SNPs. [Figure built based on UCSC Genome Browser graphics (26).] (B) Box plots of the log₂ expression ratios in the Japanese HapMap population according to the rs7740958 genotype.

As another example, human probes ILMN_1692545 and ILMN_1670800 (WG versions 2 and 3) differ in only one nucleotide (C/A at position 19) that corresponds to SNP rs13082444 (G/T for the genomic Watson strand). A striking difference in expression between the two is found for all the samples in the Miranda data set (median of 48 fold). According to the current human genome annotation, these probes primarily target an intergenic region and, therefore, no known transcript. However, these probes were designed to target transcripts XM_933970.1 and XM_944901.1 that perfectly map to the two genomic variants and have since been removed from the databases (Supplementary Figure S4).

We have listed, for the three Human WG-6 platforms, all the pairs of probes differing in 1, 2 and 3 nucleotides (Supplementary Table S2). There is a substantial increase in the number of probe pairs differing in one nucleotide from V1 (9) to V2 (177), suggesting that attention to SNPs and allele specific expression has been given in Illumina’s redesigning of probes.

Repeat sequences

It is not expected that intergenic and intronic regions of the genome will feature in gene expression studies. However, there are a few examples of probes that fall into these categories yet are highly ranked on all GEO arrays (outliers in categories ‘Bad’ and ‘Intergenic’ in Figure 2A). Most of such probes include repeat sequences and are partially or totally ‘masked’ by RepeatMasker (http://www.repeatmasker.org). They are likely to cross-hybridize and provide non-specific signal. We have, therefore, classified all ‘masked’ probes as ‘Bad’. A good example of how these probes can impact on the analysis of gene expression is provided by the study in (45), where our annotation has been used. As illustrated in Figure 4, all human transcripts but one overexpressed in Tc1 mouse livers (whose cells contain one human chromosome 21, apart from the normal mouse karyotype), when compared to their wild-type litter-mates, are either from chromosome 21 (as expected) or targeted by probes containing repeat sequences. Moreover, the number of overexpressed ‘masked’ transcripts is statistically significant. These results suggest that repetitive probes, although not originally annotated as mapping to chromosome 21, are non-specifically binding to repetitive transcripts from that chromosome. Assuming transcripts from human chromosome 21 are those truly differentially expressed between Tc1 and normal litter-mates, ignoring the information about repeat sequences would have led us to a false discovery rate of 37%, instead of the current 2%.

Open in a separate window

Figure 4.

Effect of repetitive sequences on gene expression measurements. Volcano plot (y-axis: empirical Bayes log-odds of differential expression; x-axis: log₂ fold change in expression) comparing the expression of human transcripts in livers between Tc1 mice carrying a human chromosome 21 and their wild-type litter-mates (Tc0). Blue dots depict probes targeting genes on human chromosome 21 and red dots probes comprising human RepeatMasker sequences. The dashed line is an arbitrary cut-off for differential expression based on the assumption that there are no human sequences transcribed in the wild-type mouse and, therefore, no human transcript should be overexpressed in Tc0. The P-values (Fisher’s exact test) show that the numbers of differentially expressed transcripts from human chromosome 21 (blue) and comprising human repeats (red) are both highly significant.

Another interesting example is provided by a set of six genomically clustered human probes, two of which were designed to target the PSMC1 gene while the other four match spurious transcripts that have since been removed from the databases. The six probes exhibit a good match not only with the PSMC1 locus but also with several other regions in different chromosomes (including an intronic sequence of the PLCZ1 locus) and the relative position of the matches is generally conserved (Supplementary Figure S5).

An even more striking example of potential non-specificity of hybridization is given by a set of 15 probes, each targeting transcripts from multiple genes from the GAGE cluster on chromosome X (Supplementary Figure S6).

Alternative splicing

Although for Illumina the majority of genes are targeted by one unique reliable probe type, there are still a few thousand genes covered by more than one probe type (Supplementary Figure S7). Almost all multi-exon genes undergo alternative splicing (58,59) and, therefore, comprise multiple isoforms. Probes targeting different transcripts for the same gene can generate ambiguity in the definition of gene expression but they also provide an opportunity for the detection of alternative splicing. This is well illustrated in Figure 5 and Supplementary Figure S8 for the human CAST (calpastatin) gene, covered by at least six different probe types for any of the platforms. The probes target the constitutive last exon, several alternative first exons and alternative splice junctions. The differences between probes, not only in expression levels but also in differential expression ratios, show the extra insight brought by a transcript/exon-centric analysis and how misinformative a gene-centric summary of data can be. This example is also an illustration of the importance of strand specificity of probes, as ILMN_1752145 targets the 3′-most exon of an antisense ERAP1 (endoplasmic reticulum aminopeptidase 1) transcript and genomically overlaps with a constitutive exon of CAST.

Open in a separate window

Figure 5.

Importance of alternative splicing on the interpretation of gene expression data. Schematics of the structure of the human CAST gene, with exons depicted by numbered boxes in light blue and all annotated alternative splicing events represented by the associated thin black dashed lines. The three upper tracks represent the logged expression levels for all the probes targeting the CAST locus for each of the three human Illumina WG platforms. Bars are positioned according to the relative position of the respective probes in the locus. Bar height is proportional to the average log gene expression in the MAQC data: red for brain samples and blue for reference samples. The coloured dashed lines indicate the respective average expression levels of probes targeting CAST for each sample type. The short coloured full lines in the middle of the V2 track indicate the expression level of CAST estimated by the original Illumina gene summarization procedure for each sample type; black stars identify V2 probes targeting CAST according to the original Illumina annotation. Black arrows identify probes mapping to the reverse strand. The three lower tracks represent the corresponding log-ratios. This figure shows that a gene-centric analysis indicating underexpression of the CAST gene in brain (dashed blue line) might be biased by a brain-specific skipping of exon 30 and/or an alternative first exon.

Probes matching splice junctions can be particularly effective in the detection of alternative splicing events. We clearly show, for the DASL platform, the increase in relative signal intensity for such probes from hybridization against genomic DNA to hybridization against transcriptomic RNA (Supplementary Figure S9).

The interpretation of expression data can also be affected by ambiguity in the definition of a gene or, in other words, by distinct genes sharing the same locus. For example, two different human genes, CPNE1 (copine I) and RBM12 (RNA-binding motif protein 12), share their most 5′ exons (as well as the promoter region) (60). Three Illumina probes (ILMN_1701229, ILMN_2276000, ILMN_2276002) target both an alternative non-coding first exon of CPNE1 and the last and only coding exon of RBM12 (Supplementary Figure S10). An even more striking example is the PCDHG locus on human chromosome 5, which seems to cluster several distinct genes and different isoforms of the same gene (generally characterized by alternative first exons) (Supplementary Figure S11). Illumina designed probes to target virtually all the different transcripts and annotates each of them as a distinct gene. In contrast, UniGene (34) considers them to be all part of the PCDHGC3 gene cluster. This illustrates well the difficulty in choosing a suitable universal gene identifier and that even a system like UniGene, specifically designed for organizing transcripts into non-redundant gene-oriented clusters and, particularly, popular in meta-analyses (61), cannot deal with all the idiosyncrasies associated with defining what a gene is.

Another example of the particular importance of analyzing expression at the transcript level for genes containing alternative promoters and first exons is given by the CDKN2A (cyclin-dependent kinase inhibitor 2A) locus for both human and mouse. The gene has two alternative first exons, separated by over 10 Kb, that give rise to two structurally and functionally distinct isoforms: p16INK4a (inhibitor of CDK4 kinase) and p14ARF (alternate open reading frame, involved in the stabilization of p53 by sequestering of MDM2). All the Human and Mouse WG-6 platforms have three probes targeting CDKN2A: one for each of the alternative first exons and one for the common last exon (Supplementary Figure S12). Annotating the trio of probes as having the same gene target would, for most experiments, lead to the observation of inconsistency between the associated three gene expression measurements.

Any signal given by probes matching an intron or exon junction and, therefore, partially intronic, like human ILMN_1690644 (Supplementary Figure S13), is very difficult to interpret. In particular, this holds when there is no evidence for splice variants involving an alternative splice site or a retained intron.

Generally, probes targeting constitutively transcribed sequences are expected to exhibit stronger signal than probes matching only alternative transcripts. The former are, therefore, more suitable for a gene-centric analysis of data. We support this assumption by showing that there is a positive correlation between the proportion of isoforms in a gene targeted by a probe and its intensity (Supplementary Figure S14).

We acknowledge the support of the University of Cambridge, Cancer Research UK, and Hutchison Whampoa Limited. We thank Illumina, Bin Liu, Carlos Caldas, Natalie Thorne, Leonard Goldstein, Inma Spiteri, David Allard, Michelle Osborne, James Hadfield, Mahesh Iddawela, Ana Teresa Maia, Ashraf Ibrahim, Roslin Russell, Christina Curtis, Mike Wilson, Mike Smith, Dan Tomlinson, Gos Micklem, Gordon Smyth and Lynn Amon for their support and feedback.