Exon Array Analysis using re-defined probe sets results in reliable identification of alternatively spliced genes in non-small cell lung cancer

Data set and specimens

All investigations were performed in accordance with the Declaration of Helsinki. After informed consent had been obtained, biopsies of tumour tissue and of healthy mucosa from patients with lung cancer undergoing elective curative surgery from 2001 to 2003 at the Department of Surgery, Charité, Campus Benjamin Franklin, Berlin, Germany, were collected by the clinical partner with approval by the local ethics committee (reference number EA4/006/05). Cancer and corresponding epithelial normal tissue specimens were snap frozen in liquid nitrogen within 20 min following excision and stored at -80°C. All tissue samples were evaluated by a pathologist before and during macrodissection to ensure an enrichment of vital tumour cells and normal epithelial of 80 - 90%. Total RNA preparation was performed with the RNeasy^® kit (Qiagen, Hilden, Germany). In addition, total RNA isolated from clinical samples was obtained from Cambridge Biosciences, Ltd. (Cambridge, UK). In total, 11 paired samples of adenocarcinoma of NSCLC and NAT (n = 22) and seven paired samples of squamous cell lung carcinoma and NAT (n = 14) were analysed (for details see additional file 1). All NSCLC samples have at least 60% tumour content or were enriched by macrodissection. Concentration of total RNA was quantified using a NanoDrop spectrometer and quality of RNA was measured on an Agilent Bioanalyzer (Agilent RNA 6000 Pico and Nano Kit; Agilent Technologies, Santa Clara, CA, USA). All samples had an RIN (RNA integrity number) ≥ 6.9.

Synthesis of fragmented cDNA and microarray hybridisation

Sample amplification and preparation for microarray hybridisation was performed according to the GeneChip^® whole transcript sense target labelling assay (manual rev. 4, Affymetrix). In brief, in this procedure 2 μg total RNA of each sample were depleted from rRNA (RiboMinus™ Transcriptome Isolation Kit (Human/Mouse), Invitrogen, Carlsbad, CA, USA), reverse transcribed to cDNA, amplified by in vitro transcription, and reverse transcribed to cDNA. Fragments between 40 and 70 bp were generated enzymatically from 5.5 μg cDNA, labelled terminally, hybridised onto the Affymetrix Human Exon 1.0 ST array, stained and washed, and finally scanned using an Affymetrix GeneChip^® Scanner 3000 7G. Quality of scanned array images was assessed visually using GeneChip^® operating software (GCOS). The raw data (CEL-files) have been deposited in the ArrayExpress Archive under accession number E-MEXP-2644.

Pre-processing and data analysis

In the analysis of the exon array data, we included only probe sets derived from the core set. These probe sets correspond to exons that are part of well-annotated transcripts, e.g. RefSeq transcripts [27]. However, in addition, the exon array contains probe sets that are backed by less evidence such as EST data and Genscan predictions, and thus we omitted these probe sets. All pre-processing operations were conducted in Partek^® Genomics Suite version 6.5 (Partek Inc., St. Louis, MO, USA). Probe intensities were quantile normalised, a background adjustment was performed as implemented in Partek^® Genomics Suite using a model based on the probe sequence, and finally probe intensities were summarised using the RMA algorithm. All probe set intensities were normalised using quantile normalisation. For quality assurance and in order to detect outlier samples, the overall gene expression was analysed by means of a principal component analysis (PCA) and with hierarchical clustering based on Pearson's correlation (average linkage). Quality assurance with PCA and hierarchical clustering was also performed using gene level expression values. To this end, probe set expression values were summarised to gene level expression values using Tukey's biweight algorithm as implemented in Partek^® Genomics Suite.

where cg is the contrast of the pathology group (G) effect between tumour and NAT derived from the respective MLM.

