We searched the GEO database (https://www.ncbi.nlm.nih.gov/geo/) (Barrett et al., 2013) for human AAV-, ANCA-associated glomerulonephritis (AAGN), and ANCA-associated tubulointerstitial nephritis (AATN)- related gene profiles using the keywords “ANCA,” “AAV,” “AAGN,” and “AATN.” Data from 30 AAV patients in the GSE129752 dataset (Gill et al., 2021) based on the GPL10999 platform (Illumina Genome Analyzer IIx, Homo sapiens) were used for subsequent analysis.The GSE158163 dataset (Grayson et al., 2018) is based on the GPL16791 platform (Illumina HiSeq 2500, H. sapiens), from which we extracted data from 16 control samples (virus negative and positive control groups) for subsequent analysis. The GSE104948 dataset is based on the GPL22945 platform (Affymetrix Human Genome U133 Plus 2.0 Array) and the GPL24120 platform (Affymetrix Human Genome U133A Array), and we extracted glomerular expression data from 21 normal and 22 AAGN patient samples. The GSE104954 dataset (Grayson et al., 2018) was also based on the GPL22945 platform and the GPL24120 platform, and we extracted renal tubular expression data from 21 control and 21 AATN patients (Figure 1).

For the RNA sequencing data, the sequencing count values of the GSE129752 and GSE158163 datasets were downloaded and converted to Fragments Per Kilobase of exon model per Million mapped fragments (FPKM) values using the biomaRt package in R. Since the optimal input matrix type for WGCNA was FPKM, the FPKM values were no longer converted to TPM values. Next, the FPKM matrices of the two datasets were subjected to missing value completion and data cleaning. We removed rows (genes) with NA values greater than 50% and columns (samples) with NA values greater than 50%, and furthermore, we used the impute.knn function of the R package impute to perform missing value completion by setting the Number of neighbors to K to complete the missing data. The expression values were taken as log2 (Exp+1) and the mean value of the duplicated genes was taken as their expression. We then combined the two datasets by limma and removed the batch effect, and analyzed the data distribution before and after combining the two datasets by principal component analysis (PCA). The anosim function in the vegan package of the R language was used to calculate the P and R values. For the microarray data, we used Perl scripts to preprocess the raw matrices of GSE104948 and GSE104954. Probe identification numbers were converted to gene symbols and empty probes were removed based on the annotation information contained in each platform file. When multiple probes were matched to the same gene, the average expression value was used to determine the expression level of the gene. If there is a large variation in the expression of a gene, log (Exp+1) is used to measure the expression level.

The “WGCNA” R package (Langfelder and Horvath, 2008) was used to identify the key modules most closely related to AAV in the matrix of merged GSE129752 and GSE158163 expression profiles. We extracted the AAV group’s pediatric vasculitis activity score (PVAS) scores from the original dataset as input data for WGCNA. A sample clustering dendrogram was constructed using the cluster function, and the topological overlap matrix (TOM) was constructed using the pickSoftThreshold function by WGCNA to determine the optimal soft threshold. We used candidate power values (1–20) to determine the average connectivity and independence of various modules. Dynamic shear trees were used to identify gene modules. Next, we measured the association between modules and sample traits using gene significance (GS) values and module membership (MM) values, and key modules were identified. We set |GS| to >0.3 and |MM| to >0.7 to filter hub genes in accordance with the official WGCNA guidelines and prior application examples so as to obtain the genes most relevant to the traits in the key module (Langfelder and Horvath, 2008; Tang et al., 2018). We performed GSEA (Wu et al., 2021b) using the KEGG and REACTOME gene sets in the GSEA dataset [c2.cp.v7.2.symbols.gmt (curated)]. A false discovery rate of <0.25 and p.adjust <0.05 (Subramanian et al., 2005) were set as the screening criteria for enrichment pathways.

The module hub genes of WGCNA were intersected with the AAV-related key genes and nephropathy-related key genes in the Genecards database (Supplementary Table S1) to identify co-key genes of AAV. The “bioconductor” package “org.Hs.e.g.,.db” (Carlson, 2022) and “clusterProfiler” package were installed in R (R Foundation for Statistical Computing, Vienna, Austria), which was used for GO (Ashburner et al., 2000) analysis and KEGG (Kanehisa and Goto, 2000) pathway enrichment analysis while adopting p < 0.05 and p.adjust <0.05 as criteria to screen the co-key genes.

The expression levels of co-key genes in the glomerular dataset GSE104948 and the tubular dataset GSE104954 were identified. We extracted the gene expression information after preprocessing the data and performed the Wilcoxon rank sum test. Then, we used the “PRROC” (Gatto et al., 2014) package in R to examine the diagnostic efficacy of co-key genes by performing receiver operating characteristic (ROC) curve analysis.

NMF (Gaujoux and Seoighe, 2010) clustering was performed on the eight co-key genes in the GSE104948 and GSE104954 datasets to identify potential features of interest in the gene-expression profiles. The distinct molecular subtypes of AAV-related renal injury were identified based on a cophenetic and consensus matrix. These distinct molecular subtypes were systematically assessed by ssGSEA analysis using the “GSVA” (Hänzelmann et al., 2013) R package to explore the immune infiltration and mechanisms of different subtypes of AAV-related renal injury.

The Nephroseq database (http://v5.nephroseq.org/) (Zheng et al., 2017) is an online, free public platform. A correlation analysis of co-key genes and clinical features using the Nephroseq database was performed to evaluate the clinical prognostic value of key genes for AAV-related renal injury.

The WGCNA package was performed to identify the key modules most closely related to AAV in the matrix of merged GSE129752 and GSE158163 expression profiles (which contained a total of 21,113 genes). Differences in the integrated dataset after removal of batch effects were mainly driven by intrinsic heterogeneity between the Disease and Control groups (Supplementary Figure S1). The control and AAV sets in WGCNA were automatically clustered into different branches (Figure 2A). The merged shear height of the dynamic shear tree was 0.85 for module identification. The minimum number of genes in each network module was set to 50, resulting in a total of 21 gene modules (Figure 2B). The most strongly correlated positive and negative modules were chosen as critical modules for discerning the pathogenesis of AAV traits. Our results indicate that the brown module (containing a total of 7,290 genes) was significantly negatively correlated with the AAV trait (cor = −0.29, p = 0.05) (Figure 2C). By setting |GS| > 0.3 and |MM| > 0.7 as criteria to screen the critical module genes for the pathogenesis of the AAV traits, we identified 462 genes associated with AAV pathogenesis (Figure 2D).

We performed a comprehensive GSEA analysis of the screened brown module genes to discover the function and pathways of the key modules and explore the pathogenesis of AAV. The results showed (Figure 3A) that module genes were significantly enriched in the following immune- and inflammation-related pathways: antigen processing and presentation, autoimmune thyroid disease, immunomodulatory interactions, natural killer cell–mediated cytotoxicity, MAPK signaling pathways, and Leishmania protozoa infection. It has been suggested that immune, infection, and cellular responses are critical for the progression of AAV-related renal injury, which was consistent with the existing knowledge about AAV, and the accuracy of the module genes was also verified.

The WGCNA module hub genes were then intersected with the AAV-related key genes and nephropathy-related key genes in the Genecards database to identify key genes of AAV-related renal injury. We identified a total of 11 co-key genes, including PGR, AGTR2, BDKRB1, FGA, 1L1RAPL2, ANGPTL2, PCDH11Y, CSF2, ADAMTSL1, MIRLET7E, and SAA4 (Figure 3B); then, we validated the expression levels of these co-key genes in the ANCA glomerular dataset GSE104948 and the ANCA tubular dataset GSE104954 to verify that they have similar expression profiles in independent datasets. Since the expression profiles of MIRLET7E, ADAMTSL1 and SAA4 were missing in the independent dataset, we only validated the expression profiles of the remaining eight genes. Our results indicated that the expression levels of PGR, AGTR2, BDKRB1, FGA, ANGPTL2, and PCDH11Y in the GSE104948 expression profile data were not significantly different between the control and AAGN groups (Figure 3B). However, the expression levels of 1L1RAPL2 (p < 0.001) and CSF2 (p < 0.05) were significantly lower in the AAGN group than in the control group (Figure 3C).

In the tubular dataset GSE104954, the expression levels of PGR, BDKRB1, FGA, 1L1RAPL2, ANGPTL2, and CSF2 were significantly downregulated in patients with AATN (Figure 3E). We examined the diagnostic sensitivity and specificity of the eight hub genes by performing ROC curve analysis. In the glomerular dataset GSE104948, the area under the ROC curve (AUC) values for 1L1RAPL2 and CSF2 were both >0.7, suggesting that these genes have good diagnostic efficacy as AAGN markers, with the larger AUC value being 0.816 (1L1RAPL2) (Figure 3D). In the tubular dataset GSE104954, as shown in Figure 3F, the AUC values of all six genes were >0.7 (PGR: 0.850, BDKRB1: 0.706, FGA: 0.796, 1L1RAPL2: 0.808, ANGPTL2: 0.780, and CSF2: 0.730), suggesting that these six genes could be AATN biomarkers with high specificity. These genes have good diagnostic sensitivity and specificity for AAV-related renal injury.

The biological functions and key signaling pathways of the eight candidate key genes were analyzed using GO and KEGG enrichment analyses. Our results showed that the co-key genes were mainly enriched for GO entries, including G protein–coupled peptide receptor activity, transcription factor activity, chemoattractant activity, receptor antagonist activity, endoplasmic reticulum lumen, platelet alpha granule lumen, extrinsic apoptotic signaling pathway, hormone secretion, vascular processes in the circulatory system, heterotypic cell–cell adhesion regulation, and other processes (Figure 4A). The KEGG pathways involved those designated for complement and coagulation cascades, the renin–angiotensin system, neuroactive ligand–receptor interaction, acute myeloid leukemia, the Fc epsilon RI signaling pathway, and others (Figure 4B).

To further explore whether the identified co-key genes are meaningful across renal injury subtypes, we performed NMF clustering with the eight co-key genes for molecular subtype identification using the glomerular dataset GSE104948 and the tubular dataset GSE104954. In the glomerular dataset GSE104948, a sharp decrease in the cophenetic correlation coefficient at rank r = 2 indicated that the co-key genes could classify glomerular disease in AAV-related kidney injury into a pair of subgroups (Figures 5A, C). Meanwhile, in the tubular dataset GSE104954, the highest cophenetic correlation coefficient at rank r = 3 suggested that the co-key genes could classify renal tubular disease in AAV-related kidney injury into three subgroups (Figures 5B, D). Importantly, these findings may offer a novel approach to the classification of AAV-related kidney injury. Next, we further investigated the immune infiltration among these five subtypes. The results of ssGSEA showed that, between both subtypes of AAGN, subtype 1 had a higher level of infiltration by effector memory CD4 T-cells and memory B-cells than subtype 2, suggesting that each cell type has a crucial role in the progression of subtype 1 AAGN renal injury (Figure 5E). Separately, the distribution of the immune infiltration of 28 immune cells in the three subtypes of AATN is shown in Supplementary Figure S2; notably, CD56dim natural killer infiltration in subtype 1 of AATN was significant (p = 0.048), and neutrophil cell infiltration in subtype 2 was significant (p = 0.0097), while memory B-cell infiltration was significant in subtype 3 (p = 0.0024). This indicated that these three cell types have a significant effect on the progression of the three subtypes of AATN.

In summary, we determined that AGTR2, ANGPTL2, BDKRB1, CSF2, FGA, IL1RAPL2, PCDH11Y, and PGR are key genes in patients with AAV-related kidney injury, with diagnostic sensitivity and specificity. We further validated the relationship between co-key genes and clinical traits in the Nephroseq database to assess the association between co-key genes and prognostic factors of AAV-related kidney injury. Low expression of co-key genes was associated with a low glomerular filtration rate (BDKRB1: r = 0.338, p < 0.001; PCDH11Y: r = 0.287, p < 0.001; IL1RAPL2: r = 0.243, p < 0.001) (Figures 6A–H), suggesting that the expression of these genes may participate in the process of AAV-related kidney injury.But it is not clear whether its correlation can predict its poor prognosis, which requires further experiments and large clinical cohorts for validation.