A novel approach for human whole transcriptome analysis based on absolute gene expression of microarray data

Ethics statement

Blood tissues were obtained from volunteers recruited from Hospital Infantil de México Federico Gomez. Parents or guardians of each child signed the consent form for participation in the project. The study was approved by the Ethics Committee of Hospital Infantil de México Federico Gomez (4000/536/2012).

RNA extraction, microarray hybridization and data analysis

Total RNA from leukocytes of sixteen children (nine males and seven females, ages 6–10 years) was isolated using TRIzol (Invitrogen, Carlsbad, CA, USA). To avoid DNA carryover, samples were treated with recombinant DNAse I (Roche, Basel, Switzerland). The initial concentration of RNA was 1 µg for each sample and the quality was assessed using the Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). Polyadenylated controls were prepared and RNA was converted to complementary DNA (cDNA). In vitro transcription was performed for generating complementary RNA (cRNA). The cRNA was purified and quantified; cDNA was synthesized again and the remaining cRNA was hydrolyzed with RNase H. The cDNA was purified and quantified subsequently it was fragmented into 40–70 base pairs (bp) fragments. The cDNA fragments were labeled and hybridized to an Affymetrix Gene Chip Human Gene 1.0 ST Array for 16 h at 45 °C. This microarray contains probes for the analysis of about 28,869 genes, with an average of 26 probes of 25 oligonucleotides for each gene. After hybridization, microarrays were washed in the Fluidics Station (Affymetrix Gene Chip Fluidics Station 450 FS; Affymetrix, Santa Clara, CA, USA) following the manufacturer’s instructions and scanned using the System Affymetrix Gene Chip Scanner 3000 7G. This analysis was performed on the Affymetrix platform in the Genotyping and Expression Analysis Unit of the National Institute of Genomic Medicine (INMEGEN). Quality control of the microarray data was done using RMA Express (http://rmaexpress.bmbolstad.com/). The original array data was preprocessed with background adjustment, quartile normalization and PLM summarization method by robust multiarray average (RMA) using RMA Express (Irizarry et al., 2003). The raw expression data is available on the NCBI Gene Expression Omnibus (GEO) under accession number GSE89571.

Selection of the Y chromosome genes for absolute gene expression threshold

The 124 genes of the Y chromosome were obtained using the Genome Browser bioinformatics tool (https://genome.ucsc.edu/). The fluorescence intensities of these genes were analyzed in a set of nine males and seven female microarrays (GSE89571) from Affymetrix Gene Chip Human Gene 1.0 ST Array, which were previously normalized by RMA (Irizarry et al., 2003). The non-parametric Mann Whitney test was performed to select the Y chromosome genes with fluorescence significantly different between the sexes. To calculate the absolute gene expression threshold, we used the Z-test to the Z-score at 95% confidence interval.

RNA extraction and Real-time PCR

To perform the RT-qPCR expression analysis, 1 µg of total RNA was reverse-transcribed with TaqMan Reverse Transcription Reagents (Applied Biosystems, Foster City, CA, USA) using random hexamers according to the manufacturer’s protocol. One-tenth volume of each RT reaction was amplified with SYBR Green master mix (Roche, Basel, Switzerland) containing 0.5 mM of intron-spanning specific primers as well as Fast Start enzyme, PCR buffer and 3.5 mM MgCl₂ (Lipshutz et al., 1999), in a final volume of 10 µl, the reactions were measured in a Light Cycler 2.0 real-time PCR detection system (Roche, Basel, Switzerland). Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) expression was measured as a reference. PCR was conducted using the following cycling conditions: pre-incubation and denaturation at 95 °C during 10 min. Amplification with 35 or 40 cycles that included: denaturation at 95 °C for 10 s with a thermal ramp rate at 20 °C/s; annealing at 61 °C for 7 s with a thermal ramp rate at 20 °C/s; amplification at 72 °C for 10 s with a thermal ramp rate at 20 °C/s using GAPDH primers. The primer design for the selected genes for RT-qPCR analysis was performed using the Primer-BLAST bioinformatics tool (Table S1).

KEGG pathway enrichment

The KEGG enrichment analysis was conducted using WebGestalt (Zhang, Kirov & Snoddy, 2005). To this end, the 5% most expressed genes shared between the 16 arrays were selected to identify the KEGG pathways significantly enriched. Thus, the pathway enrichment was conducted using 1,539 genes. All human genes were used as background to calculate statistical significance and setting the hypergeometric test for enrichment evaluation analysis. The p values were adjusted using the Benjamini–Hochberg multiple test adjustment. Finally, we identified the enriched KEGG pathways with an adjusted P-value (p ≤ e⁻⁹).

Selection of the Y chromosome genes to establish the absolute gene expression threshold

The Y chromosome is the sex-determining genetic element in higher mammals, and it is absent in female subjects. Therefore, the fluorescence of the Y chromosome genes in female subjects could be used as the background fluorescence for all the probes in the microarray. To test this hypothesis, the transcriptome of leukocytes from 16 human samples were analyzed. After microarray data normalization, we obtained the fluorescence intensities for the 124 genes of the Y chromosome (Table S2). We found several genes with similar fluorescence intensity between male and female subjects (gene A and gene C in Fig. 1A), while other genes had a sex-dependent differential fluorescence (gene B in Fig. 1A). The histogram with the fluorescence of the 124 genes did not show a differential fluorescence distribution between male and female subjects (Fig. 1B). To discern the genes with a sex-dependent expression, only the genes with a significant difference (p ≤ 0.05) in fluorescence between the sex groups were selected (Fig. 1C). However, several genes showed an overlapped fluorescence between male and female (gridded region in Fig. 1C). Thus, we eliminated the genes which fluorescence distribution was overlapped between sexes. After that, only three genes named DDX3Y (probe: 8176624), EIF1AY (probe: 8176719), and TXLNG2P (probes: 8176709 and 8176698) showed a fluorescence significantly different between sexes (p = 4.10⁻³²) (Fig. 1D, Table S3).

Thus, the fluorescence emitted by these three genes in female subjects (Fig. 1E) was used to calculate the upper limit of fluorescence with a 95% confidence. The value found was 4 and it was selected as the absolute gene expression threshold. Hence, the genes with fluorescence values ≥4 were considered as expressed, while the genes with fluorescence <4 were considered as non-expressed genes (background fluorescence).

We selected 11 housekeeping genes from the Cheng-Wei (2011) study (Chang et al., 2011) and their fluorescence intensity (Fi) was analyzed in our 16 microarray samples. These genes have high (>10), medium (6–10)or low (4–6)fluorescence intensity (Fig. 2, Table S4) and according to the absolute gene expression threshold (≥4), the 11 housekeeping genes were expressed in all samples (green lines in Fig. 2). Importantly, their fluorescence was not sex-dependent (p = 0.81) (blue and pink zones in Fig. 2). This result was expected because housekeeping genes are typically expressed under normal conditions. Also, there were housekeeping genes with high and low fluorescence values as compared to the fluorescence of the absolute gene expression threshold (red dotted line in Fig. 2).

The SRY gene was analyzed as a negative control of the gene expression. We found that SRY fluorescence was under the absolute gene expression threshold, indicating that this gene was not expressed in either sample (blue line in Fig. 2). Furthermore, the fluorescence of the SRY gene was not significantly different between sexes (p = 0.93) (blue line in Fig. 2). This result was expected because it has been reported that SRY gene, unique for the Y chromosome, is only expressed during the human embryonic development in male subjects (Berta et al., 1990).

Validation of the absolute gene expression threshold by RT-qPCR

The consistency of the absolute gene expression threshold in assessing the expression of a gene was also evaluated using RT-qPCR. To this end, we selected ten genes that exhibited the lowest variability of fluorescence intensity in all the microarrays, and their expression levels were measured by RT-qPCR. These genes were separated into three groups according to their fluorescence level (Fig. 3, Table S5). The first group contained six expressed genes with fluorescence intensity above the absolute gene expression threshold (Fi ≥ 4). This group included two genes with high fluorescence intensity (Fi ≥ 10, red lines in Fig. 3), two genes with medium fluorescence intensity (Fi 6–10, orange lines in Fig. 3) and two genes with low fluorescence intensity (Fi 4–6, green lines in Fig. 3). The second group contained a set of non-expressed genes (fluorescence intensities below the absolute gene expression threshold of 4). This second group included two genes with fluorescence close to the expression threshold (3–4, dark blue lines in Fig. 3) and two genes with the lowest fluorescence (Fi < 3, light blue in Fig. 3). The third group contained the three genes used to define the absolute gene expression threshold (DDX3Y, EIF1AY and TXLNG2P) by RT-qPCR. Importantly, the expression level observed between RT-qPCR (measured by Ct: Ct threshold cycle) and microarray (measured by fluorescence intensity) had a good Spearman correlation (r = 0.882, p = 2.84⁻⁶⁹) on the 13 genes (Fig. S1). The fluorescence of the genes DDX3Y, EIF1AY and TXLNG2P shown a significant difference between male and female subjects (p = 6.74⁻²⁰) (Fig. 4A, Table S6), while, the fluorescence level of the other 10 genes was not significantly different between sexes (p = 0.98) (Table S7).

According to the Ct and fluorescence levels (Fig. 4B), three clusters of genes were observed: one cluster containing highly expressed genes (cluster 1), another cluster containing medium expressed genes (cluster 2) and a third cluster containing low and the non-expressed genes (cluster 3). These results suggest that the absolute gene expression threshold not only assessed what genes are expressed, but also allowed the level of expression according to the fluorescence intensity. However, the genes BEND7 and CAPN11 have low expression and their fluorescence values overlapped with the non-expressed genes, suggesting that genes with low fluorescence (close to the absolute gene expression threshold) can be really not expressed genes (cluster 3 in Fig. 4B).

Human transcriptome profiling using the most expressed genes

We selected the 5% of genes with the highest fluorescence intensity (most expressed genes) shared between the 16 microarrays to identify the KEGG pathways significantly enriched, resulting in 1,539 genes that were expressed in all the leukocyte samples. The pathway enrichment analysis showed 45 significantly enriched pathways (p ≤ E−9) (Table S9). Interestingly, the most enriched pathways were those expected for a leukocyte cell type (Table 1). For example, Phagosome (p = 5.67E−35), chemokine signaling pathway (p = 1.62E−21), Fc gamma R-mediated phagocytosis (p = 6.49E−23), Natural Killer cell mediated cytotoxicity (p = 1.5E−19), antigen processing and presentation (p = 9.22E−24), T cell receptor signaling pathway (p = 9.08E−27), B cell receptor signaling pathway (p = 5.77E−14) and leukocyte transendothelial migration (p = 5.74E−17).

Application of the absolute gene expression threshold to another microarray platform

We calculated the absolute gene expression threshold for previously published microarray data of liver cells (GSE18269) using our approach (Fig. 1). After that, we selected the 5% most expressed genes shared between all samples, resulting in a set of 1,547 expressed genes which were used to identify the significantly enriched pathways. We found that 37 KEGG pathways (p ≤ E−9) were significantly enriched (Table S10). The majority of these enriched pathways were functional related with liver cells (Table 2). For example, Complement and coagulation cascades (p = 9.19E−37), Drug metabolism-cytochrome P450 (p = 3.5E−24), Fatty acid metabolism (6.2E−22), Metabolism of xenobiotics by cytochrome P450 (8.03E−21) and Glycolysis/Gluconeogenesis (1.8E−19).

Additionally, we calculated the absolute gene expression threshold for another microarray platform. To this end, we selected the previously reported GSE22974 data in which yeast RNA was hybridized to a human microarray (U133 plus 2.0) and the absolute gene expression threshold was calculated for this dataset. After that, the fluorescence value for the absolute gene expression threshold was 9.17. Thus, all the probes with fluorescence above this fluorescence value were considered as expressed. According to our expression threshold, the 97.9% (n = 53,466) of the probes were not expressed (Fig. S2). This result was expected, because the hybridized RNA was from yeast and the probes on the microarray were designed for human genes.

The Gene Expression Barcode algorithm provides absolute measures of expression for human genes using an expression threshold for each probe according to the probe effect of each one (McCall et al., 2013). We applied this algorithm on the same data (GSE22974) to compare their results with the obtained using our expression threshold. The Barcode algorithm found that 94% (n = 51,196) of probes also were not expressed. Importantly, the 93.7% (50,638) of these non-expressed probes were shared between Barcode and our approach (Fig. S3). These results suggest a good consistency to call expressed genes between the Barcode algorithm and the threshold approach reported here.