The role and robustness of the Gini coefficient as an unbiased tool for the selection of Gini genes for normalising expression profiling data

We recently introduced the Gini coefficient (GC) for assessing the expression variation of a particular gene in a dataset, as a means of selecting improved reference genes over the cohort (‘housekeeping genes’) typically used for normalisation in expression profiling studies. Those genes (transcripts) that we determined to be useable as reference genes differed greatly from previous suggestions based on hypothesis-driven approaches. A limitation of this initial study is that a single (albeit large) dataset was employed for both tissues and cell lines. We here extend this analysis to encompass seven other large datasets. Although their absolute values differ a little, the Gini values and median expression levels of the various genes are well correlated with each other between the various cell line datasets, implying that our original choice of the more ubiquitously expressed low-Gini-coefficient genes was indeed sound. In tissues, the Gini values and median expression levels of genes showed a greater variation, with the GC of genes changing with the number and types of tissues in the data sets. In all data sets, regardless of whether this was derived from tissues or cell lines, we also show that the GC is a robust measure of gene expression stability. Using the GC as a measure of expression stability we illustrate its utility to find tissue- and cell line-optimised housekeeping genes without any prior bias, that again include only a small number of previously reported housekeeping genes. We also independently confirmed this experimentally using RT-qPCR with 40 candidate GC genes in a panel of 10 cell lines. These were termed the Gini Genes. In many cases, the variation in the expression levels of classical reference genes is really quite huge (e.g. 44 fold for GAPDH in one data set), suggesting that the cure (of using them as normalising genes) may in some cases be worse than the disease (of not doing so). We recommend the present data-driven approach for the selection of reference genes by using the easy-to-calculate and robust GC.

