RNA-seq reveals more consistent reference genes for gene expression studies in human non-melanoma skin cancers

Patient samples

Skin lesions and non-photodamaged skin tissue samples were collected from patients at the Dermatology Department at the Princess Alexandra Hospital. The study was approved by Metro South Human Research Ethics Committee and The University of Queensland Human Research Ethics Committee (HREC-11-QPAH-236, HREC-11-QPAH-477, HREC-12-QPAH-217, and HREC-12-QPAH-25). Written, informed consent was obtained from all patients prior to participation. Following biopsy, tissues were immersed in RNA later (Life Technologies, Carlsbad, CA) and stored at −80 °C until required. All samples were sectioned and processed according to routine protocol in the Department of Anatomical Pathology located in Princess Alexandra Hospital.

RNA isolation and cDNA synthesis

RNA isolation was performed using the Qiagen RNeasy Plus Mini kit (Qiagen GmbH, Hilden, Germany). Briefly, tissue samples were cut into small pieces, and transferred into 1.5 mL tube containing lysing matrix D (MP Biomedicals, Santa Ana, CA, USA) and 600 uL buffer RLT containing 1% beta-mercaptoethanol and homogenised using a Fast Prep benchtop homogeniser (MP Biomedicals, Santa Ana, CA, USA). Samples were spun five times using setting 6.5 for 30 s each, and chilled on ice between spins. Lysate was removed and transferred to a fresh tube. For unfixed BCC samples embedded in OCT, 20 × 10 µm sections were cut and placed into 600 µL buffer RLT, and homogenised using a 18.5 gauge blunt needle attached to an RNAse free syringe and resuspended five times. The remaining RNA extraction steps were performed as above. RNA concentration was measured using the Qubit fluorometer (Life Technologies, Carlsbad, CA, USA) and RNA integrity determined using the 2100 Bioanalyser (Agilent Technologies, Palo Alto, CA, USA) on RNA Pico chips. The minimum acceptable quality for RNA for analysis by qPCR was RIN >6. Complementary DNA (cDNA) was synthesised from 200 ng total RNA using the Sensifast cDNA kit (Bioline, London, UK) as per the manufacturer’s instructions.

Library preparation and RNA-seq

RNA-seq libraries of poly (A) RNA from 500 ng total RNA obtained from AK, IEC and SCC and SK samples, were generated using the TruSeq unstranded mRNA library prep KIT for Illumina multiplexed sequencing (Illumina, San Diego, CA, USA). TruSeq stranded mRNA library preparation kit was used to generate poly (A) RNA libraries for RNA obtained from normal skin samples (Illumina, San Diego, CA, USA). Libraries were sequenced (100 bp, paired-end) on the Illumina HiSeq 2500 platform and FASTQ files were analysed.

Bioinformatics pipeline

Sequencing data were checked for sequencing quality by FASTQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Adaptors and poor quality sequences were then removed using Trim Galore v0.3.7 (∼6% of reads removed) (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/). An average of 39.1 million trimmed reads pairs were then aligned using Tophat (V 2.1), using the unstranded protocol, against the Human Genome (hg19 build) guided with a human transcriptome generated from the GENCODE Gene annotation version 19 (Table S1) (Harrow et al., 2012; Kim et al., 2013). Quantification of gene expression based on counting the overlaps of mapped reads with genes annotated in the GENCODE gene annotation v19 using HTSeq (Version 0.6.1p2) (Anders, Pyl & Huber, 2015; Harrow et al., 2012) (Table S2). Read counts were then normalised using trimmed mean normalization method (TMM) using bioconductor’s edgeR package (Robinson, McCarthy & Smyth, 2010; Robinson & Oshlack, 2010). Based on the Gencode gene model, an expression level of 57,278 genes was counted. Read counts for the remaining 57,278 genes were then transformed into Transcripts Per Kilobase Million (TPM) values as previously described (Wagner, Kin & Lynch, 2012) (Gene length information found in Table S3).

To ensure that RGs can be detected in any of the samples, genes with any samples without detectable expression level (counts < 0) were discarded (38,670 genes discarded).

Statistical analysis for identification of RGs from RNA-seq

