Microarray datasets from kidney transplant patients with AMR were identified by searching the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo). The search was conducted in September 2017. The key word used in the search was ‘kidney transplantation’ and the Homo sapiens search filter was applied. Each dataset was manually curated to select biopsy or peripheral blood samples from Homo sapiens. Only original experimental articles that compared the expression levels of mRNAs between patients with AMR and those without AMR (controls) were retained. All eligible publications met the following inclusion criteria: i) The studies were associated with the diagnostic value of mRNA for a diagnosis of AMR; and ii) the studies provided sufficient data with which to assess the diagnostic value of mRNA in AMR. The exclusion criteria included: i) Duplicate publications; ii) studies without sufficient data; and iii) letters, reviews, editorials, meeting abstracts and case reports. Two researchers independently screened the list of publications and evaluated the possibility of inclusion. Inter-rater agreement was assessed with Cohen's κ statistic.

The following information was extracted from each identified study: The GEO accession number, platform, sample type, number of cases and controls, references and expression data. Two independent reviewers extracted data from the original studies and a consensus was reached, or a third reviewer would resolve any discrepancies between the two reviewers. The studies used the gene expression platforms HG-U133 Plus 2.0, HuGene-1_0-st (both Affymetrix; Thermo Fisher Scientific, Inc., Waltham, MA, USA) and Agilent-014850 (Agilent Technologies, Inc., Santa Clara, CA, USA). In total, 9 studies were initially identified. The remaining 6 studies (GSE36059, GSE44131, GSE50084, GSE51675, GSE64261 and GSE93658) were included in the final analysis. The detailed information about the downloaded datasets is summarized in Table I. GSE36059, GSE44131, GSE50084 and GSE51675 were used for the discovery set, while GSE64261 and GSE93658 were used for the replication set. After analyzing each set, the combined dataset was used for the analysis.

Each individual dataset was preprocessed using the log2 transformation and normalization approach. To combine the results of the individual studies and to obtain a list of more robust DEGs between the control and AMR cases, the guidelines provided by Ramasamy et al (2008) (19) for a meta-analysis of gene expression microarray datasets were followed. The R packages MetaQC (20) and MetaDE (17,21–23) were used for quality control (QC) and for the identification of DEGs, respectively. MetaQC implements the 6 quantitative QC measures of internal QC (IQC), external QC (EQC), accuracy QC of the featured genes (AQCg), accuracy QC of the pathway (AQCp), consistency QC in the ranking of featured genes (CQCg) and consistency QC in the ranking of the pathway (CQCp). In addition, the mean rank of all QC measures in each dataset was computed as a quantitative summary score by calculating the ranks of each QC measure among all included datasets. All probe sets on the three different platforms were re-annotated to the most recent National Center for Biotechnology Information Entrez Gene Identifiers (Gene IDs), and the Gene IDs were used to cross-map genes among the three different platforms. When multiple probes matched the same gene symbol, probes presenting the greatest inter-quartile range were selected. Only genes present in all of the selected platforms were considered. The moderated t-statistic was used to calculate the P-values in each dataset, and a meta-analysis was conducted with the MetaDE package to identify DEGs using the Fisher's (24), maximum P-value (maxP) (25), rth ordered P-value (roP) (26), Stouffer (27) and naive rank summation (SR) (28) methods.

In order to select genes and pathways associated with AMR, genes were annotated using GSA, as conducted by GSA-SNP software version 1.0 (29). The GSA-SNP analysis uses Gene Ontology (GO) (30), the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (31,32) and the Molecular Signatures Database (MSigDB) (33). Genes that showed nominal significance levels of P<0.05 were selected.

To identify the phenotype-related modules and protein-protein interactions (PPIs), the PPI databases BioGRID (34) and Agile Protein Interactomes Data Server (35) were downloaded, and the list of the identified genes from the present study was imported into Cytoscape version 3.5.1 (36). The molecular complex detection (MCODE) clustering algorithm was used to identify sub-network modules (37). A network score was calculated based on the complexity and density of each sub-graph. A module with an MCODE score >2 was considered significant. Post-filtering was performed to remove low-quality modules. During the filtering process, the parts of each module that showed consistent expression and high connectivity levels were selected to constitute the final module through a manual review.