abscissa in increasing order of the size of their contribution, and the cumulative contribution is plotted against this on the ordinate. The GC is given by the fractional area mapped out by the resulting 'Lorenz' curve ( Fig. 1). For a purely 'socialist' system in which all contributions are equal (GC = 0), the curve joins the normalised 0,0 and 1,1 axes, while for a complete 'autocracy' , in which the resource or expression is held or manifest by only a single individual (GC = 1), the 'curve' follows the two axes (0,0 → 1,0 → 1,1).
Since the early origins of large-scale nucleic acid expression profiling, especially those using microarrays [26][27][28] , it has been clear that expression profiling methods are susceptible to a variety of more or less systematic artefacts within an experiment, whose resolution would require or benefit from some kind of normalisation (e.g. [29][30][31][32][33][34][35][36][37][38][39]. By this ('normalisation of the first kind'), and what is typically done, we mean the smoothing out of genuine artefacts within an array or a run, that occur simply due to differences in temperature or melting temperature or dye binding or hybridisation and cross-hybridisation efficiency (and so on) across the surface of the array. This process can in principle use reference genes, but usually exploits smoothing methods that normalise geographically local subsets of the genes to a presumed distribution.
Even after this is done, there is a second level of normalisation, that between chips or experiments, that is usually done separately, not least because it is typically much larger and more systematic, especially because of variations in the total amount of material in the sample analysed or of the overall sensitivity of the detector (much as is true of the within-run versus between-run variations observed in mass spectrometry experiments 40,41 ). This kind of normalising always requires 'reference' genes whose expression varies as little as possible in response to any changes in experimental conditions. The same is true for expression profiling as performed by qPCR [42][43][44][45][46][47] , where the situation is more acute regarding the choice of reference genes since primers must be selected for these a priori. Commonly, the geometric mean of the expression levels of that or those that vary the least is selected as the 'reference' . The question then arises as to which are the premium 'reference' genes to choose.
It became obvious that an analysis of the GC of the various genes was actually precisely what was required to assess those 'housekeeping' (or any other) genes that varied least across a set of expression profiles, and we found 35 transcripts for which the GC was 0.15 or below when assessing 56 mammalian cell lines taken from a wide variety of tissues 1 . These we refer to as the 'Gini genes' . Most of these were 'novel' as they had never previously been considered as reference genes, and we noted that their Gini indices were significantly smaller (they were more stably expressed) than were those of the more commonly used reference genes 66 . However, this analysis was done on only one (albeit large) dataset of gene expression profiles. While some of the compilations (e.g. 65,80 ) contain massive amounts of expression profiling data, many of these, especially the older ones, may well be of uncertain quality. Thus, especially since the GC is very prone to being raised by small numbers of large outliers, we decided for present purposes that we should compare our analyses of candidate Gini genes using a smaller but carefully chosen set of expression profiling experiments. The more modern RNA-seq (e.g. [81][82][83][84][85], in which individual transcripts are simply counted digitally via direct sequencing, is seen as considerably more robust 81,86,87 and sensitive 88,89 , and so we selected additional large and recent datasets that used RNA-seq in cell lines and tissues (Table 1). We note too that the precision of these digital methods (as with other, digital, single-molecule strategies [90][91][92], means that the requirement for reasonably high-level expression levels is much less acute. In a similar vein (Table 2), we selected a small number of reasonably detailed studies in which particular housekeeping genes had been proposed as reference genes.
To our knowledge, there are no large-scale studies to determine housekeeping genes in large, cell-line cohorts; the present paper serves to provide one. In addition, we include an experimental RT-qPCR analysis of a subset of the Gini genes.

Results
The Gini Coefficient as a robust measure of gene expression stability in multiple cell-line data sets. We previously identified a number of genes in the Human Protein Atlas (HPA) cell line data set 93 with very low expression variability and thus potential for use as reference genes 1 . However, we did not compare these Gini genes to other genes that have previously been proposed as housekeeping genes. We therefore performed a similar analysis using the potential housekeeping genes we proposed in 1 as well as other reference genes proposed in other studies (Table 2) with additional large RNA-Seq cell line data sets (Table 1). Figure 2A shows a plot of the GC of a variety of candidate Gini genes versus their median expression level in the HPA cell lines dataset set 93 . It is clear that genes we identified previously have much lower GC values in the HPA dataset than do any of the others (just two, VPS29 and CHMP2A, were also identified by Eisenberg and Levanon and another, RPL41, by Caracausi). This is not at the expense of an unusually low expression ( Fig. 2A), a finding broadly confirmed when we look at the median expression levels for the CCLE dataset ( Fig. 2B) and of the Klijn dataset (Fig. 2C). Figure 3 shows the GC values for the various genes in two other datasets, viz CCLE and Klijn. Our previous Gini genes have a lower GC than that of any of the other housekeeping genes in 25 out of 38 cases in Klijn (all under 0.2) and in 26 out of 40 cases for CCLE (all under 0.22). In confirmation of this, and of the correlation found above between the median expression levels in CCLE and Klijn, the GC values are also well correlated with each other for the two datasets (Fig. 3). Thus, although the absolute numbers are slightly larger than are those for the HPA dataset (unsurprisingly, given the much larger number of examples), the trend is still very clear: the GiniGenes remain the best among those variously proposed as reference genes in a variety of large and quite independent datasets. It also suggests that variations in the total amount of mRNA are not an issue either.
Another common statistical measure, more resistant to individual outliers, is the interquartile ratio (the ratio between the 25 th and 75 th percentile when expression levels are ranked); by this measure too, the Gini genes that we uncovered previously stand out as being the least varying ( Fig. 4 A,B). This suggests that, as a measure of gene expression stability, the GC is robust: the GiniGenes have the lowest ratio between their maximum and minimum expression values in the HPA dataset (Fig. 4C) and also the lowest interquartile ratio in their levels of expression in all three cell line data sets explored here (Fig. 4B,C) with good correlation between these two datasets.
Use of the Gini Coefficient to find GiniGenes in an unbiased manner in cell-line data sets. Up to now, our analyses of these data sets have used a set of predefined genes to look at expression stability. We next sought to investigate whether the GC would highlight genes with high expression stability that have been reported by others or by ourselves when performing this analysis in a data-driven manner. To that end, we found 115 genes shared between the three data sets with a GC ≤ 0. 2 (Figs. 5,6). This value for the GC was chosen since reducing this to ≤0.15 meant no or very few genes were found in some data sets (e.g. no genes in the CCLE data set had a GC ≤ 0.15) and going above this meant the number of genes were unmanageable (e.g. 1051 genes with a GC ≤ 0.21 in the Klijn data set). Of the 115 genes shared between the datasets with GC < 0.2, 13 were GiniGenes

Study short name Comments Reference
GiniGene Study presenting novel potential housekeeping genes in cells and tissues from the HPA project cell and tissue RNA-seq data. 1 geNorm or Vandesompele Classic set of reference genes in tissues and a means of analysing them 66 Eisenberg Very detailed analysis of housekeeping/ reference genes in tissues using the Illumina Body Map study of RNA-seq of 16 Human Tissues. E-MTAB-513. 49 Lee Two novel reference genes from a detailed analysis of 281 normal tissue samples from 17 different organs then compares between disease states m and cell lines. 131 Caracausi 646 expression profile data sets from 54 different human tissues. and two were housekeeping genes defined by Caracausi and colleagues (Fig. 5B). When we selected the top 20 expressing genes in each data set, only 13 of these were common across these data sets; Table 3 shows some descriptive statistics of 13 of these, with descriptive statistics of all 115 genes found in Supplementary Table S1  share important roles in RNA transcription, translation and stability [94][95][96][97][98][99][100][101][102] , are implicated in a number of diseases, including cancer 94,97,[103][104][105][106][107][108][109][110][111][112][113] , and some, such as SRSF3 are essential for embryo development 114 . Given their pivotal functions, it may be unsurprising that the expression of these genes are tightly regulated across cell lines of different tissue origins, even where these are cancer cell lines. Overall, the distribution, expression stability and important functional roles of these genes suggest that these are excellent housekeeping genes across different cell types.
Of particular interest to us was finding one gene encoding a mitochondrial phosphate transporter protein (SLC25A3 115 ) to be within this list of the top expressing stably expressed genes. This might seem logical since mitochondrial ATP synthesis is required by all cell types and tissues. Figure 7 shows the robustness of the GC for the subset of 115 genes common between the three data sets studied here with a low GC (<0.2). Lower Gini coefficients correlate with lower IQR and Max:median ratios ( Fig. 7: only results for the Klijn data set are shown). The range of IQR values of these genes was smaller in the larger two data sets (CCLE, 1.42-1.67; Klijn, 1.30-1.64) than in the HPA data set (1.26-1.84) suggesting the measured expression values were more stable in the larger data sets (Supplementary Table S1). This may, however, be due to a larger number of cell lines in these two large datasets (934 and 622 in CCLE and Klijn) compared with the HPA data set (56 cell lines).
Application of the Gini coefficient to human tissue RNA-Seq data sets. The results presented thus far are representative of human cell lines. Most reports in the literature regarding housekeeping genes refer to tissue expression data. This may be due to the cell lines being "dedifferentiated" with respect to the tissues from which they are derived 116 .
In our previous report 1 we also analysed RNA-Seq data from tissues 93 and found 22 genes with a GC < 0.15, of which 3 (CHMP2A, VPS29 and PCBP1) were also found in cell line data with a GC < 0.15. The median expression level and GC of these and other candidate GiniGenes in this tissue data set are shown in Fig. 8. As with cell line data, the genes we previously identified (GGs, green dots in Fig. 8) have much lower GCs in this tissue data set than do any of the other candidate GiniGenes, with only two of these genes (VPS29 and CHMP2A) identified previously by Eisenberg & Levanon 49 . The low GC value of these GiniGenes is not at the expense of low expression: of the 22 GiniGenes, 13 are expressed at a median level of between 40 and 200 TPM (see Supplementary  Table S2). Moreover, the GC was also representative of the variation in expression of these genes (albeit influenced to a lesser extent by outliers), as shown in Fig. 9A,B, with all GiniGenes having a GC < 0.15 and the lowest RSD (relative standard deviation), ranging from 24.096% to 28.66% and IQR (1.26 to 1.44) of this list of housekeeping genes. The expression of other housekeeping genes such as GAPDH, ACTB, RPL13A, SDHA, B2M was quite varied according to these measures. For example, the GC of GAPDH (a commonly used HKG) was 0.33, with a RSD of 72.4% and IQR of 2.24, and for ACTB (another commonly used HKG) these values were 0.29, 55.24%, and 2.11.
The median expression levels of the proposed reference genes show a similar level of correlation between the data sets as was found with the cell line data ( Fig. S1A-C), and GiniGenes displayed a mid-range level of expression. The GC of the tissue GiniGenes we proposed however, tended to be higher and more variable in their GC values than in the HPA dataset (Fig. S2,A-C) suggesting that those genes may be representative of the HPA data set only. As an example, in the GTEx dataset only 28 genes had a GC < 0.2, of which the majority (17) were those reported by Caracausi and colleagues, and 7 were GiniGenes. The results here are likely influenced by the number and status (disease or normal) of the tissues analysed in the various data sets compared; for example, the GTEx data come from 53 different, normal human tissues, whereas the HPA tissue data include a mixture of disease and normal tissue samples. In addition, compared to the cell line data where hundreds (in the case of the Cancer Cell Line Encyclopedia) of cell lines were analysed, the number of tissues in these data sets was fewer than 100. In the case of the data set used by Eisenberg and Levanon 49 , viz. the Illumina Human Body Map (E-MTAB-513), 10 of the 11 housekeeping genes proposed here (which included 2 Gini Genes, CHMP2A and VPS29) had a GC ≤ 0.2 and were reasonably well expressed (with median expression levels between 50-270 TPM, see Supplementary Table S2 and Supplementary Fig. S4). This may be compared to the 5 other GGs with GC < 0.2 in this data set whose expression value was lower, with median expression between 19-35 TPM. This suggests that finding suitable HKGs may be dependent on the data set itself, and the type of tissue under investigation.
We next sought to perform a more comprehensive and integrative analysis by filtering the tissue data sets to only include genes with a GC ≤ 0.2 to find common genes across these data sets with reasonable expression stability (Supplementary Table S3). As shown in Fig. 10 only 15 genes were shared between the four data sets with a GC ≤ 0.2, none of which has been reported previously as a housekeeping genes. Table 4 shows some descriptive statistics of these genes. In any case, the names of the proteins encoded by these 15 genes suggest these play important and essential roles. The median expression values of these genes varied from around 10-450 TPM, with SNX3 www.nature.com/scientificreports www.nature.com/scientificreports/ (Sorting nexin-3 (Protein SDP3)) and COX4I1 (Cytochrome c oxidase subunit 4 isoform 1) being consistently the two highest-expressing genes.
Sorting nexins are a group of cytoplasmic and membrane-associated proteins involved in the regulation of intracellular trafficking 117 . SNX3 has been reported to play a role in receptor recycling and formation of multivesicular bodies 118 , and its dysregulation has been implicated in disorders of iron metabolism and the pathogenesis of some neurodegenerative diseases 119,120 .
The COX4I gene encodes the nuclear-encoded cytochrome c oxidase subunit 4 isoform 1, the terminal enzyme in the mitochondrial respiratory chain. Given the key role of the mitochondrial respiratory chain in all human cells (except red blood cells), stable expression of such a gene in all tissues may not be a surprising result. Increased RNA COX4I1 levels have been reported in sperm of an obese male rat model 121 and thus may play a role in obesity-related fertility problems, and reduced expression of this subunit leads to a reduction in mitochondrial respiration as well as sensitising cells to apoptosis 122 .
The small number of genes shared between these data sets with a GC < 0.2 indicates that the data in these studies are more variable compared to cell lines alone. The cause of this variation may be due to the tissue data having been obtained from different subjects 123 . Moreover, tissues are themselves a mixture of cell types with varying levels of gene expression in each cell type 124 , while cell lines are nominally clonal.
Our results suggest that in the case of RNA-seq tissue data sets, where gene expression tends to be more variable, an unbiased approach, using the Gini coefficient, may be more fruitful in the search for stably expressed genes with which to perform normalisation, than the other commonly used methods used until now 123,125 . RT-qPCR analysis of gene expression stability of some housekeeping genes in 10 cell lines. In order to illustrate the utility of the GC to find suitable housekeeping genes, we next chose to assess this experimentally by RT-qPCR using a small subset of candidate reference genes (40; top 32 genes from genes ordered by GC and expression value from 94 , plus 8 of the most commonly used from the literature, including seven from 66 and one (RPL32) from 126,127 , and 10 cell lines from a range of tissues (see Tables 5 and 6). We first set a Cq value (which is inversely proportional to expression level) cut-off of 32, above which no expression is observed, and subsequently used the Cq values of genes in cell lines as a relative expression level (Cq cut off/Cq value of gene). Descriptive statistics of the expression of each gene in individual cell lines were then calculated. As a final step, the median expression value of each gene in individual cell lines was used to calculate descriptive statistics, including the GC, of gene expression across these cell lines. Figure 11 illustrates a KNIME workflow 128-130 that we wrote for this purpose. The raw data and descriptive statistics extracted are provided in Supplementary Tables S5 and S6 respectively, and the KNIME analysis workflow in Supplementary File 1. Figure 12 uses RT-qPCR data to plot the GC of the candidate reference genes analysed here versus their relative median expression level. Three GiniGenes 94 (RBM45, TRNT1 and CNOT2) had very low and variable expression. Most of the other genes analysed showed low GC values with a range of (relative) expression values; the inset in Fig. 12 shows genes with a GC < 0.2 including a mix of 35 genes: 26 GiniGenes and 6 housekeeping genes referenced by Vandesompele and colleagues 66 , one referenced by Caracausi 65 and one by Lee et al. 131 . Two of these GiniGenes, HNRNPK and PCBP1, which we also found to be stably expressed in the cell line data suggesting these may be potential stable housekeeping genes. As shown in Fig. 13 and inset, the GC is well correlated with the % RSD.
More importantly, the GC of our GiniGenes was particularly low (Fig. 12). The low absolute magnitude reflected the fact that Cq value is based on a logarithmic scale. Various commonly used housekeeping genes (HPRT1, GAPDH, ACTB, SDHA, HMBS and B2M) displayed higher % RSDs and GC than other genes studied here in spite of their higher relative expression levels. This was also the case when inspecting the interquartile ratio against the GC of these (Fig. S3).
The above results suggest that the GC is also applicable to RT-qPCR data, with GiniGenes having good potential (as novel "housekeeping" genes) for the normalisation of such data.  Table 2  www.nature.com/scientificreports www.nature.com/scientificreports/

Discussion
Reference genes are commonly used to normalise gene expression data, so as to account for bias resulting from both biological and technical variability, and to enable quantification of gene expression changes or differences in the system under study. It is generally considered that such reference genes should come from pathways that are required for general metabolism, using only one gene per 'pathway' to avoid co-regulation which might make the gene expressions look very stable.
Various tools have been developed to evaluate and screen reference genes from experimental datasets; these include geNorm 66   www.nature.com/scientificreports www.nature.com/scientificreports/

O00442
Catalyzes the conversion of 3'-phosphate to a 2' ,3'-cyclic phosphodiester at the end of RNA. The mechanism of action of the enzyme occurs in 3 steps: (A) adenylation of the enzyme by ATP; (B) transfer of adenylate to an RNA-N3'P to produce RNA-N3'PP5' A; (C) and attack of the adjacent 2'-hydroxyl on the 3'-phosphorus in the diester linkage to produce the cyclic end product. The biological role of this enzyme is unknown but it is likely to function in some aspects of cellular RNA processing.  www.nature.com/scientificreports www.nature.com/scientificreports/ These tools assess expression stability of genes in different ways: • geNorm determines gene stability through a stepwise exclusion or ranking process followed by averaging the geometric mean of the most stable genes from a chosen set. Python implementation: https://eleven.readthedocs.io/en/latest/. • BestKeeper also uses the geometric mean but using raw data rather than copy numbers. BestKeeper 136 can be used as an Excel-based tool. It can accommodate up to 10 housekeeping genes in up to 100 biological samples. Optimal HKGs are determined by pairwise correlation analysis of all pairs of candidate genes, and the geometric mean of the top-ranking ones. http://www.gene-quantification.info. • NormFinder measures variation, and ranks potential reference genes between study groups. NormFinder 135 has an add-in for Microsoft Excel and is available as an R programme. It recommends analysis of 5-10 candidate genes and at least 8 samples per group. https://moma.dk/normfinder-software. • The comparative ΔCT finder requires no specialist programmes since this involves comparison of comparisons of ΔCTs between pairs of genes to find a set of genes that show least variability. • RefGenes allows one to find genes that are stably expressed across tissue types and experimental conditions based on microarray data, and a comparison of results from geNorm, NormFinder and Best Keeper to find a set of reference genes. However, this is not a free service unless one searches for one gene at a time. Furthermore, the site for this tool is no longer available. Moreover, all these tools require the user to make a prior selection of such HKGs (introducing bias and potential errors) and most are cumbersome to understand and calculate.
We have here shown how via a simple calculation, the GC, we can find potential reference genes, and illustrated its utility in large-scale cell-line, tissue RNA-Seq data sets and RT-qPCR data. The expression of a number of classical HKGs from a number of carefully selected publications do in fact vary much more substantially between large RNA-Seq data sets, both for tissues and cell lines.
Whilst not all studies will involve large data sets such as those we have analysed here; the GC should also be of use for smaller-scale studies to select a subset of genes in a panel of cell lines or tissues relevant to the study in question.
Overall we find that (i) two of these genes, HNRNPK and PCBP1, seemed to be particularly robustly and stably expressed at reasonable levels in all cell lines studied, and (ii) a data-driven strategy based on the GC represents a useful and convenient method for normalisation in gene expression profiling and related studies.

Methods
The datasets used are described and referenced below. The data, in transcripts per million (TPM) units were downloaded from the EBI expression atlas as a .tsv file. As previously 1 , the Gini Index was calculated using the ineq package (Achim Zeileis (2014). ineq: Measuring Inequality, Concentration, and Poverty. R package version 0.2-13. https://CRAN.R-project.org/package=ineq) in R (https://www.R-project.org/). These calculations were incorporated into KNIME via KNIME's R integration R Snippet node. A spreadsheet giving the extracted analyses is provided as Supplementary Tables (Tables S7 and S8 Table 5. Details of human cell lines used for the assessment of expression of candidate reference genes by RT-qPCR. Figure 11. The KNIME workflow described here to calculate descriptive statistics and the Gini coefficient from RT-qPCR data. This workflow can be adapted for use with large RNA-Seq Data sets.   www.nature.com/scientificreports www.nature.com/scientificreports/ Harvesting cells for RnA extraction. Cells from adherent cell lines were harvested by removing growth media and washing twice with 5 mL of pre-warmed phosphate buffered saline (PBS) (Sigma, Cat No. D8537), then incubated in 3 mL of 0.025% trypsin-EDTA solution (Sigma Cat No. T4049) for 2-5 min at 37 °C. At the end of incubation cells were resuspended in 5-7 mL of respective media when cells appeared detached to dilute trypsin treatment. The cell suspension was transferred to 15 mL centrifuge tubes and immediately centrifuged at 300 × g for 5 min. Suspended cell lines were centrifuged directly from cultures in 50 mL centrifuge tubes and washed with PBS as above. The cell pellets were resuspended in 10-15 mL media and cell count and viability was determined using a Nexcellom Cellometer Auto 1000 Cell Viability Counter (Nexcellom Bioscience) set for Trypan Blue membrane exclusion method. Cells with >95% viability were used for downstream total RNA extraction.
RnA extraction. Total RNA was extracted from 2-5 × 10 6 cells using the Qiagen RNeasy Mini Kit (Cat No. 74104) and DNAse treated using Turbo DNA-free kit (Invitrogen, Cat No. AM1907) according to the manufacturer's instructions. Briefly, 1 X DNA buffer was added to the extracted RNA prior to adding 2U (1 µL) of DNAse enzyme. The reaction mixture was incubated at 37 °C for 30 min and inactivated for 2 min at room temperature using DNAse inactivating reagent. The mixture was centrifuged at 10,000 × g for 1.5 min and the RNA from the supernatant was transferred to a clean tube. The RNA concentration was determined using a NanoDrop ® ND-1000 spectrophotometer and further validated using an Agilent 2100 bio-analyser coupled with 2100 Expert software system. Only RNA samples with an RIN (RNA Integrity Number) between 9-10 were selected for cDNA synthesis.
Reverse transcription and cDnA Synthesis. 1 µg of RNA was reverse transcribed into cDNA. Briefly, a 20 µL reaction was setup by adding 1 µL each of oligodT (50 µM, Invitrogen, cat No. 18418020) and dNTP mix (10 mM, Invitrogen, Cat No. 18427-013) followed by adding an appropriate volume for 1 µg of RNA. Nuclease free water (Ambion, Cat No. AM9937) was then added to make the volume up to 13 µL and incubated at 65 °C for 5 min then cooled on ice for 1 min. To initiate transcription 4 µL of 5 X first strand buffer (Invitrogen, Cat No. Validation of gene expression by genorm. A set of candidate reference genes (40; top 32 genes from genes ordered by GC and expression value from 94 , plus 8 of the most commonly used from the literature including seven from 66 ). RNAseq data were selected for validation of stable gene expression using geNorm 66 . First, a typical qPCR protocol was prepared from a master mix for each gene to be tested per cell line in triplicate. This consisted of 10 µL/well made by adding 0.8 µL of nuclease free water (Ambion), 5 µL of LC480 SYBR Green I Master (2 X conc. Roche, Product No. 04887352001), 0.1 µL each of forward and reverse primers (20 µM) (for primer and amplicon sequences see Supplementary Table S9) and 4 µL of 1:100 diluted cDNA in a 384 well qPCR plate (Starlab Cat. No. E1042-9909-C). The no template controls (NTC) for each gene were produced by replacing cDNA with 4 µL of nuclease free water. Thermal cycling conditions used were: one cycle of 95 °C for 10 min followed by 40 cycles of 95 °C for 10 sec and 60 °C for 30 sec. qPCR was performed using Roche LightCycler LC480 qPCR platform. The fluorescence signals were measured in real time during amplification cycle (Cq) and also during temperature transition for melt curve analysis.
The mean Cq values were converted into relative values for a gene across all cell lines using ΔCq method 138 . Briefly, the lowest Cq value in a panel of cell lines for a gene was subtracted from all the values in that panel using q control , where C q sample is the mean Cq value obtained for a gene in each of the cell lines and C q control is the lowest Cq value in that panel. The relative values for each gene in a panel were then obtained by applying = −Δ R 2 C q . These relative values were applied in geNorm Visual Basic applet for Microsoft Excel ® 66 that determines the most stable reference genes from a set of genes in a given panel of cell lines.

Validation of gene expression using the Gini coefficient.
To the raw RT-qPCR data a Cq value (which is inversely proportional to expression level) cut-off of 32 was set, above which no expression is observed. The Cq values of genes in cell lines were subsequently converted to a relative expression level (Cq cut off/Cq value of gene). Descriptive statistics of the expression of each gene in individual cell lines were then calculated. As a final step, the median expression value of each gene in individual cell lines was used to calculate descriptive statistics, including the GC, of gene expression across these cell lines. Figure 11 illustrates a KNIME workflow [128][129][130] for this purpose. The raw data and descriptive statistics extracted are provided in Supplementary Tables S5 and S6 respectively, and the KMNIME analysis workflow in Supplementary File 1.

Data availability
All data generated or analysed during this study are included in this published article (and its Supplementary  Information Files). The original datasets used are referenced throughout and are summarised in Table 2 www.nature.com/scientificreports www.nature.com/scientificreports/