Explanted lung tissue samples for RNA sequencing analysis were obtained from the Pusan National University Yangsan Hospital (PNUYH) biobank. These included lung tissues from three patients with IPF and lung tissues from three patients with AEs of IPF. The criteria for the diagnosis of IPF were those of the American Thoracic Society and European Respiratory Society [14]. An AE of IPF was defined as an acute, clinically significant respiratory deterioration characterized by evidence of new widespread alveolar abnormality according to the 2016 revised criteria [4]. All cases were reviewed by expert pulmonologists and pathologists. Detailed clinical information regarding the subjects is provided in Table 1. The mean forced vital capacity (expressed as a percentage of the normal expected value; FVC %) and diffusing capacity for carbon monoxide (expressed as a percentage of the normal expected value; DLCO %) of patients are provided in Table 1. Lung tissue acquisition and analysis were conducted with the appropriate approval from the Institutional Review Board of Pusan National University Yangsan Hospital (PNUYH IRB No 05-2020-019).

Total RNA from cell samples was extracted using a TRIzol reagent kit (Invitrogen), according to the manufacturer’s protocol, and RNA integrity was assessed using a TapStation RNA screentape. After total RNA extraction, RNA libraries were independently prepared using the Illumina TruSeq Stranded Total RNA Library Prep Gold Kit (Illumina, Inc., San Diego, CA, USA). The cleaved RNA fragments were copied into first-strand cDNA using SuperScript II reverse transcriptase (Invitrogen, Carlsbad, CA, USA) using random primers. Qualified libraries were sequenced on an Illumina NovaSeq platform (Illumina, Inc., San Diego, CA, USA). Total RNA concentration was calculated using the Quant-IT RiboGreen Assay Kit (Invitrogen, Waltham, MA, USA). To determine the values of DV200 (percentage of RNA fragments $>$200 bp), samples were run on the TapeStation RNA screentape (Agilent, Wilmington, DE, USA). A total of 100 ng of total RNA was subjected to sequencing library construction using the Agilent SureSelect RNA Direct kit (Agilent, Wilmington, DE, USA) according to the manufacturer’s protocol. Briefly, the total RNA was first fragmented into small pieces using divalent cations at elevated temperatures. The cleaved RNA fragments were copied into first-strand cDNA using random primers. This was followed by second-strand cDNA synthesis. These cDNA fragments then underwent end repair, the addition of a single “A” base, and the ligation of the adapters. The products were purified and enriched by PCR to create a cDNA library. To capture the human exonic region, the Agilent SureSelect XT Human All Exon v6+UTRs Kit (Agilent, Wilmington, DE, USA) was used according to the standard Agilent SureSelect Target Enrichment protocol. A cDNA library (25 ng) was mixed with hybridization buffers, blocking mixes, RNase block, and 5 $\mu$L of SureSelect XT Human All Exon v6+UTRs capture library (Agilent, Wilmington, DE, USA). Hybridization to the capture baits was conducted at 65 ${}^{\circ }$C using the heated thermal cycler lid option at 105 ${}^{\circ }$C for 24 h in a PCR machine. The captured library was washed and subjected to a second round of PCR amplification. The final purified product was then quantified using qPCR according to the qPCR Quantification Protocol Guide (KAPA Library Quantification kits for Illumina Sequencing platforms) and qualified using the TapeStation DNA screentape D1000 (Agilent, Wilmington, DE, USA). The indexed libraries were analyzed using the NovaSeq system (Daejeon, Korea). Paired-end (2 $\times$ 100 bp) sequencing was performed by Macrogen Inc. The gene expression data obtained through RNA-sequencing is provided in the Supplementary file 1.

Expression profiling array data for IPF was obtained from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/; accession numbers GSE44723 and GSE10667). All datasets included data related to gene expression in the lung tissue of the AE group and/or that of the stable IPF group. The data were downloaded for DEG analysis. The characteristics of the included datasets are listed in Table 2 (Ref. [15, 16]).

DEGs of each dataset were identified using the GEO2R online analytic tool provided by GEO. Adjusted p 0.05 and $|$log2FC$|$$>$ 0 were the cutoff criteria. For the adjusted p-value, false discovery rate methods were used. The probe identifier was converted to a genetic symbol based on the annotation file downloaded from the GEO. When multiple probes corresponded with a single gene, the average level of the probes was calculated as the expression value of the gene. Probes that did not correspond to the genetic symbols were removed. PNUYH RNA-seq data were quality checked using the fastqc (v 0.11.9) tool [17]. This check was performed after removing the sequence with a thread score 30 and Illumina adaptor sequence. The reference genome used for alignment was GRCH38 using the hisat2 tool (v 2.2.1) [18]. Stringtie (v 2.1.5) was used to obtain the read count for the DEG analysis [19]. DEGs were confirmed using the limma R package. Data preprocessing was performed using a voom provided in the limma package [20]. To screen the intersectional genes that were significantly expressed in each dataset, the R package “venn” was used to plot Venn diagrams [21]. Among the DEGs from the three datasets, only DEGs that commonly appeared in all three datasets were considered as significant DEGs and were selected for the subsequent analysis.