The study selection process is presented in Fig. 1. Following the search and selection process, the 6 studies of GEO36059, GEO44131, GSE50084, GSE51675, GSE64261 and GSE93658 from whole blood or peripheral blood mononuclear cells and biopsy samples of renal transplant patients met the inclusion criteria. The value of Cohen's κ coefficient was 0.81, indicating a sufficient level of agreement (coefficient value >0.80). All probe sets were re-annotated with the most recent Gene IDs and then mapped, yielding 27,047 common genes across three different platforms. The resulting dataset contained 328 samples from the control and 132 samples from the AMR cases for the discovery set, and 21 samples from the control and 38 samples from the AMR cases for the replication set. Principal component analysis plots for the discovery set were plotted with MetaQC to visualize the quality of the studies via a systematic analysis (Fig. 2). The first two PCs captured 99% of the variance.

Five systematic analysis methods were performed by combining the P-values in the MetaDE package. These were the Fisher's, maxP, roP, Stouffer and SR methods. A total of 608 DEGs were selected by the maxP and roP methods (P<0.05) in the discovery set (Fig. 3); 291 genes were upregulated in AMR cases compared with those in the control cases, while 317 genes were downregu-lated in AMR cases compared with those in the control cases. C-X-C motif chemokine ligand 10 (CXCL10), CXCL9 and guanylate binding protein 1 (GBP1) were most significantly upregulated among the AMR cases in the discovery set. A total of 317 DEGs were selected by the maxP and roP methods (P<0.05) in the replication set; 219 genes were upregulated in AMR cases compared with those in the control cases, while 98 genes were downregulated in AMR cases compared with those in the control cases. CXCL10, CXCL9 and GBP1 were most significantly associated with AMR in the replication set. The 20 top-ranked differentially downregulated genes associated with AMR in the replication set are listed in Table II. In total, 608 DEGs were identified by the maxP and roP methods (P<0.05) in the combined discovery and replication sets; 248 genes were upregulated in AMR cases compared with those in the control cases, while 360 genes were downregulated in AMR cases compared with those in the control cases. Among the genes that showed a significant association with AMR in the present study, CXCL10, CXCL9, GBP1, complement C1q A chain (C1QA), complement C2 and interferon-induced protein with tetratricopeptide repeats 2 (IFIT2) were most significantly upregulated, while renin binding protein (RENBP), kynurenine 3-monooxygenase (KMO), thyroid hormone receptor interactor 6 (TRIP6) and glutathione S-transferase θ1 (GSTT1) were most significantly downregulated (Table III).

In the combined discovery and replication sets, the GO algorithm identified >231 gene sets associated with AMR (all P<0.05). The highly ranked gene sets are shown in Table IV. The KEGG pathway identified 31 gene sets associated with AMR (all P<0.05). The MSigDB identified >1,720 gene sets associated with AMR (all P<0.05). The regulation of the immune response gene set, in this case β2-microglobulin (B2M), baculoviral IAP repeat-containing 3 (BIRC3), C1QA, CD3e molecule (CD3E), CD79A, major histocompatibility complex class II DM α (HLA-DMA), HLA-DOA, HLA-DOB, HLA-DPB1, HLA-DQA1 and HLA-DRA, was significantly associated with AMR (P<1×10⁻⁵). Allograft rejection genes, in this case HLA-DMA, HLA-DOA, HLA-DOB, HLA-DPB1, HLA-DQA1, HLA-DRA, HLA-DMB, tumor necrosis factor (TNF) and granzyme B, were also significantly associated with AMR (P<1×10⁻⁵). Hallmark interferon-γ response genes, specifically transporter 1, ATP binding cassette subfamily B member (TAP1), MX dynamin-like GTPase 2, CXCL10, GBP4, CXCL9, interleukin 10 receptor subunit α, SAM and HD domain containing deoxynucleoside triphosphate triphosphohydrolase 1 and intercellular adhesion molecule 1, were also significantly associated with AMR (P<1×10⁻⁵).

There were two modules with a MCODE score of >2 in AMR cases. Fig. 4 shows the sub-network modules. One module consisted of immunoglobulin heavy constant γ1 (G1m marker) (IGHG1), Fc fragment of IgG receptor IIIa (FCGR3A), FCGR3B, amyloid P component serum (APCS) and glycerol kinase (GK), and the other module consisted of CD74, HLA-DMA, HLA-DMB, HLA-DRA, HLA-DRB1, GBAA0374, proteasome subunit β5 (PSMB5), PSMB8, PSMB9, TAP1 and TAP2.