A standardized fold change method for microarray differential expression analysis used to reveal genes involved in acute rejection in murine allograft models

Murine transplantation models are used extensively to research immunological rejection and tolerance. Here we studied both murine heart and liver allograft models using microarray technology. We had difficulty in identifying genes related to acute rejections expressed in both heart and liver transplantation models using two standard methodologies: Student's t test and linear models for microarray data (Limma). Here we describe a new method, standardized fold change (SFC), for differential analysis of microarray data. We estimated the performance of SFC, the t test and Limma by generating simulated microarray data 100 times. SFC performed better than the t test and showed a higher sensitivity than Limma where there is a larger value for fold change of expression. SFC gave better reproducibility than Limma and the t test with real experimental data from the MicroArray Quality Control platform and expression data from a mouse cardiac allograft. Eventually, a group of significant overlapping genes was detected by SFC in the expression data of mouse cardiac and hepatic allografts and further validated with the quantitative RT‐PCR assay. The group included genes for important reactions of transplantation rejection and revealed functional changes of the immune system in both heart and liver of the mouse model. We suggest that SFC can be utilized to stably and effectively detect differential gene expression and to explore microarray data in further studies.

Murine transplantation models are used extensively to research immunological rejection and tolerance. Here we studied both murine heart and liver allograft models using microarray technology. We had difficulty in identifying genes related to acute rejections expressed in both heart and liver transplantation models using two standard methodologies: Student's t test and linear models for microarray data (Limma). Here we describe a new method, standardized fold change (SFC), for differential analysis of microarray data. We estimated the performance of SFC, the t test and Limma by generating simulated microarray data 100 times. SFC performed better than the t test and showed a higher sensitivity than Limma where there is a larger value for fold change of expression. SFC gave better reproducibility than Limma and the t test with real experimental data from the MicroArray Quality Control platform and expression data from a mouse cardiac allograft. Eventually, a group of significant overlapping genes was detected by SFC in the expression data of mouse cardiac and hepatic allografts and further validated with the quantitative RT-PCR assay. The group included genes for important reactions of transplantation rejection and revealed functional changes of the immune system in both heart and liver of the mouse model. We suggest that SFC can be utilized to stably and effectively detect differential gene expression and to explore microarray data in further studies.
At the stage of organ failure, organ transplantation is a life-saving medical procedure, though it still has some problems, e.g. transplant rejection and the requirement for life-long immunosuppressive drugs. Transplantation models without immunosuppression are important and the mechanisms of rejection and tolerance in these models need to be revealed.
Microarray technology is well established and widely used, providing a picture of gene expression or RNA profiling in different tissues [1]. To identify differential expression, Student's t test and linear models for microarray data (Limma) are two popular choices [2][3][4]. The t test utilizes information for all the samples (or standard deviations) in one microarray probe and is conducted independently among different probes [4], while Limma uses the empirical Bayesian approach of shrinkage of the estimated sample variances towards a pooled estimate. The information (means and standard deviations) from all the probes in a replicate set of experiments is combined and used at the level of one probe to detect differential expression in Limma [2].
In the present study, we established murine heart and liver allograft models and used microarray technology to reveal the significant genes that related to transplant rejection. By using the t test and Limma, no significant intersecting genes were obtained in these models. Therefore, we developed a new method, named standardized fold change (SFC), to detect differential expression by taking information from the neighbors of one probe with an adjustable bin size. To compare SFC with the t test and Limma, we generated a simulated data set to estimate the performance and used the real experimental datasets from the MicroArray Quality Control (MAQC) platform and the transplantation model to estimate the reproducibility. We concluded that SFC can be applied as a new and effective approach to detect differential expression and contribute more reliable results in microarray studies. Then, SFC reported a set of significant genes from expression data from the murine heart and liver allograft, and we further validated them by qRT-PCR. Gene expression changes revealed functional reactions and pathway activities in the early stage of allograft in both heart and liver.

Transplantation and RNA extraction
Cardiac transplantation was performed from a sex-matched B10 donor to a CBA recipient by microsurgical techniques. Intra-abdominal vascularized heterotopic mouse cardiac transplantation was performed [5]. The cardiac graft survival was determined using daily palpation of the recipient's abdomen. Three case samples on the fifth day were obtained. BR mice were used as donors and D2 mice were used as recipients in the orthotopic hepatic transplantation. We performed transplantation surgery on the mice [6] in which for orthotopic liver transplantation, BR mice were used as donors and D2 mice were used as recipients. We subsequently transplanted the liver into the recipient mice using the cuff technique [6]. Grafts were harvested at postoperative day 5 (POD5) or at POD8 after transplantation and were submerged in RNAlaterÒ stabilization solution (Life Technologies, Carlsbad, CA, USA) for freezing. Total RNA was extracted from frozen tissue samples using ISO-GEN (NipponGene, Tokyo, Japan). We also designed control groups of three normal cardiac tissues and three hepatic tissues.