To discover tumour sub-type dependent differences in the alternative splicing patterns, the levels of the factors pathology group G were changed from tumour and NAT to AdCa, SCC, and NAT. Expression ratios were obtained as the contrast of the pathology group (G) effect between tumour and NAT from a MLM ANOVA based on samples of each subtype, respectively. Likewise, SI values were calculated for each subtype. Information about the location of alternative splicing events identified from well-known transcript variants was obtained from UCSC genome browser (database hg19, table knownAlt, table schema described under http://genome.ucsc.edu/cgi-bin/hgTables?db=hg19&hgta_group=genes&hgta_track=knownAlt&hgta_table=knownAlt&hgta_doSchema=describe+table+schema, data available for download under http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/knownAlt.txt.gz).

Background filter

Probe sets with a signal near the background noise are termed absent. We estimated the expression above background noise using the MAT algorithm [39] and utilised the implementation of MAT in ProbeEffects version 2.0.1. The software was modified in such a way that detection above background p values are generated per probe set instead of per gene as in the original implementation. The MAT background model is fitted per sample to genomic and anti-genomic background probes (GC-bins). Probe sets are considered to measure a signal above the background if p ≤ 0.01. A probe set is treated as present if p ≤ 0.01 in at least 75% of the samples of the respective pathology group (tumour or NAT). Otherwise, a probe set is treated as absent in the respective sample group. Only probe sets that are absent in both sample groups, i.e. absent both in tumour as well as in NAT, are filtered out. Genes containing fewer than five present probe sets are also removed.

Re-definition of probe sets

New core level probe sets were defined using most up-to-date transcript data. From publicly available sources (Ensembl and Vega version 56.37a, RefSeq release 38), in total 269611 transcript variants related to 42429 gene loci were obtained. Non-redundant exons of all transcript variants were collected.

Overlapping exons cause splitting of such exons into distinct probe selection regions (PSRs). In contrast to the approach followed by Affymetrix, we do not include less reliable EST data and Genscan predictions for the generation of PSRs. Thus, we try to avoid extensive fragmentation of exons into small PSRs. It can be expected that PSRs will be larger and that each will be covered by more probes than in the Affymetrix probe set definition. Mapping information of exon array probes to the genome was obtained from X:Map version 56.37a [46]. Probe sets were created for PSRs that are covered by at least one probe. Probes mapping to multiple locations on the genome as well as PSRs from repetitive regions were filtered out (RepBase version 9.11, RM database version 20050112 [50]). The new probe set definition together with annotation are available as additional file 2.

Laboratory validation of exon array results

Of each RNA sample, 50 ng total RNA were reverse transcribed in cDNA and amplified linearly according to manufacturer's instructions (WT-Ovation™ RNA Amplification System, NuGen Technologies, Inc., San Carlos, CA, USA). Endpoint reverse transcription PCR (RT-PCR) was carried out for AdCa/NAT paired samples in 20 μl reactions starting from 100 ng NuGen amplified template cDNA. Amplification using Immolase™ DNA polymerase (ImmoMix and ImmoMix Red, Bioline GmbH, Luckenwalde, Germany) was performed for 30 cycles at 60°C annealing temperature. Assays were designed to yield two distinct amplicons, one of smaller size skipping the cassette exon, the larger one including the enframed cassette exon (see additional file 3). HPLC purified primers were obtained from TIB MOLBIOL GmbH (Berlin, Germany). PCR products were separated by electrophoresis using MetaPhor^® Agarose (Lonza, Basel, Switzerland), cut out of the gel, and extracted (QIAEX^® II Gel Extraction Kit, Qiagen, Hilden, Germany). The extracted products were ligated into a TOPO^® pCR4^® vector (TOPO TA Cloning^® Kit For Sequencing, Invitrogen, Carlsbad, CA, USA) and transferred into chemically competent TOP10 cells. Plasmid DNA was isolated and sequenced with M13-reverse and M13(-20)-forward primers at AGOWA GmbH (Berlin, Germany). Subsequently, the PCR amplicon sequences were assembled and checked against the transcript sequences annotated in Ensembl 56.37a.

Transcript variants were quantified using quantitative reverse transcription PCR (qRT-PCR) in paired samples of AdCa or SCC and NAT. qRT-PCR assays were designed using the Universal Probe Library Assay Design Center (Roche Diagnostics GmbH, Mannheim, Germany; https://qpcr.probefinder.com/) or the software Primer3 [51]. HPLC purified primers were obtained from TIB MolBiol GmbH. Assays were designed either to yield transcript variant-specific amplicons (transcript variant expression) or to yield amplicons in exons common to all known transcripts (gene level expression). Assay details are given in additional file 4. Where possible, assays were designed to span an exon-exon junction such that contaminations of genomic DNA have no influence. qRT-PCR was carried out in 10 μl reactions starting from 10 ng NuGen amplified template cDNA. Amplification was performed with MESA BLUE qPCR MasterMix Plus for SYBR^® Assays (Eurogentec, Seraing, Belgium). All qRT-PCR reactions were carried out in an Applied Biosystems ABI Prism^® 7900 HT cycler with 40 cycles at 60°C annealing temperature. A dissociation curve was recorded after the run. All assays were checked with regard to efficiency and PCR product size using cDNA from cell lines. After each run, obvious outliers (no amplification, multiple products) were manually excluded from the subsequent analysis. Ct values were set in the logarithmic phase of the PCR-reaction. Averaging of technical replicates and normalisation to the geometric mean [52] of the endogenous control gene Ct values (ESD, POLR2A [53]) results in ΔCt values. Both fold-changes (FC) as well as splicing index (SI) values were calculated from ΔΔCt values (see equations 5 and 6). Significance was determined by means of a paired t-test on the SI values.

We consider exon array results to be validated successfully if qRT-PCR results show a significant (paired t-test p ≤ 0.05) or high (SI ≥ 2.0 or SI ≤ 2.0^-1) difference of transcript variant expression that is in agreement with our hypothesis of AS and the exon array results. Validation is successful if qRT-PCR results fulfil these criteria for at least one subtype of NSCLC (AdCa or SCC). In addition, if an RT-PCR result is available, it must exhibit the AS pattern in the majority of the samples analysed.

Expression of FOX-genes

Gene expression levels of FOX1 (A2BP1) and FOX2 (RBM9) in different types of cancer and healthy body tissues were obtained from a collection of a total of 1015 samples hybridised on the Affymetrix GeneChip^® array HG-U133_Plus_2.0. Signal intensities were summarised and normalised with the MAS5.0 algorithm [54]. Mean values and standard deviations per sample group were calculated with log-transformed probe set intensities.

Expression of FOX1 and FOX2 transcripts was detected with RT-PCR as described above. Primer sequences were: FOX1, forward 5'-CCAGTTGTGGGTGCAGTCTA-3' and reverse 5'-AGCTTCCTTTCTCCCCACAT-3'; FOX2, forward 5'-GCGGACAGTATATGGTGCAGT-3' and reverse 5'-GGTTGCAGTAGCAGGCTGTG-3'. Alternatively spliced exons in RT-PCR products of FOX2 were identified after cleavage with exon specific restriction enzymes (FastDigest^® Hpy8I, AfeI, and HinfI; Fermentas, Glen Burnie, MD, USA). DNA fragments obtained from RT-PCR and optionally after digestion with restriction enzymes were separated by agarose gel electrophoresis as described above or using the Agilent Bioanalyzer (Agilent DNA 1000 Kit; Agilent Technologies, Santa Clara, CA, USA).

FOX1 gene expression, FOX2 gene expression, and FOX2 transcript variant expression was quantified using qRT-PCR as described above. Primers were: FOX1 gene level, forward 5'-CCTTACCTTCCTGGACTGATTC-3' and reverse 5'-GTAAGGCTGAGCCATTGTGTC-3'; FOX2 gene level, forward 5'-TGGAAATTAAGCCCAGTAGTTG-3' and reverse 5'-TGATACCCCCTCTTCCTGA-3'; FOX2 transcript variants, common forward primer 5'-GCGGACAGTATATGGTGCAGT-3'; FOX2 cassette exon, reverse primer 5'-TAGAGGTCAGCACCGTAAAATCC-3'; FOX2 exon skipping, reverse junction primer 5'-CATATCCACCCCTGGATAGG-3'.

Background filtering reduces the number of false positive results

We utilised the generally accepted analysis of variance (ANOVA) method in order to identify gene loci affected by differential splicing. A false discovery rate (FDR) of 0.05 corresponds to an ANOVA p value of 0.018 in the NSCLC data set. According to this analysis, 5340 candidate genes are affected by alternative splicing. Of the genes showing a p value close to zero (p < 1.4 · 10^-45), we manually inspected the top 100 list with the most extreme SI, and assigned them to one of six classes according to their expression profile (Figure 2). Although this classification has not been verified and may contain some errors, it will help us to detect key features of an analysis based on ANOVA and SI alone. Representative gene profiles are shown in Figure 2. It became evident that only 30% of all gene loci in the top 100 list are true positives (Figure 2c and 2d). All of the other candidates appear to be false positives (Figure 2b, e, and 2f).

In particular, more than half of all gene loci in the top 100 list exhibited probe sets with a low expression value in both pathology groups (Figure 2e). We assume that these probe sets are absent in both pathology groups, i.e. the corresponding exon is expressed neither in tumour nor in NAT. These probe sets will only measure the background signal in the respective sample group and thus are non-informative. Still, their expression value affects the statistical analysis: the FC of absent probe sets does not follow the gene level FC. The statistical ANOVA method scores genes containing such background level probe sets with a low p value which leads to the high rate of 57% false positives. Therefore, we introduced a data set-specific background filter that identifies and removes probe sets that are absent in both sample groups before starting any statistical analysis (see Material and Methods).

After applying our background filter, we repeated the ANOVA analysis for the identification of candidates differentially spliced between NSCLC and NAT (FDR := 0.05, p = 0.02). This time, we identified 3414 candidate genes. We compared this result set with the result set obtained from the analysis without background filter (5340 candidate genes). There is an overlap of 2505 genes between both sets, i.e. most genes were found with both analysis methods. First, we inspected the set of genes that were removed in the analysis with background filter. Then, we inspected candidate genes that were found only in the analysis with background filter.

In total 2835 genes were removed by the background filter. Of these, 1965 genes were completely removed since less than five of their probe sets were found to be present. Of the remaining 870 genes, in total 4388 absent probe sets were removed by the background filter (31% of all core probe sets of these genes). As a consequence, these genes are no longer identified by the ANOVA method as being differentially spliced. In summary, the background filter removes 53% of all candidate genes from the result list.

However, also 689 additional genes were found to be differentially spliced only after applying the ANOVA analysis with background filter. We hypothesise that the majority are statistically borderline such that their inclusion is due to the effect of setting a cut-off value (here: FDR := 0.05). To test this hypothesis, we analysed the AS ANOVA p values of genes removed by the background filter and genes that were only found with the background filter (in total 1559 genes). p values were taken from the ANOVA analysis both without and with the background filter. For each gene individually, the difference of both log-transformed p values was calculated (Δ log p). Under the assumption that our background filter has no effect and all differences are merely an effect of setting a cut-off value, all Δlog p values should be small and evenly distributed. We observe a significant skewness in the Δlog p value distribution (skewness of -1.77, D'Agostino skewness test p < 2.2 · 10^-16). This means that only genes removed by the background filter show highly different p values. Apparently, the background filter has the expected effect and other changes in the result list can be attributed to setting a cut-off value.

Inappropriate use of annotation file leads to a high false negative rate

Affymetrix provides meta probe set (MPS) files that contain relationships between probe sets and genes. MPS files are commonly used (e.g. several implementations in commercial analysis software) to define the set of probe sets constituting a gene, but this is not intended by Affymetrix and leads to a loss of many candidates for differential splicing. We compared this file with supposedly redundant information from more comprehensive chip design annotation files. It became evident that a gene does not contain all available core level probe sets in the MPS file. Affymetrix provides these annotation files for the analysis of differential gene expression and reports that probe sets less suitable for this purpose were left out. We found that by inappropriate use of the annotation file in total 220 differentially spliced genes were lost as false negatives. Our results suggest that by using the original MPS files for the analysis of differential splicing, the false negative rate is unnecessarily increased.

Further analysis improvement by chip re-annotation with most up-to-date transcript data

As outlined above, it is important to base any analysis of differential splicing upon a proper definition of probe sets and gene loci. In contrast to the identification of novel splice variants, our focus is on finding AS events that are well-annotated by transcript variants in public databases and that are differentially expressed in NSCLC versus NAT. It is intuitive to make use of the best-available transcript annotations. Thus, we created a new chip definition based on recent transcript data (Figure 3).

In order to assess the improvements by the new chip definition, we investigated whether events of differential splicing can be reliably identified. We chose a straightforward approach of manually inspecting gene expression profiles that were generated on the basis of the original or the modified chip definition. The impact on the detection of differential splicing can be most easily assessed with genes where an AS pattern is not overlaid by a differential gene expression pattern. Therefore, we focused on genes that are only affected by differential splicing but that do not exhibit differential expression on the gene level (e.g. Figure 2c). A list of genes was created based on the difference in the absolute SI value: the top end shows results that were gained and the bottom end shows results that were lost by the new chip definition. Expression profiles of the top 10 and the bottom 10 candidates were inspected manually.

It became apparent that almost all of the candidates only identified with the new chip definitions show differential splicing for exons not covered by probe sets in the Affymetrix chip definition. These exons were not annotated as well-known exons before, but are meanwhile supported by new transcript annotations. In most cases, the new probe set contains at least two probes, in one case even nine probes. Thus, additional information on transcript annotations also leads to additional identification of differential splicing events. The probe set of another gene was found to be present in both chip definitions. Only with the new chip definition does this probe set indicate differential splicing. Closer investigation reveals that this probe set contains four probes in the Affymetrix chip annotation. Three of these probes map to multiple locations on the genome and thus were left out in the new chip definition. The new probe set contains only a single, but more reliable probe that now indicates differential splicing. On the top 10 list, we found one case where five genes were combined into a single transcript cluster according to the new chip definition. A single RefSeq transcript annotation spans over five gene loci and thus causes their combination. Differences in the expression levels of the five genes leads to extreme SI values. This example demonstrates that an exon array analysis is sensitive to the definition of transcript clusters. Errors in transcript annotations can lead to artefacts and the example here could represent a false positive result.

For candidate genes that are no longer identified with the new chip definition, we found that in most cases adjacent probe sets were combined to form a single new probe set. The new probe set measures an averaged signal which leads to a reduced SI or no evidence for differential splicing. In two cases, the evidence for differential splicing based on the Affymetrix chip annotation was provided by probe sets containing a single probe. In three cases, probe sets of the Affymetrix chip definition were completely removed. These probe sets contain either a single probe or overlapping probes and furthermore, almost all probes map to repeat masked regions and thus were left out in the new chip definition. Together, combination of probe sets and removal of less reliable probes leads to a new chip definition that is best suited to detect differential splicing of known exons at the expense of the ability to detect novel exons. Our improved workflow now consists of a background filter and a new chip definition which is based on up-to-date transcript annotations. At the same time, the number of artefacts due to cross-hybridisation of probes is reduced.

Application of the new workflow to reliably identify alternative splicing events

From a previous laboratory validation approach using candidates from the original workflow (Affymetrix chip definition, MLM ANOVA p value) we achieved only a low validation rate. In order to reliably identify genes that are differentially spliced in NSCLC versus NAT, we utilised our enhanced workflow (background filter, summarisation using the new chip definition, MLM ANOVA). Based on all available data, differential splicing candidates were selected by an informed decision-making approach in order to generate a result list for laboratory validation. In particular, we inspected known transcript variants.

Starting with our new chip definition, we omitted all genes containing fewer than five probe sets (leaving 268132 probe sets that correspond to 16171 genes). After applying the background filter, we generated a primary result list based on the AS ANOVA p value (FDR := 0.05, p = 0.02). This generated a primary result list of 3096 genes (additional file 6).

Because this list is too large for manual inspection, sub-lists were generated based on additional criteria. Sub-list A focuses on genes that are differentially spliced but that are not affected by differential gene expression (see category in Figure 2c). Here, we filtered by pathology group p > 0.05 and took the top 100 genes with most extreme SI (additional file 7). As we know that most of the alternatively spliced genes are also affected by differential gene expression (see category in Figure 2d), we generated sub-list B by filtering pathology group p ≤ 0.05. Only genes showing high homogeneity in all tumour and all NAT samples (i.e. no dependency on the patient, patient p > 0.2) were included in this sub-list. Furthermore, we focused on genes where the gene expression FC was not higher than the SI. We took the top 100 genes with most extreme SI (additional file 8). Thirdly, we generated sub-list C by focusing on probe sets that are known from transcript annotations (UCSC genome browser, database hg19, table knownAlt) to be affected by AS. Genes were included in this sub-list if these specific probe sets had a SI ≤ 1.4^-1 or SI ≥ 1.4 (in total 194 genes; additional file 9). All expression profiles of genes in any of these sub-lists (Venn diagram Figure 4; additional file 10) were inspected manually in order to formulate a hypothesis on the mode of AS (e.g. cassette exon). Based on all available information, 14 promising candidates were selected for validation by independent methods in the laboratory.

Laboratory validation of exon array results

Firstly, we tried to validate all of the cassette exon events, the intron retention event, and the alternative 5'-splice site event using RT-PCR in AdCa samples. In most cases, we identified two distinct PCR products, one originating from the short and the other originating from the long transcript variant of the respective gene (see additional file 11). Both the product size as well as the estimated abundance in tumour versus NAT were in concordance with the exon array results and with the AS event hypothesised by us. However, in two events (CTNND1-CE-3', intron retention event of KIAA1217) only two of six patients exhibited the expected result. In two cassette exon events (CTNND1-CE-5', MYH14) merely a single PCR product was found. Together, this leads to eight positive and four negative results. In all cases, sequencing of representative PCR products confirmed the expected exon-exon junction sequence according to our AS hypothesis.

Secondly, we tried to quantify the expression of different transcript variants with qRT-PCR. For all of the 16 differential splicing events, transcript variant specific primers were designed. The SI was calculated separately for AdCa and for SCC (based on the qRT-PCR results; see additional file 11). Four events of differential splicing could not be confirmed in AdCa (ANTXR1, CTNND1-CE-3', intron retention event in KIAA1217, MYH14). The remaining events showed a difference of transcript variant expression that is in agreement with the exon array results. High difference (SI ≥ 2.0 or SI ≤ 2.0^-1) was found for ten events (ADD3, CLSTN1, CTNND1-CE-5', FN1, alternative transcription start site in KIAA1217, MYO18A, NCOR2, NUMB, SYNE2, TPM1-MX-5') and of these, the results of six events achieved significance (ADD3, FN1, alternative transcription start site in KIAA1217, NCOR2, NUMB, TPM1-MX-5'). A smaller yet significant difference was found in two events (SLK, TPM1-MX-3'). In total, 12 events of differential splicing were confirmed in AdCa. In SCC, eight events exhibited differential splicing that is in agreement with the exon array results (ADD3, CTNND1-CE-5', FN1, alternative transcription start site in KIAA1217, NUMB, SLK, TPM1-MX-5', TPM1-MX-3'). Although the difference in transcript variant expression was high in all cases (SI ≥ 2.0 or SI ≤ 2.0^-1) the significance level was not reached. Three results are borderline (ANTXR1, CLSTN1, and MYO18A) as there was evidence for differential splicing, but neither with significance nor with high magnitude. qRT-PCR results of NCOR2 were inverse compared to exon array results (details see additional file 12). Four events of differential splicing were not confirmed with qRT-PCR in SCC (CTNND1-CE-3', intron retention event in KIAA1217, MYH14, SYNE2).

We consider exon array results to be validated if both RT-PCR and qRT-PCR results are in agreement with our hypothesis of AS and with the exon array results (details see Materials and Methods). Based on our validation criteria, 11 events were validated successfully whereas five events were devalidated (see Table 1). Of the latter, two events (CTNND1-CE-3', intron retention event of KIAA1217) showed evidence for differential splicing in some patients but failed to meet our relatively strict criteria. Likewise, differential splicing of one cassette exon event of CTNND1 (CTNND1-CE-5') was confirmed with qRT-PCR but not with RT-PCR and thus did not meet the validation criteria.

Expression analysis of FOX splicing factors in NSCLC

Interestingly, many of our validated CE events were previously described as cancer-specific splice variants in breast and ovarian cancer [55]. It was shown that these and AS events of other genes depend on the splicing factor FOX2, which in cancer is either down-regulated or expressed as an inactive splice variant. We tried to elucidate whether FOX2 regulates cancer-specific AS also in other types of cancer such as NSCLC.

First, we analysed expression of the splicing factors FOX1 and FOX2 using a collection of in total 1015 gene expression microarrays. This database consists of twelve types of cancer and corresponding normal tissue as well as 50 other healthy body tissues. FOX1 was only expressed in brain, heart muscle, and skeletal muscle, but neither in cancer nor in other healthy tissues (additional file 13). RT-PCR and qRT-PCR results confirmed no FOX1 expression in NSCLC and NAT (data not shown). FOX2 was down-regulated in tumours of endocrine origin (thyroid, breast, fallopian tube, ovary, cervix, and prostate) but up-regulated in kidney, oesophagus, stomach, and colon cancer. In brain and lung cancer, the FOX2 expression did not change compared to normal tissue (Figure 5, additional file 14).

To confirm the expression levels and to investigate AS of FOX genes in lung cancer, we examined the respective exon expression profiles in our NSCLC exon array data set. All probe sets corresponding to exons of FOX1 were found to be absent (data not shown). None of the FOX2 probe sets showed any indication of differential expression or differential splicing (Figure 6). The exon array results confirmed the data from our gene expression microarray database.

Although no evidence for a change between NSCLC and NAT was obtained, we wanted to prevent a false negative result and analysed FOX2 expression using PCR based methods. For the analysis of AS of FOX2, we focused on exons in the C-terminal region whose skipping was shown to generate an inactive form of the splicing factor [56]. We noted that two cassette exons are annotated in this region. In a previous study of the mouse orthologue of FOX2, yet another cassette exon designated as M43 was described [57]. An orthologous open reading frame in the intron sequence of human FOX2 is also present and we included it as a putative exon in our analysis. Primers for RT-PCR assays were placed in the flanking constitutive exons. The experimental results showed that two transcript variants were present both in NSCLC as well as in normal lung tissue (Figure 7). The smaller RT-PCR product resembled an amplicon lacking all cassette exons whereas the longer product allowed the inclusion of only one of the three cassette exons. All cassette exons were of similar size, hence exon-specific restriction enzymes were used in order to identify the exons involved by cleavage. DNA fragments after restriction digestion indicated that only the first of the three cassette exons was expressed. This was the same exon that had previously been reported as being skipped in breast cancer [55].

The overall FOX2 expression levels and expression levels of each of the two transcript variants were quantified using specific qRT-PCR assays. It became evident that the longer FOX2 transcript variant, which includes the first cassette exon, was the dominating form and was expressed up to six-fold higher than the smaller transcript variant. But neither the FOX2 gene expression nor expression of any transcript variant differed significantly between NSCLC and NAT (Figure 8).

Quantification of ADD3 transcript variants in NSCLC

One of the successfully validated genes that showed a significant validation result is now described in more detail: ADD3. ADD3 contains at least 16 exons of which exon number 15 (ENSE00000986819) is a known cassette exon of 96 bp. In Ensembl, 22 transcript variants of ADD3 are annotated (Figure 9a). The exon array probe set that maps to the cassette exon exhibited a high SI in the NSCLC exon array data set (Figure 9b-d). This indicates that the cassette exon is preferably included in cancer, but skipped in normal tissue. This differential splicing pattern was detected in 13 of 18 patients (72%). Four patients showed an inverse pattern and one patient did not exhibit any indication of differential splicing. There was no correlation between these differences and factors such as subtype or staging. RT-PCR of AdCa samples confirmed the differential splicing pattern measured using the exon array in four cases (Figure 9e samples 1 - 4). Furthermore, in two patients (samples 5 and 6) AS clearly took place, but no difference was observed between tumour and NAT. These were the same patients that showed an inverse result or no differential splicing on the exon array. Each transcript variant was quantified with qRT-PCR using specific primer pairs. Moreover, the overall gene expression of ADD3 was measured using a primer pair targeting a constitutive gene region. qRT-PCR results (Figure 9f-g) confirmed that the overall amount of ADD3 did not change between tumour and NAT; thus, FC and SI values were similar. Both in AdCa as well as in SCC a switch from exon skipping to exon inclusion was observed such that the cassette exon was expressed more than 2-fold higher in tumour compared to NAT. Likewise, the exon skipping transcript variant was expressed about 2-fold lower in tumour versus NAT; however, the latter difference was only observed in AdCa, not in SCC. Based on these results, we postulate that ADD3 lacking the cassette exon is expressed in normal lung tissue. In NSCLC, a switch occurs to exon inclusion while the overall gene expression remains unchanged.

Alternative splicing patterns that correlate with the NSCLC sub-type

The laboratory validation results demonstrated that we could reliably identify genes that are differentially spliced from an exon array data set using the improved workflow. Finally, we investigated whether there are genes whose splicing pattern shows a clear difference that correlates with the NSCLC subtype, i.e. AdCa and SCC. Our data set contains a balanced set of AdCa and SCC samples. In the ANOVA model, the factor subtype (levels: AdCa, SCC, and NAT) was used instead of factor pathology (levels: tumour and NAT). In total, 4285 genes showed evidence of differential splicing in at least one subtype (FDR := 0.05, p = 0.026; additional file 15). For all samples and subtypes, we calculated separate SI values. A high difference of SI values between AdCa and SCC is an indication of a subtype-specific pattern of AS. There were 16 genes that exhibited a subtype-specific AS pattern (Δ log₂SI ≥ 0.5 or Δ log₂SI ≤ -0.5). In each gene, we analysed the probe set with the highest difference in SI values between AdCa versus NAT compared to SCC versus NAT. With regard to these probe sets, we found that three genes exhibited a more pronounced differential splicing pattern in AdCa (GTPBP10, DUOX1, and KIAA1217), whereas 13 genes had a more pronounced pattern in SCC (e.g. CD44, FN1). In the case of KIAA1217, the respective probe set covered the exon that had been quantified using qRT-PCR. In agreement with the exon array expression values, also qRT-PCR results showed an increased expression of the short transcript variant that was more pronounced in AdCa versus NAT compared to SCC versus NAT.

Data set

Clinical samples show high variability. This can be attributed to patient-to-patient variability and in the case of cancer, to differences in tumour content and especially to considerable differences between individual tumours. We only used matched pairs of NSCLC and NAT in order to level out patient-to-patient variations. Tumour content was high in all NSCLC samples and a sub-set of the samples was enriched by macrodissection. The data set is balanced such that both major sub-types of NSCLC (AdCa and SCC) and staging are present in similar proportions. Although detailed sample annotations are available, finding AS patterns specific to each factor remains a question of power. With 18 sample pairs our NSCLC data set has the power to discriminate between tumour and NAT on a level similar to other cancer data sets with matched samples. While Xi et al. focused on the AdCa sub-type of NSCLC [34], our data set also contains the SCC sub-type. This allows identification of AS patterns common to both AdCa and SCC. In addition, one can find sub-type specific AS patterns.

Workflow improvements

It is equally important to use an analysis workflow with a high detection rate and a justifiable false positive rate. We found that the standard workflow leads to a low true positive rate of 30% caused primarily by artefacts. We introduced an enhanced analysis workflow consisting of (1) filtering of probe sets that are absent in both of the sample groups, (2) an updated chip annotation, and (3) a statistical analysis based on a mixed linear model (MLM) ANOVA and SI.

A modified version of the model-based algorithm MAT was used in order to identify absent probe sets. Rather stringent cut-off values were set with regard to the MAT algorithm p value as well as homogeneity in one sample group in order to flag a probe set as absent. Thus, genes with low expression and signal intensities close to the background noise might have got lost as false negatives. Nevertheless, these background filter settings resulted in a pronounced decrease of artefacts which we consider to be essential for a reliable result list.

We showed that the analysis can be further improved by an updated chip annotation based on recent transcript data. Thus, differential splicing of new exons that were not annotated as well-known at chip design time can now be detected. In contrast to the chip definition created by Affymetrix, we did not include EST data or predicted exons in order to avoid excessive clustering of exons. Probes that map to multiple locations to the genome as well as probes mapping to repeat-masked regions were left out during the reannotation process in order to avoid artefacts caused by cross-hybridisation. We found evidence that this approach leads to probe sets with a higher number of probes and thus a more reliable measurement of exon expression. Different parts of an exon that has a predicted internal splice site can obviously not be discriminated. Our chip definition is best suited to analyse differential splicing of known transcript variants in contrast to the detection of novel transcript variants. It should be noted that the new chip definition can be used in combination with any algorithm of choice and within other analysis workflows than the one described in this study. All files are compatible with the respective Affymetrix file format and no changes in software implementations are necessary.

Thirdly, we introduced an informed decision-making approach for the identification of genes that are differentially spliced. A MLM ANOVA yielded the significance of AS, differential gene expression, and information about heterogeneity between the patients. In addition, we estimated the gene level expression value directly from the MLM by calculating a contrast. We assume that this is a straight yet robust approach for reliably calculating the SI. Potential candidate genes that are affected by differential splicing were identified using all information available including SI, significance, known events of alternative splicing, and transcript annotations. In effect, we used advanced statistical tools to generate a primary result list. Based on all information available, we selected genes from this list and inspected their exon expression profile manually.

The enhanced workflow was applied to our NSCLC exon array data set. A validation of the representative candidates using independent laboratory methods led to a 69% confirmation rate of exon array data. In addition, we found genes that are fairly known to be affected by AS in cancer, e.g. FN1. We did not adjust our workflow to give better results for well-known cancer genes. That these kind of genes can be found on the result list confirms that enhancements led to a more reliable identification of differential splicing in general.

Differential splicing in non-small cell lung cancer

In this study, we found that 330 genes are affected by differential splicing in NSCLC compared to NAT. We investigated their expression profiles and moved on to the next step, the formulation of a hypothesis of AS, i.e. definition of the exons involved and the mode of AS (e.g. cassette exon). Of the primary result list, 14 candidate genes were selected and validated by the independent laboratory methods RT-PCR and sequencing. In total, validation results confirmed 11 events of differential splicing in NSCLC. One of these events represents usage of an alternative start of transcription which is not strictly a mode of AS and might be regulated at a different level than splicing. The other events, however, are genuine events of AS. Several of the successfully validated genes have already been shown to be differentially spliced in cancer. Our data confirmed a switch to the onco-foetal EDB antigen in FN1. CLSTN1, which is a transmembrane and cell adhesion protein, was found to be affected by AS in colon, breast, bladder, and prostate cancer before [32, 58]. Dutertre et al. showed that expression of CLSTN1 transcript variants can be utilised in order to give a prognosis of metastasis-free survival in breast cancer [58]. We found that exon skipping in CLSTN1, which Dutertre et al. demonstrated as being associated with metastasis in breast cancer, occurs in NSCLC. FN1 and CLSTN1 were also reported in another study as being differentially spliced in breast and ovarian cancer versus normal tissue [55]. Closer investigation of the primers used in the RT-PCR revealed, however, that exons other than the cassette exons that we found being differentially spliced in NSCLC were affected. CEACAM1 is a transmembrane protein that is involved in signal transduction. One alternatively spliced cassette exon leads to a short isoform. Xi et al. showed qRT-PCR results indicating that the short isoform of CEACAM1 is over-expressed in both AdCa and SCC of NSCLC compared to NAT [34]. In our exon array data set, the signal intensity of the respective probe set was at the background level and hence no conclusion could be made. For another cassette exon in CEACAM1, however, our exon array data set indicated a higher skipping rate in NSCLC compared to NAT. Tropomyosin 1 (TPM1) contains two pairs of mutually exclusive exons. Recent transcript annotations in Ensembl 56.37a provide evidence that the first pair is not strictly mutually exclusive since both exons were found together in two transcript variants. TPM1 is a tumour suppressor that is involved in cytoskeleton remodelling and that has pro-apoptotic properties [59]. With regard to the first exon pair (exon 2 and exon 3), our results provided evidence for a switch from exon 2 to exon 3 in NSCLC. The exon pair 7A/7B was found to be differentially spliced in different kinds of cancer (colon, bladder, kidney, lung, prostate, and uterus [60, 27, 32]). Our results confirmed that the mutually exclusive exon pair 7A/7B of TPM1 is differentially spliced in NSCLC such that there is a switch from exon 7B to exon 7A in tumour. The alternative 5'-splice site in one exon of NCOR2 leads to different gene products that are known as SMRT-α (long exon) and SMRT-τ (short exon) [61]. In this study we demonstrated that in normal lung tissue both transcript variants of NCOR2 are equally expressed; however, at least in NSCLC AdCa a strong shift towards the longer transcript variant SMRT-α was observed. In KIAA1217, usage of an alternative start of transcription leads to a shorter transcript that is in-frame. To our knowledge, cancer-specific usage of an alternative start of transcription in KIAA1217 has not been reported before. In the genes ADD3, MYO18A, NUMB, SLK, and SYNE2, AS events lead to in-frame variations in the protein sequence affecting from 15 to 48 amino acids.

FOX2 activity is not regulated at the transcript level in NSCLC

Venables et al. found that many cancer-specific AS events in breast and ovarian tumours are caused only by alterations of the splicing factor FOX2 (RMB9) [55]. Among other genes, this was shown for ADD3, MYO18A, NUMB, and SYNE2. FOX2-dependent AS was demonstrated for SLK in a study of human embryonic stem cells [62]. Thus, many of the cassette exon events of which we confirmed cancer-specific AS in NSCLC are known to be regulated by FOX2. This led us to analyse whether FOX splicing factors (FOX1 or FOX2) could possibly explain the changes in splicing patterns between NSCLC and normal lung tissue. FOX1 shows only limited expression in healthy tissue [63, 64] and is expressed neither in breast nor in ovarian tissue [55]. Our results demonstrated that FOX1 is not expressed in normal lung tissue or in NSCLC, either. Two distinct mechanisms were identified by which FOX2 modifies the splicing pattern of target genes in cancer [55]. In ovarian cancer, FOX2 expression is significantly down-regulated compared to normal ovarian tissue. FOX2 itself is also affected by AS in breast cancer and a shift towards a functionally inactive splice variant was observed before. Our gene expression microarray data confirmed that the overall FOX2 expression level is down-regulated in ovarian cancer compared to normal ovary. Also other tumours of endocrine origin in our microarray database showed this kind of down-regulation (thyroid, breast, fallopian tube, cervix, and prostate). Expression of FOX2 was, however, not changed between lung cancer and normal lung tissue. Furthermore, other types of cancer (kidney, oesophagus, stomach, and colon cancer) showed an even higher expression of FOX2. Our analysis of splice variants revealed that the functional variant as well as the inactive form of FOX2 is expressed in lung tissue. A change of the splicing pattern between NSCLC and normal lung tissue was not observed. This means that there is no shift to the inactive FOX2 transcript variant in NSCLC in contrast to breast cancer. These findings suggest that FOX2 activity is not regulated at the transcript level in NSCLC. With the present data, regulation of FOX2 activity at the protein level cannot be excluded. Another explanation is that FOX splicing factors, unlike in other types of cancer, are not responsible for the cancer-specific splicing patterns in NSCLC. FOX2-dependent regulation of AS might be linked to steroid hormone systems that have a significant influence on reproductive system tissues like ovary, cervix, breast, and prostate. This can be explained for instance by the known interaction of FOX2 protein with the oestrogen receptor and transcription factor ER (ESR1) [65].

The long isoform of the cytoskeleton protein ADD3 is cancer-specific

It was also noted by Venables et al. that many of the cancer-associated genes found in their study are functionally related to remodelling of the cytoskeleton and cell movement [55]. Most of the candidate genes of which we confirmed cancer-specific AS in NSCLC are also involved in these processes (ADD3, CLSTN1, FN1, MYO18A, NUMB, SYNE2, and TPM1). This is another indication that AS potentially plays a major role in tumour progression, especially in tumour invasion and metastasis. We analysed AS of one gene that is functionally involved in the cytoskeleton in detail.

Gamma-adducin (ADD3) is a structural constituent of the spectrin-actin cytoskeleton that contains one known cassette exon. We found that this cassette exon has higher inclusion in NSCLC compared to NAT for the majority of patients. Dutertre et al. found that AS of the same cassette exon is associated with metastasis: they analysed AS in primary tumours generated from four different mouse mammary tumour cell lines with different potential for metastasis. Cassette exon inclusion in the mouse ortholog of ADD3 is high only in tumours originating from cell lines with a high potential for metastasis [58]. Interestingly, ADD3 was also detected in acute lymphoblastic leukemia as a fusion protein with NUP98[66]. We investigated the gene structure of ADD3 in more detail in order to examine potential functional consequences of AS. ADD3 encodes the γ polypeptide chain of adducin. Adducin occurs as a heterodimer or heterotetramer of α adducin together with either β adducin or γ adducin. The function of adducin oligomers in cells is capping (+)-ends of actin filaments (f-actin) and recruitment of spectrin [67]. All adducin monomers have a similar structure: a globular head domain that participates in subunit interaction, a neck domain, and an unstructured tail domain. At the C-terminus, a sequence with homology to myristoylated alanine-rich C kinase substrate (MARCKS) protein can be found that contains the spectrin/actin interaction site as well as calmodulin and protein kinase C (PKC) targets [68]. In ADD3, the cassette exon can generate an insert of 32 amino acids in the tail domain that is upstream of the MARCKS-related domain. Protein structure predictions from Ensembl 56.37a provided information that this insert leads to a small coiled coil formation followed by a low complexity region. With the present data, we cannot estimate what effect the amino acid sequence insert generated by AS has on either the oligomerisation of adducin monomers or the binding of other proteins like f-actin, spectrin, calmodulin, and PKC.

While our exon array analysis and the validation process focused on finding events of AS that were present in the majority of patients in our study, we also found evidence for considerable heterogeneity, e.g. four of 18 patients showed an inverse AS pattern for ADD3 and an inconclusive result for NCOR2 in SCC. We gave an outlook as to how additional factors like tumour subtype can be considered in our workflow. In total, 16 genes showed an AS pattern that is differentially pronounced in AdCa versus SCC. Yet, inclusion of additional factors is a matter of power and hence study design.