The mean TPM value of a gene is calculated by taking the average of the gene expression across all samples. The CoV was measured by taking the standard deviation of expression value (TPM) for a given gene divided by its mean. The ratio of the maximum to minimum (MFC) was calculated using the largest TPM value divided by the smallest TPM value. A product score (MFC-CoV) was calculated based on the multiplication of CoV and MFC value for each gene. Biological Coefficient of Variance (BCV) was obtained by calculating the square root of the tag wise dispersion measured using edgeR bioconductor package (Robinson, McCarthy & Smyth, 2010) normalized counts per million. To facilitate exploration of the statistical analysis, data visualisation of the analysis can be found on http://skinref-dev.dingerlab.org/ and the raw data on ArrayExpress (https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-5678/).

qPCR

The primers for qPCR reactions were designed using NCBI Primer BLAST (http://www.ncbi.nlm.nih.gov/tools/primer-blast/) (refer to Table 1). Primers were designed to span intron boundaries to avoid amplification of genomic DNA and to amplify all isoforms known to each gene based on the NCBI Reference Sequence Database (RefSeq). Primers were synthesized by Sigma-Aldrich (Castle Hill, Australia). qPCR reactions were performed in triplicate using 1 µL diluted cDNA template in a 10 µL total volume. Reactions were performed in 384-well plate format on the ABI Viia7 Real-Time PCR system (Life Technologies, Carlsbad, CA, USA) using Sensifast SYBR Lo-rox mastermix (Bioline, London, UK). A 2-step cycling protocol was performed, comprising an initial 95-degree polymerase activation for 2 min, followed by 40 cycles of 95 degrees for 5 s, then 60 degrees for 20 s. The comparative Ct (ΔΔCt) method was used for data normalisation.

Measurement of novel RG stability

The RG stability was calculated using the geNORM algorithm (Andersen, Jensen & Ørntoft, 2004), which is integrated with qbase+ software (Biogazelle, Gent, Belgium) and Normfinder software (Vandesompele et al., 2002) (available as an add-in for Microsoft Word at http://moma.dk/normfinder-software).

Statistical analysis of qPCR data

Statistical tests used have been described in each figure legend. Data analysis was performed using GraphPad Prism version 5.04 for Windows (GraphPad Software, Inc., La Jolla, CA, USA).

Identification of novel candidate RGs

To identify RG specific for NMSC and precancerous lesions, we first performed gene expression profiling using RNA-seq on 13 AK, seven IEC, five SCC lesions and four non-photo damaged skin samples. As previously described, a RG should show similar expression across samples, expressed at detectable levels and not display any exceptional expression in any of the samples (De Jonge et al., 2007; Eisenberg & Levanon, 2013). As with previous studies on identification of housekeeping genes, to identify genes that fall within these criteria, we measured the mean expression, CoV and the MFC for each gene within the dataset (De Jonge et al., 2007; Eisenberg & Levanon, 2013) (Table S4). The CoV measures variability within samples and the MFC is the ratio of the maximum and minimum TPM values for a given gene.

Ideally, a RG candidate should have a low CoV and MFC value and is expressed at detectable levels. As such, we generated a score based on the product of CoV and MFC value for each gene (MFC-CoV). We then obtained an initial list of 4643 gene candidates with a MFC-CoV that falls below the lower quartile value of MFC-CoV ≤1.63084, (Table S7). KEGG Pathway enrichment analysis on this initial gene list using geneSCF (Subhash & Kanduri, 2016) revealed that most of these stably expressed genes were involved in maintaining critical cellular activities including the citrate cycle, proteasome, RNA polymerase, protein export, RNA transport and Ribosome activities (B.H FDR Value < 1.8e–30) (Huang da, Sherman & Lempicki, 2009) (Fig. 1C and Tables S5–S6).

To identify RGs specific for NMSC and precancerous lesions, we shortlisted 10 candidate genes from the enriched pathways involved in critical cellular activity for further validation with qPCR. In addition, we selected candidates with functions well described in literature. For instance, we chose the RPLP0 gene, whose function is not only well known for different cell and tissue types but also shown to be a suitable RG for research in the differentiation of human epidermal keratinocytes (Bar, Bar & Lehmann, 2009). Finally, to distinguish mRNA from genomic DNA, we selected multi-exonic genes as candidates to aid design of primers across intron boundaries.

To determine the performance of our gene candidates, we compared the CoV, MFC and mean expression value of our 10 RG candidates with three commonly used RGs in qPCR analysis of skin—ACTB, GAPDH and HPRT1 (Table 2) (De Kok et al., 2005). With the exception of RPL23, our candidate RGs had a lower CoV-MFC compared to the traditional RGs ACTB, HPRT1 and GAPDH (Table 2 and Figs. 1A and 1B). RPL23 had a lower MFC but a higher CoV than GAPDH. We also calculated the BCV and found that our candidate genes fall below the mean (BCV is similar to that of the traditional house keeping genes (Fig. S1B and Table S2).

Data visualisation of TPM, CoV and MFC metrics of the calculated genes in this study have been made available online (http://skinref-dev.dingerlab.org/) as a resource to facilitate other investigators in the community to explore the datasets and identify their candidate RG of interest (Fig. S2).

qPCR validation of new RGs

To validate and extend our findings from the RNA-seq analysis, we conducted qPCR on our 10 candidate RGs in addition to the three commonly used skin RGs ACTB, GAPDH, and HPRT1 on samples derived from a diversity of skin conditions within the spectrum of NMSC (SCC and precursors and BCC) (refer to Table 1 for list of gene primers). A total of 24 samples were tested and they include AK (4), SCC (3), SK (3), BCC (4), IEC (5), and non-photodamaged skin (5). The use of independent samples across a variety of NMSC skin lesions also minimises any patient specific expression biases of our RGs. This was important as some of the lesions used in RNA-seq was derived from the same patient and readings of MFC, CoV and BCV among our samples could exhibit patient-specific biases. Results from the qPCR were analysed using geNorm (Andersen, Jensen & Ørntoft, 2004) within Qbase+ software (Biogazelle) and Normfinder to determine the consistency of expression values among the samples for each candidate RG (Vandesompele et al., 2002).

Statistically, geNorm conduct pairwise variation (V) analysis to identify genes with the least variance between samples and is denoted as ‘stability’ (M) values. In general, lower M values indicate lower variance in expression value among samples and genes with M values ≤0.5 are associated with homogeneous samples. Remarkably, all of our 10 RG candidates had M values ≤0.5 with RPLP0, RPL7A, RPL23, RPS27A and RPL38 ranked in the top five genes for geNorm M value/Stability value. In addition, to eliminate errors related to the usage of a single housekeeping gene, it is common practice to use two or more housekeeping genes. By calculating the normalization factor based on the geometric mean of multiple control genes, we identified that we need only two of our RG candidates for accurate normalization (geNorm V, V2∕3 = 0.084). V values of <0.15 indicate acceptable stability of the RG combination, indicating no further need for additional RGs. Amongst our RG candidates, the pair of genes with optimal normalization factor was RPL38 and RPS27A, which demonstrated the lowest M values (0.257 and 0.265 respectively) (Fig. 2A).

In addition, Normfinder analysis was performed for the same dataset. Normfinder analysis performs estimation of both intra- and intergroup expression variation for each subgroup of samples (lesion types), with output given as a Stability Value. The most stable candidate was RPL7A, and the best combination of genes was RPL7A and RPLP0 (Fig. 2B). Overall trends between geNorm and Normfinder analyses were similar. In both formats, the traditional RGs ACTB, GAPDH and HPRT1 were ranked as having the most variability in gene expression across the groups (increased stability values), and the genes RPLP0, RPL7A, RPL23, RPS27A and RPL38 ranked in the top five genes for geNorm M value/Stability value.

To demonstrate the significance of our findings in NMSC research, we investigated the difference in expression of keratin 17 (KRT17) in AK between normalization using our candidate RGs and normalization with GAPHD (Fig. 3). When using GAPDH as calibrator, there was an approximate 2-fold increase in levels of KTR17 when comparing non-photodamaged skin to AK (Fig. 3A) or approximate 3-fold increase in SCC (Fig. 3B). However, this fold change was significantly higher at approximately 7-fold for AK and approximately 12-fold for SCC when using the combination of RPL32 and RPS27A or either one of them as RG (P < 0.05). There was no statistical difference between data normalised with RPLP0, RPS7A or both of the two candidate genes. This result demonstrates that use of a RG, which is not stably expressed, can lead to inaccurate data, particularly in instances where the relative fold change is subtle.