Standardized fold change method
The probe signals from microarray data were firstly natural log-transformed and then manipulated with quantitative normalization. To assess the differential expressions among cases and controls, the statistic SFC is defined as: For the variance of each probe, we ranked all probes by the mean values of signals from all samples and then took the median value of its b nearest neighbors as the variance, where the default bin size of b here is 1000. The SFC software now implements this algorithm in the Linux system and can be found at https://github.com/WeichenZhou/SFC.

Simulation data study
We generated the simulated data from simple formulas with the Gaussian noise (mean = 0, variance = 1) as a default distribution for gene expression data [7]. The control and case samples in the null hypothesis are shown as follows: where h represented the differential expression underlying cases versus controls and we defined h 0 as 0% and k is 1.
The control and case samples in the alternative hypothesis are shown as follows: We defined h 1 as 10%, 25% and 50%, respectively. The size of real positive calls consists of 1%, 5% and 10% of the whole simulated data, respectively. Following these, a 100-time simulation was conducted to assess the false positive rate (FPR) and the false negative rate (FNR).

MAQC data and the reproducibility analysis
The MAQC project was developed by the US Food and Drug Administration (FDA) to provide standards and quality control metrics and involved six centers [Applied Biosystems (Thermo Fisher Scientific, Waltham, MA, USA), Affymetrix (Santa Clara, CA, USA), Agilent Technologies (Santa Clara, CA, USA), GE Healthcare (Chicago, IL, USA), Illumina (San Diego, CA, USA) and Eppendorf (Hamburg, Germany)] that are major providers of microarray platforms and RNA samples [1,8]. The reproducibility of the top 100 and 1000 significant genes was estimated inter-and intra-platform by the three statistical methods, and heatmaps were drawn with the matrix of each batch. For the expression data from the mouse transplant model, we picked up two out of three cases and controls to build one batch and made a 9 9 9 matrix heatmap to estimate the reproducibility. The significance level of mouse microarray data was 0.05.

