A Critical Evaluation of Analytic Aspects of Gene Expression Profiling in Lymphoid Leukemias with Broad Applications to Cancer Genomics

In cancer research, transcriptional aberrations are often deduced from mRNA-based gene expression profiling (GEP). Although transcriptome sequencing (RNA-seq) has gained ground in the recent past, mRNA-based microarrays remain a useful asset for high-throughput experiments in many laboratories. Possible reasons are the lower per-sample costs and the opportunity to analyze obtained GEP data in association with published data sets. There are established and widely used methods for the analysis of microarray data, which increase the comparability of different GEP data sets and facilitate data-mining approaches. However, analytic pitfalls, such as batch effects and issues of sample purity, e.g. by complex tissue composition, are often not properly addressed by these standard approaches. Moreover, most of these tools do not capitalize on the full range of public data sources or do not take advantage of the analytic possibilities for functional interpretation or of comprehensive meta-analyses. We present an overview of the most critical steps in the analysis of microarray-based GEP data. We discuss software and database query solutions that may be useful for


Introduction
Traditionally, gene expression analysis includes reverse transcription of mRNA into cDNA and probing of gene transcripts of interest by specific primers designed for target PCR amplification (gold standard), followed by quantitative, semi-quantitative (e.g.qRT-PCR), or electrophoresis (e.g.Southern blotting) detection methods.Based on efforts provided by the Human Genome Project [1,2] and studies on expressed sequence tags (ESTs) in mammalian genomes, cDNA hybridization array chips have originally been designed to investigate deregulated mRNA expression of distinct and well-characterized gene transcripts in various diseases.Modern mRNA-microarray platforms apply one or two-color fluorescence labeling (i.e.Cyanine3 / Cy3 for green and Cy5 for red dye fluorescence) for one or two samples to be loaded on the chip, respectively, and allow the detection of more than 47 000 transcripts.In contrast to two-color arrays (e.g.HuA1 by Agilent Technologies, Santa Clara, CA, USA), one-color arrays, are most commonly used today (e.g.HG-U133 Plus 2.0 by Affymetrix, Inc., Santa Clara, CA, USA, or BeadArray HT-12v4, Illumina, Inc., San Diego, CA, USA) and represent the focus of this review.
The past few years have seen the advent of transcriptome sequencing (RNA-seq) based on the next-generation sequencing (NGS) technology using high-throughput platforms, such as the GA IIx or HiSeq2000 sequencer from Illumina.RNA-seq does not require the prior design of specific probes, rendering it a highly versatile approach for gene expression profiling (GEP).Accordingly, a number of publications on the genomic landscape of various neoplasms have applied RNA-seq to investigate gene-specific aspects such as differential splicing and exon usage [3], hidden viral transcripts [4], and cancer-specific fusion transcripts [5].However, published reports using RNA-seq in cancer often lack statistical power for comprehensive gene expression analyses due to a limited sample size.In contrast, mRNA-based microarrays have remained the initial method of choice for high-throughput analyses of gene expression in many laboratories.Reasons for this include the associated lower per-sample costs as well as the availability of already published microarray-derived GEP data in
public databases.Many of these data sets were processed by established and widely used methods, thereby improving their comparability and the suitability for data-mining approaches.
Within this review, we present an overview of critical steps in the analysis of microarray-based GEP data (see overview in Figure 1) and the corresponding library and code information (summarized in Table 1 and 2).We will discuss step-by-step software and database query solutions that may be useful for data analysis, to avoid analytic pitfalls, and to provide an increased capability for clinical and biological interpretation of data.To illustrate the proposed analytic steps, we present analyses on exemplary data of previously published and own GEP data, all obtained in patients with B-and T-cell leukemias or lymphomas.