A PPI network was constructed using the Search Tool for the Retrieval of Interacting Genes (STRING) database [22] to identify the relevant pathways and functions of common DEGs. The minimum required interaction score used to construct the PPI network was 0.4. The molecular complex detection (MCODE) [23] was used for optimal module selection. The parameters used in MCODE were: degree cutoff = 2; cluster finding, node score cutoff = 0.2; k-core = 2; and maximum depth = 100. CytoHubba [24] was used to identify the highest-linkage hub genes in the network. The parameters used in CytoHubba were the top 10 nodes ranked by degree and the display option in which the expanded subnetwork was displayed.

GO analyses were performed using the R package “clusterProfiler” (v 3.14.3) [25]. The common DEGs identified in each dataset were included in the GO term enrichment analysis, which was performed separately for up- and downregulated genes. Terms with five or more associated genes and two or more DEGs in each experiment were included. The adjusted p-values were calculated using the Benjamini-Hochberg method. Statistical significance was set at p 0.05.

The flowchart for this study is shown in Fig. 1. Detailed information on the datasets is presented in Table 1. DEG analysis of GSE10667 identified 4155 DEGs, of which 1832 genes were upregulated and 2325 genes were downregulated (Fig. 2). DEG analysis of PNUH RNA-seq identified 1819 DEGs, of which 897 genes were upregulated and 922 genes were downregulated (Fig. 2). DEG analysis results for each dataset are in Supplementary files2.1, 2.2, and 2.3. We performed GO analysis of common DEGs and the results are in Supplementary files 2.4 and 2.5. Overlapping genes in each dataset included 76 upregulated genes and 135 downregulated genes, and the included genes list is in Supplementary files 2.6.

Fig. 1.

Flow chart of the study.

Fig. 2.

Analysis of the differentially expressed genes. (A) Venn diagram of upregulated genes. (B) Venn diagram of downregulated genes. (C) Volcano plot of each cohort. In the plot, genes that met the p-value and log${}_2$ fold change cutoff value are denoted in red; genes that passed only the p-value are denoted in blue; a gene that passed only the log${}_2$ fold change is denoted in green; and genes that passed both the p-value and log${}_2$ fold change are denoted in black. (D) Heatmap of common DEGs for each cohort. The common upregulated genes are denoted in red and the common downregulated genes are denoted in green.

To identify the PPIs of common DEGs, the STRING network-based protein interaction assay was used to create a PPI network. To confirm the interaction between common DEGs, a PPI network was constructed. Three significant modules (confidence score $>$7) containing 24 out of 211 genes were identified (Fig. 3). The largest module comprising the upregulated genes consisted of CCNA2, PTTG1, SKA3, SPAG5, TROAP, KIF14, ASPM, CENPE, and E2F7 (Fig. 3C). The second module comprising the upregulated genes consisted of histone-related genes (HIST1H2AG, HIST3H2A, HIST1H2BJ, HIST1H2BI, HIST1H2BE, and HIST1H3B; Fig. 3D). The third module comprising the downregulated genes contained six genes (EGR1, NR4A1, FOSB, FOS, CYR61, and BTG2; Fig. 3E).

Fig. 3.

Protein-protein interaction (PPI) network construction and significant module selection involving the common differentially expressed genes (DEGs). (A) PPI network of upregulated genes. (B) PPI network of downregulated genes. (C) The module with the highest confidence score among the PPI network of upregulated genes. (D) Module with the second-highest confidence score among the PPI network of upregulated genes. (E) The module with the highest confidence score among the PPI network of downregulated genes.

To validate the genes included in significant modules, the independent cohort GSE44723 was used. Eleven upregulated genes (ASPM, CCNA2, CDC25B, CENPE, CLSPN, KIF14, PTTG1, RECQL4,SKA3,SPAG5, and TROAP) and six downregulated genes (BTG2, CYR61, EGR1, FOS, FOSB, and NR4A1) were significantly altered in the cases presenting AEs of IPF compared to those presenting stable IPF (Fig. 4). Interestingly, the histone-associated genes and E2F7 gene were not significantly increased in AEs of IPF patients.

Fig. 4.

Boxplots of 17 genes with statistically significant differences in expression values in GSE44723 among 24 genes included in significant modules. (A) Upregulated genes. (B) Downregulated genes. The blue boxplot indicates stable IPF and the red boxplot indicates acute exacerbation of IPF.

To investigate the functional role and the associated signaling pathways of the common DEGs, GO analysis was performed. The total significantly enriched GO terms from the up- and downregulated common DEGs are shown in Fig. 5. The top biological processes (BPs) of the eleven upregulated and six downregulated genes associated with the AEs of IPF were nuclear division, organelle fission, skeletal muscle cell differentiation, skeletal muscle tissue development, skeletal muscle organ development, response to mechanical stimulus, regulation of neuron death, and neuron death, respectively. The top cellular components of the upregulated and downregulated genes were spindle and transcription regulator complexes. The top molecular functions were microtubule binding, tubulin binding, and RNA polymerase II-specific DNA binding transcription factor binding for upregulated proteins, and DNA binding transcription factor binding for downregulated proteins.

Fig. 5.

Gene ontology assay of 17 genes identified through protein–protein interaction network and module selection. (A) Upregulated genes. (B) Downregulated genes. The gene ontology analysis result is composed of biological process (BP), cellular component (CC), and molecular function (MF) terms. The top 10 of each term is shown. The size of the circle indicates the number of genes corresponding to each term. The closer to red, the lower the p-value.