Application on mouse transplantation data
We detected differential expression of genes between cases and controls in three phases: POD5 of cardiac transplantation, POD5 of hepatic transplantation and POD8 of hepatic transplantation. All P-values from expression data were adjusted by the Bonferroni correction. After getting all significant probes from SFC, we converted the probe level significance to gene level using an annotation file. Venn diagrams showed the significant genes with differential expression. Pathway and gene ontology (GO) enrichment analyses were performed by using the Database for Annotation, Visualization and Integrated Discovery (DAVID; http://david.abcc.ncifcrf.gov/) with the Bonferroni correction-adjusted P-values < 0.05 [9]. Mouse transplantation data have been deposited in NCBI's Gene Expression Omnibus [10] and are accessible through GEO Series accession no. GSE89340. All data were conducted by quantile normalization before processing by different methods. Limma can be found as the R package LIMMA [2,3] and the heatmaps were created by GPLOTS. All R packages can be downloaded from Bioconductor (www.biocond uctor.org).

Quantitative RT-PCR (qRT-PCR)
The RNA was reverse-transcribed to cDNA using a Prime-ScriptÒ RT Reagent Kit (Takara Bio, Shiga, Japan) as described previously [11]. The sequences used in our study are shown in Table S3. Quantitative RT-PCR (qRT-PCR) was performed using a SYBR Green system on the Applied Biosystems PRISM7700 instrument (Thermo Fisher Scientific), and experiments were conducted using 0.4 lM of each primer in a final reaction volume of 20 lL of KAPA SYBRÒ FAST qPCR kit (Kapa Biosystems, Cape Town, South Africa). The PCR cycling conditions were as follows: 95°C for 30 s, and 50 cycles of 95°C for 5 s, 60°C for 1 min. The normalized threshold cycle (C t ) value of each gene was obtained by subtracting the C t value obtained for 18S rRNA. The cardiac mRNA levels were analyzed on POD5. Figure 4 indicates the number of copies of each of the three representative mRNAs measured in the syngeneic grafts or allografts obtained from three individuals. The relative amount of each mRNA was normalized to that of 18S rRNA. All experiments were analyzed in three mice per time point and expressed as the mean AE SEM. The significance level was set as P < 0.05 compared with syngeneic grafts on day 5.

The SFC method
We observed that the distribution of the mean value and variance of one probe signal is non-linear (Fig. S1). The information from neighboring probes can usually be borrowed to improve the statistical power [2]. SFC was introduced to estimate variance for each probe, rather than obtaining this from all samples; it takes information from the neighbors of that probe with an adjustable bin size b. As we set up the default value of b as 1000, the variance of cases and controls in one probe can be obtained by calculating the median for those probes separately. Eventually, by following Eqn (1), we can obtain the statistic SFC for every probe, and the P-value can be further estimated based from these.

SFC had a better sensitivity and specificity based on simulation data
We investigated the FPR and the FNR of the three methods under the null hypothesis and alternative hypothesis. As indicated in Eqn (2), signals of the null hypothesis were generated by a simple formula, y = x, with the Gaussian noise added. The basic formulas are adjustable with the parameters k. The signals of the alternative hypothesis were described by Eqn (3), with variable values of h and the portion of real positive calls. We calculated the FPR and FNR for every different h and portion of real positive calls with a 0.05 significance threshold and 100-times simulation ( Table 1).
Under the null hypothesis, the rates of the three methods are all near the significance threshold between 5% and 6% (Fig. 1A). Under the alternative hypothesis, SFC had a better performance for FPR than the other two methods generally (Fig. 1B). With an increasing h and portion of real positive calls, the FPR of SFC showed a decreasing bias, whereas Limma and the t test showed a positive bias with these parameters (Table 1). For the FNR, as the h and portion of real positive calls increased, Limma showed a faster decline than the t test, while SFC had a lower FNR than Limma and performed better with larger h and portion of real positive calls. Interestingly, SFC shows a relatively small number of calls (from 4.9% to 10.5%, Table 1), while Limma and the t test calls a larger set in this situation. In sum, comparing with Limma and the t test at the significance threshold of 0.05, SFC had a better sensitivity and specificity, especially with a larger value of differential expression fold change (h = 50%).

Reproducibility of SFC is better than Limma and the t test based on MAQC and mouse transplantation data
Reproducibility is an indispensable estimator for the experiments and algorithms [12,13]. We chose both the MAQC dataset and the mouse cardiac transplantation data to assess the reproducibility of SFC, Limma and the t test.
We calculated the reproducibility of the top 100 and top 1000 genes for MAQC by using the three methods. For the interplatform comparison, the heatmap shows that SFC performed a better reproducibility than Limma and the t test among six platforms when detecting both the top 100 and the top 1000, while for intra-platform reproducibility, all three methods did not perform well in detecting either the top 100 or the top 1000 significant genes ( Fig. 2A,B). The same operations were conducted in the mouse cardiac transplantation data, where SFC also showed a better performance than the others (Fig. 2C). Therefore, according to better performances of reproducibility for both the MAQC data and the mouse transplantation data, SFC is more stable than Limma or the t test.
Intersected significances from mouse transplantation data were found by SFC and validated by qRT-PCR We further utilized the three methods to analyze the mouse organ transplantation data and validated the results. After the experimental process generating CEL files from mouse tissues, we conducted these methods on the expression data of POD5 of cardiac transplantation and POD5 and POD8 of hepatic transplantation. According to SFC, 178 significant genes were differentially expressed in the cardiac allografts compared with isografts, including 158 overexpressed genes and 20 underexpressed genes (Fig. 3). There were also 362 genes (263 overexpression and 99 underexpression) having significantly different expression in the hepatic POD5 allografts compared with isografts, and 389 genes (258 overexpression and 131 underexpression) having significantly different expression in the hepatic POD8 allografts compared with isografts. Based on these, an intersection of these three groups was obtained that included 52 important genes, in which they are all overexpressed for cardiac transplantation and 51 overexpressed and one underexpressed for hepatic transplantation (Fig. 3). At the same time, the calling sets of significant genes underlying Limma and the t test (Fig. S4A,B) showed no intersecting ones.
We further performed qRT-PCR for the calls derived from SFC to validate the fold change of the mRNA expression. Nineteen mRNAs, which were upregulated in both the cardiac and the hepatic allografts compared with isografts, were randomly selected (Tables 2 and S3). Being consistent with the results of microarray, a significantly higher amount of mRNA expression was detected in allografts versus isografts in cardiac (Fig. 4A) and hepatic (Fig. 4B) allografts.

Discussion
Microarray is widely used and accepted as a stable, well established and less costly technology to investigate gene expression data [1,8,14,15]. In this study based on microarray data, we established a novel method, SFC, to detect differential expression and compared it with the t test and Limma. According to Eqn (1), the parameter b can be adjusted to control the nearby number of probes, which contribute the variance of the central probe. We set 1000 as default, and users are able to customize this value based on a different number size of microarray probe. For the simulation data, the parameter configurations (h and k) of the null hypothesis and alternative hypothesis also can be adjusted (Eqns 2 and 3) [7]. Moreover, we calculated the FPR and FNR based on different significance levels (P = 0.01 and 0.001) for different values of h and k. With a more stringent significance level (0.05, 0.001), the FPRs were decreased while the FNRs

0%
50% 100% were increased, which was observed by all three methods (Figs 1, S2 and S3, Tables 1, S1 and S2). Notably, when P = 0.001, the t test gave a high FPR (48%, h = 50%, and the true positive gene percentage was 10%) and Limma performed with a high FNR (sometimes more than 90%). This suggests the t test will give more positive hits with a high FPR, while Limma will report fewer hits to reduce the FPR but miss some true positive ones. Importantly, SFC can give a good balance of FPR and FNR, and perform well for both FPR and FNR with a stringent significance level. Statistical correction (e.g. the Bonferroni correction) is often introduced for multiple comparisons to adjust the P-value and control the false discovery rate [16]. We also analyzed the mouse transplantation data by the other two methods (Limma and the t test) with different significance levels (P = 0.05, 0.001 and 0.05 with the Bonferroni correction). Limma and the t test had a large number of positive hits when the P-value was 0.05 in three phases (Figs S5 and S6). When the level of significance was P < 0.001, the positive hits by Limma and the t test decreased a lot while by SFC the number stayed relatively stable. When the P-value was stringent at 0.05 with the Bonferroni correction (Figs 3 and S4), SFC still reported 52 significances overlapping with three phases, but Limma and the t test showed no overlapping significance. The results of the t test showed no shared significance with SFC. Intriguingly, in 67 significances for cardiac POD5 reported by Limma (Fig. S4), 30 genes showed in the cardiac POD5 result for SFC, and 16 showed in the 52 significances. Besides, for hepatic POD5 and POD8 by Limma, 4 out of 7 (POD5) and 19 out of 36 (POD8) significant genes were observed in the corresponding results of SFC, and 2 out of 5 (overlapping in POD5 and POD8) significant genes appear in the 52 genes from SFC. As 19 of 52 genes from SFC were randomly selected and all passed the validation of qRT- Fig. 3. Venn diagram of significant genes analyzed by SFC with the level of significance set at P < 0.05 after the Bonferroni correction. The overall numbers of significant genes in three phases are shown outside, which are followed by numbers in parentheses showing the counts of overexpressed genes versus underexpressed ones. The circle at the top represents POD5 for heart; the circle at the bottom left represents POD5 for liver and the circle at the bottom right represents POD8 for liver. PCR, these results indicated that SFC gave a more stable result than the t test and Limma. We therefore investigated the functions of these 52 genes (Table S4), revealing the most significant pathways were graft-versus-host disease (mmu05332) and allograft rejection (mmu05330). Moreover, immune system response (e.g. mmu04612, mmu04660, GO: 0006955) and positive regulation (e.g. GO: 0050863, GO: 0051249, GO: 0050870) were also activated. All these enrichment analyses indicated a reaction of transplantation rejection in vivo and functional changes of the immune system both at the cardiac and at the hepatic level after 5 days of allografts [6,17,18].
In conclusion, based on the quality control experimental data and simulated data, SFC performed better than Limma and much better than the t test by using the nearby information of one probe in pooled probes. We utilized SFC for the real data of mouse transplantation models, and it reported a more stable and convincing set with 52 significant genes, revealing insights into pathway and gene expression changes after both cardiac and hepatic allografts. Nineteen genes were further randomly picked up and validated by qRT-PCR. We suggest SFC is a new and effective approach that can detect differential expression and help to obtain more reliable information in microarray studies.

Supporting information
Additional Supporting Information may be found online in the supporting information tab for this article: Fig. S1. Distribution of mean and variance of sample microarray signals in each probe derived from the MAQC data. Fig. S2. Bar graphs of FPR and FNR from the three methods under the null hypothesis (H0) and the alternative hypothesis (H1) with the level of significance set at P < 0.01. Fig. S3. Bar graphs of FPR and FNR from the three methods under the null hypothesis (H0) and the alternative hypothesis (H1) with the level of significance set at P < 0.001. Fig. S4. Venn diagrams of significant gene numbers analyzed by the t test and Limma with the level of significance set at P < 0.05 after the Bonferroni correction. Fig. S5. Venn diagrams of significant gene numbers analyzed by the t test, Limma and SFC with the level of significance set at P < 0.05. Fig. S6. Venn diagrams of significant gene numbers analyzed by the t test, Limma and SFC with the level of significance set at P < 0.001. Table S1. Evaluation of three methods with the level of significance set at P < 0.01. Table S2. Evaluation of three methods with the level of significance set at P < 0.001. Table S3. Primer sequences for qRT-PCR. Table S4. GO term and pathway enrichment analysis based on the 52 significant genes.