Quality Control can Greatly Differ by Platform
There are various possibilities to apply basic steps of quality control (QC) prior to or during preprocessing of GEP raw data.In order to avoid false estimates of background intensities and false inputs for normalization, removal of potential problematic samples and probes before data preprocessing is essential towards a correct interpretation of data.Problematic samples often present as outliers in density distributions or in an unsupervised cluster analysis on global gene expression values (after data preprocessing).The latter, e.g. in form of dendrograms (Code 1) or principal component analyses (PCA; Code 2), is created by using the R [7] library arrayQualityMetrics from Bioconductor [8] with its informative HTML report per array.
Numerous methods and libraries for R are available for more specific quality assessments for each of the three major microarray platforms.Affymetrix arrays can be analyzed using the affyQCReport and simpleaffy libraries (see Table 1 for all library references), which normalize expression values using housekeeping genes (e.g.calculating the actin3/actin5 ratio), while the affyPLM library allows calculation of important quality measures such as the normalized unscaled standard error (NUSE) and relative log expression (RLE) as well as their plotting across samples (Code 3).The quality of data obtained with Illumina chips can be assessed by statistical standard measurements (mean and standard deviation) or outlier detection using the lumiQ function within the lumi library (Code 4).Possible slide inhomogeneities (i.e.scratches) or contamination on two-color arrays may be detected with the imageplot function of the limma library.This package also allows the calculation of the RNA Integrity Number (RIN) as a measure of mRNA degradation with a subsequent option to remove samples below a given threshold.

Proper Preprocessing of Raw Data
A first step in the standard analysis protocol of cDNA microarrays usually is the conversion of hybridization image spots obtained by array scanners into raw gene expression values.For Affymetrix chips this is normally done either by using the freeware Affymetrix Power Tools or the R library affy.For Illumina's BeadChips the proprietary GenomeStudio software or manual decryption via the R library beadArray may be used.For two-color arrays, scanner output files, e.g. in TIF format, can easily be read with the read.maimagesfunction from the limma R library.
In a second step, background correction is conducted by subtracting technical noise from biological variation.This is accomplished by using e.g.RMA [9] for Affymetrix arrays or the bgAdjust function from the lumi R library for Illumina arrays, which employs a similar algorithm as GenomeStudio (Code 5).In order to account for outliers and to remove systematic variation, normalization of expression values is required.The most common procedures include quantile-normalization, which preserves the rank, but may eliminate small differences in expression values, and LOESS (locally weighted scatterplot smoothing)-normalization, which does the opposite.
Robust splice normalization (RSN) aims to combine the advantages of both methods through a AIMS Medical Science Volume 3, Issue 3, 248-271.
monotonic splice fit to one reference sample, while simple scaling normalization (SSN) forces samples to have the same scale and background.Both approaches are included in the lumi R library for Illumina arrays.For two-color arrays it may be essential to further account for dye biases in the normalization [10] and to normalize within the array itself (between both color-labeled samples) and between all two-color arrays of the cohort, e.g. by use of the limma R library.Variance-stabilizing normalization (VSN) constitutes another method for combining background correction and normalization [11], while preserving biological variation.It is implemented in the vsn (Code 6) library, applicable to arrays of all major platforms.Within the normalization process raw intensities are usually transformed, either into a log2 scale or glog in case of VSN, in order to smoothen extreme values.

Probe Annotation and Deconvolution
Frequent impediments for GEP data analysis are missing array annotations or outdated annotation files provided by the manufacturers (e.g.frequently old GenBank predictions are included).Data-mining tools such as biomaRt [12] can be used to acquire up-to-date probe information (Code 7).They may also be helpful in assigning probes to transcripts, thereby enabling filtering for redundancies of probes, which map primarily to transcripts that are prone to nonsense-mediated mRNA decay (NMD) or to unprocessed pseudogenes.Deconvolution of genes with known transcript variants of differential function into probed isoforms may also be important for extrapolations on biological relevance.An example is the apoptosis regulator myeloid cell leukemia sequence 1 (MCL1), of which the longer isoform (MCL1-001) has been reported to enhance survival by inhibiting apoptosis, while its shorter isoform (MCL1-002) acts as a pro-apoptotic molecule [13].

Exploring Differentially Expressed Genes Considering the "Multiple Comparisons
Problem" Raw data preprocessing and QC is followed by the actual statistical analysis, usually in the form of probe-by-probe hypothesis tests for differential expression including: ( non-parametric).However, statistical testing of all genes / transcripts detected by an array requires correction for multiple testing, in order to avoid a substantial number of false-positive findings [14,15].
For example, using a significance level of 0.05 for each of 10,000 tests would result in approximately 0.05 * 10,000 = 500 significant rejections by chance, even if all null hypotheses of no differential expression were true.To this end, we can either control the family-wise error rate (FWER) to curtail the number of statistically significant results, e.g. by use of the (conservative) Bonferroni correction, in which the significance level for each probe-specific test equals the FWER (e.g.0.05) divided by the total number of tested probes, or by some permutation / resampling approach.Furthermore, we can aim for controlling the false-discovery rate (FDR), i.e. the proportion of falsely rejected null hypotheses, e.g. using the Benjamini-Hochberg's procedure, q-values, or other approaches.It should be noted, however, that control of the FDR, while very helpful in limiting the number of erroneously followed-up probes, does not imply a notion of statistical significance.The procedures by Bonferroni and by Benjamini-Hochberg are implemented in the multtest library [16], while the qvalue library provides an implementation for the rank-preserving q-value calculation (Code 10).
Nominally differentially expressed probes (e.g. with a single-test level of p < 0.05) can also be filtered by multiple-testing correction, for example by applying a q-value / FDR cutoff (common cut-off, e.g.0.1) to ensure a low proportion of false-positives in the set of probes to be subsequently followed up.To reduce time in the analysis, it may also be useful to exclude genes / probes that are not expected to be differentially expressed either due to biologically low variability in the investigated samples, or due to technically low detectability on the array.This can be achieved either by non-specific filtering of expression values restricted to a given range (e.g. the shortest interval containing half of the data by standard deviation (sd) or interquartile range) or by setting an empirical cut-off to the coefficient of variation (sd/mean), e.g. the top 10 percent or a fixed value of 0.6.Note, however, that this may increase the rate of false-negative findings (Code 11).

Pitfalls: Batch-correction and Contamination Estimation
When comparing GEP data obtained in the same laboratory, but with two or more different batches of arrays, the results will deviate from one another beyond the expected biological and arrayspecific technical variation.Batch correction addresses this issue.Two approaches commonly considered to be performing best [17] are mean-centering and a Bayesian framework named ComBat [18] (Figure 2a-c; Code 12).
A particular problem for cancer transcriptomics / genomics is the contamination of cancer tissues by normal cells (irrespective whether to consider them as actual milieu components) and vice versa.Even in lymphomas and lymphoid leukemias, such problems are encountered in lymph-node samples or in the seemingly 'pure' blood samples, as these are also of mostly multicellular composition.Tools like ESTIMATE [19] can weigh specific markers (e.g.indicating an immune or stromal cell origin) within gene expression profiles in the form of gene set enrichment analyses and thus evaluate the degree of purity.Unfortunately, due to intrinsic aberrations of 'immune cell' genes within tumor cells of leukemias / lymphomas, the immune gene set used within ESTIMATE is not reliable for the enrichment analysis within these malignancies (Figure 2d; Code 13).An alternative approach especially for leukemias / lymphomas might be CellMix [20] which uses gene sets from  When comparing the black dot to gray dots (all other samples), one can observe that the sample is among those with highest purity.Lower panel: sample among those with lowest purity.

Making Use of Public Databases
Two public databases are commonly used for the comparison of own microarray data with independent data sets, for example in a meta-analysis, namely the GEO (gene expression omnibus) database [21] (http://www.ncbi.nlm.nih.gov/geo) and the ArrayExpress database [22] (https://www.ebi.ac.uk/arrayexpress), with GEO featuring a larger number of integrated samples.

Implemented queries include:
• "Query 1: Get experiments where the sample description contains diabetes" • "Query 2: Get differentially expressed genes where factor is asthma" • "Query 3: Show expression for ENSG00000129991 (TNNI3)" • "Query 4: Show expression for ENSG00000129991 (TNNI3) with its GO annotations from Uniprot (Federated query to http://sparql.uniprot.org/sparql)" • "Query 5: For the genes differentially expressed in asthma, get the gene products associated to a Reactome pathway" • "Query 6: Get all mappings for a given probe e.g.A-AFFY-1/661_at" Query 2 and 5 can be further modified in order to compare gene dysregulation in other types of diseases, e.g. in lymphoid leukemias, such as chronic lymphocytic leukemia (CLL; Table 3).User's familiarity with the underlying ontologies (controlled vocabulary; [24]) is, however, necessary to construct queries.

Meta Analyses: Exploring Possible Phenotypic Markers across Different Conditions
For conceptualizing a pharmacologic compound (e.g.inhibitor) acting against a specific gene product or for designing specific gene-knockouts within a model organism, it may be particularly important to know in what conditions and disease subtypes expression of a distinct gene is up-or down-regulated and to which degree (basal or extreme).Integrative analyses of expression changes within a multitude of samples of the same entity, or model organism, or any other comparable biological system as well as across initially separately analyzed (and published) series (cohorts) are often called gene expression meta-analyses.In the following we describe multiple ways to conduct a meta-analysis of GEP data with their limitations and advantages.
The first approach includes construction and sending of specific queries to the EMBL / EBI RDF platform.Querying can further be semi-automated using the SPARQL R library, which allows the investigation of different data sets in a specific condition, e.g.comparisons of CLL vs. normal B-cells, or between distinct groups of tumor samples stratified by a characteristic of interest, e.g.
Since not all 'ArrayExpress' data sets are yet integrated into the EMBL / EBI RDF platform and the GEO database contains additional data sets, the manual download, processing, and integration of such additional data is often necessary.
Therefore, a second, more hands-on approach to meta-analyses is a search by keyword, e.g.
'chronic lymphoid', within GEO and / or ArrayExpress (or any other public database).Once the data set has been picked, it is background-corrected and the annotated replicates can be combined with their original samples by calculating their mean.Afterwards all samples within the data set are normalized (e.g.quantile-normalized).
Probe sets of a gene which map to retained / dysfunctional transcripts (or which map to more retained / dysfunctional transcripts than other probe sets of the same gene) should be removed to obtain meaningful expression values (Suppl.For further evaluation of the GEP meta-analysis, three different techniques for integration can be used to observe gene expression patterns and entity clustering: 1) The first method quantile-normalizes a matrix of average gene expressions across entities from different experiments and finally gives a visual approximation.If there is also a tumor suppressor gene (very low expression) and an oncogene (very high expression) in the gene set to be evaluated, one can expect an expression range similar to the whole transcriptome.It should be noted that in previous Affymetrix sets, such as HG-U133A, some genes (e.g.BMF and BOK) are not covered by specific probes on the array and, therefore, need to be imputed by the median of the respective data set.This guarantees that in the heatmap (or PCA) these genes are not visualized as up-or down-regulated; they in fact can be manually labeled (blackened).Expression values from all data sets are merged into one matrix and again quantile-normalized to account for variability in platform specifications and noise.A more suitable approach than normalizing on each gene set separately might be to normalize on the whole combined transcriptome (intersection of all probed genes).However, this would disregard genes not covered by all platforms used.The resulting heatmap (generated by function heatmap.2,library gplots; Figure 3b) shows the expression of selected genes and transcripts in their respective data set and can be additionally subdivided by the different entities (median across samples of an entity).
2) Batch effects cannot be entirely excluded by using method 1) as may be observed by a bias in clustering of samples from the same experiment.Therefore, we recommend a novel method called inSilicoMerge [25], which combines data sets and removes their batch effect with a choice of various methods, such as the empirical Bayes method ComBat (Figure 3c).
Unfortunately, data sets from different platforms can only be combined gene-wise, meaning that e.g.MCL1 would not be deconvolutable into its isoforms MCL-001 / MCL1-long and MCL-002 / MCL1-short.
3) For an advanced evaluation, one can further perform differential expression analysis for data sets with different control samples (of varying quality, number, and specificity) available for comparison, such as 'normal' non-malignant cells or bulk tissue specimens.Fold-changes with a p-value < 0.1 (trend) or < 0.05 (significant) are extracted to compare normal-matched gene expression between different experiments and probe targets representing different gene transcripts or protein isoforms.The results are again visualized by a heatmap, either in the order obtained by hierarchical clustering (using Euclidean distance) or in order of rows sorted by gene name.
As exemplified by illustration of expression levels of Death-Associated Protein Kinase (DAPK) gene family members in subsets of CLL and normal B-cells (Figure 3d), this method allows different disease vs. 'normal' comparisons and facilitates the evaluation of which genes are exclusively down-or up-regulated and which show no clear pattern or which are specific to small subgroups.In the meta-analysis itself every differential expression analysis is further evaluated by statistical testing.
Default setting is the Student's t-test, except for low variation or non-normal distributions, for which the non-parametric Wilcoxon rank sum test is recommended.

9.
Functional Analyses: the More the Merrier In the abundance of genes obtained as significantly dysregulated, the role or function of a specific gene is often unknown and it is therefore encouraged to group them functionally by software tools often coined as 'pathway analysis' or 'enrichment' tools.One of the most user-friendly, however, costly tools is QIAGEN's Ingenuity® Pathway Analysis (IPA®, QIAGEN Redwood City, www.qiagen.com/ingenuity/).Users can upload their differential expression results in the format of Excel tables into the Java GUI (graphical user interface).Annotation in the form of chip design or symbol identifiers (such as Gene Symbol, Ensembl ID or GenBank ID) can be selected for a given column as well as statistical parameters in separate columns, such as p-values, fold-changes, q-values / FDRs or simply expression values (fluorescence in microarrays or FPKM (fragments per kilobase of exon per million reads mapped) for RNA-seq).The list can be further restricted to a given range (e.g.p-value < 0.05).The selected genes are subsequently assembled into manually curated biological or toxicological / pharmacological pathways provided with an E-value (chance of a random hit).One advantage of IPA compared to other tools is the easy visualization of results by intuitive geometric forms, i.e. nodes / genes are drawn as distinct geometric symbols and edges / protein modifications in distinct line types.Similar graphs can be drawn with igraph in R, but are restricted to users that are more experienced in bioinformatics.
Other user-friendly and open-source alternatives include DAVID [26], gene set over-representation analysis (GSOA) by ConsensusPathDB [27] (Suppl.Figure 2), and gene set enrichment analysis (GSEA; Figure 4a) by the Broad Institute [28].All three tools can be operated from web GUIs, while the first two options also offer an R implementation or in the case of GSEA, also a JAVA desktop application.
For more advanced users and those seeking to work with protein identifiers (complementary to above mentioned tools) STRINGdb10 [29] is a potential alternative.Within the R library PPI (protein-protein interaction) graphs (nodes colored according to fold-change and also reachable via  or of other molecular features (i.e.mutational or cytogenetic strata in CLL) (Figure 7c-f; Suppl.Figure 5); measured by ANOVA for numerical or by entropy for categorical values.When looking for a cut-off for adverse prognosis, they can be further used in the form of regression trees [41].Different parameters can be controlled in this approach, such as the maximum size of a tree or the number of portions / bins.It is recommended to keep these relatively low in the training set to avoid "overfitting" and thus enable re-evaluation in the test set.Random forests [42] (as an assembly of permutated decision trees) can be used to determine the chance of observing random tree branching (library randomForest) (Code 21).Both algorithms are also included in the rattle library, which offers a user-friendly GUI with interactive plots and a selection menu for class variable and co-variates as well as algorithm and parameter choices.For a more detailed review on current machine learning algorithms in GEP, we refer to [43].
AIMS Medical Science Volume 3, Issue 3, 248-271.amplification, MYC mRNA upregulation, ATM mRNA downregulation, and TCL1A mRNA upregulation.ATM deletion status seems to be the most informative co-variate, however due to the excessive size of the tree (controlled by pruning and number of bins) there is a risk of "overfitting".f) Shown is a more feasible and smaller decision tree.Again, the most informative co-variate seems to be the status of ATM gene deletion.Followed by AGO2 amplification status.This is further confirmed in random forests (permutated decision trees) in order to circumvent 'overfitting' (not shown).

Discussion
In this review we discuss procedures to optimize GEP analyses.We highlight the importance of advanced preprocessing, such as batch correction and admixture modeling, but also appraise the versatility and sophistication of analysis and classification algorithms.Many of the presented

Figure 1 .
Figure 1.Flow chart describing a suggested GEP protocol.Steps in yellow boxes are modular and may function somewhere independently downstream of the steps in grey boxes.The red text refers to those figures of this review that illustrate the respective step.

Figure 2 .
Figure 2. a) PCA (principal component analysis) of the 1000 most variable genes (by variation coefficient) within 12 distinct batches of our T-PLL (T-cell prolymphocytic leukemia) data set reveals batch-specific clustering.b) After batch correction samples do not cluster anymore due to technical bias, but rather due to biological information when annotated as in c).c) Entity information can be included in ComBat (besides batch information) to fit batches.T-PLL samples (further divided by different oncogene protein status) and normal T-cells form a cloud, natural killer T-cell lymphoma seems to be most similar to T-cell prolymphocytic leukemia (T-PLL) and hepatosplenic T-cell lymphoma (HSTL).ii) Variables factor maps (produced by libraries like FactoMineR) show what marker contributes (or correlates) the most to each principal component and thus carries the highest specificity.Platform overlap was reduced to gene level, then batch-corrected using ComBat and quantile-normalized. d) Example taken from [51].Fold-changes were calculated according to labeled comparisons for each Death-Associated Protein Kinase (DAPK) gene family member, then the range was cut off and results were visualized.Color bars used for 3 distinct comparisons: (1) CLL vs. normal B cells (various subtypes); (2) CLL with IGHV unmutated vs. mutated gene status; (3) CLL with post-to-pretreatment and other clinical comparisons.
web link) and enrichments (including p-values and number of observed and expected interactions)AIMS Medical ScienceVolume 3, Issue 3, 248-271.A popular supervised learning approach are support vector machines[40] (SVM; R libraries gmum.r or e1071).They try to separate classes by projecting features and their interactions into high-dimensional space and subsequently by searching for either linear (Figure6a-b) or non-linear (Figure6c; Suppl.Figure4) separating hyperplanes in the original feature space (Code 19).

Figure 6 .
Figure 6.a) Support vector machine (SVM) classifies samples of T-cell prolymphocytic leukemia (T-PLL) based on TCL1A protein status (positive, intermediate, negative; by flowcytometry) predicted by TCL1A and TCL1B mRNA expression.As one can see in the top left two samples are misclassified by SVM as TCL1A-negative (red, but squared symbols).b) SVM of T-PLL samples of different TCL1A protein status ("dim" being intermediate) by numerous mRNA markers performs more robust classification.c) Example of a linear (upper panel) and a non-linear, radial / polynomial fit (lower panel) of a SVM.T-PLL samples which carry the ATM gene in mutated vs. unmutated constitution are classified by their status of ATM deletion and AGO2 amplification.Results, as seen by approximate pattern in linear and more distinct pattern in non-linear classifier, elucidating that ATM unmutated samples are more likely to be biallelic for ATM and AGO2.

Figure 7 .
Figure 7. a) Example of a rpart decision tree.Chronic lymphocytic leukemia (CLL) samples are stratified according to TCL1A protein status.Model design then includes IGHV gene mutation status, mRNA markers linked to adverse prognosis (from algorithm described in Figure 5c), and further clinical features.IGHV gene mutations status, as seen in top of branch, is the most informative divider.When left out in b) an mRNA marker linked to adverse prognosis (with somewhat arbitrary cut-off for illustrative purposes) functions as the most informative divider.c-e) ctree offers more intuitive visualizations of decision trees.c) When stratifying CLL samples by TCL1A mRNA expression, IGHV mutations status is the most informative divider.d) This is confirmed when stratifying CLL samples by IGHV mutations status (switching the comparison) hence TCL1A mRNA expression is the most informative discriminator.e) T-cell prolymphocytic leukemia (T-PLL) samples stratified by ATM mutation status.Co-variates include ATM deletion, miR-34B deletion, MYC amplification, AGO2

Table 2
residual unambiguous probe sets assigned to a gene are then further summarized by calculation of average expression values per gene.AIMS Medical ScienceVolume 3, Issue 3, 248-271.