# Integrated Analysis of a Gene Correlation Network Identifies Critical Regulation of Fibrosis by lncRNAs and TFs in Idiopathic Pulmonary Fibrosis

Authors: Fan Wang; Pei Li; Feng-sen Li
Journal: BioMed Research International
 (2020)

---

## Abstract

Idiopathic pulmonary fibrosis (IPF), the most frequent form of irreversible interstitial pneumonia with unknown etiology, is characterized by massive remodeling of lung architecture and followed by progressive loss of lung function. However, the key regulatory genes and the specific signaling pathways involved in the onset and progression of IPF still remain unclear. The present study is aimed at investigating the key role of long noncoding RNAs (lncRNAs) and transcription factors (TFs) involved in the pathogenesis of IPF through the integrated analysis of three gene expression profiles from the GEO dataset (GSE2052, GSE44723, and GSE24206). A total of 8483 differentially expressed genes (DEGs) including 988 upregulated and 7495 downregulated genes were filtered. Subsequently, following the intersection of these DEGs, 29 overlapping genes were identified and further analyzed using a bioinformatics approach. Furthermore, the protein-protein interaction (PPI) network was used to obtain 18 modules of related genes. The hub genes were identified through hypergeometric testing, which were closely associated with ubiquitin-mediated proteolysis, the spliceosome, and the cell cycle. The significant difference was observed in the expression of these key genes, such as lncRNA MALAT1, E2F1, and YBX1, in the peripheral blood of IPF patients when compared with those normal control subjects by real-time polymerase chain reaction (RT-PCR) analysis. This study indicated thatlncRNA MALAT1, E2F1, and YBX1 may be key regulators for the pathogenesis of IPF.


---

## Body

## 1. Introduction

Idiopathic pulmonary fibrosis is a chronic and progressive lung tissue damage of unknown etiology, which is characterized by the abnormal proliferation of activated fibroblasts/myofibroblasts and excessive deposition of collagen in the extracellular matrix (ECM) from adjacent alveoli to the lung parenchyma. IPF has a poor prognosis and high mortality rate with the postdiagnosis median survival rate of only 20% to 30% and the median survival of approximately 3 to 5 years [1, 2]. Due to the complexity and heterogeneity of IPF, its incidence and mortality rate, which has a positive relationship with advanced age, have shown a steadily increasing trend worldwide [3]. Although the pharmacotherapy of IPF has made certain progress over the past 5 years, the therapeutic efficacy is unsatisfactory because of the variable and unpredictable course of IPF and large individual differences [4].Increasing studies related to transcriptome, including both protein-coding mRNAs and noncoding RNAs (ncRNAs), have provided novel insights into the molecular mechanism of IPF pathogenesis. Among them, ncRNAs implicated in multiple fibrotic diseases have been divided into short and long ncRNAs (lncRNAs) based on its length of nucleotide sequences. Multiple studies have shown that lncRNAs (≥200 nucleotides) contribute to the pathogenesis and progression of IPF and gain more attention [5, 6]. However, varying proportions of transcripts that can be detected and the accuracy of measurements of changes in low-abundance transcripts reduce detection accuracy of lncRNAs in transcriptome-related lung fibrosis research. In addition, transcripts detected and measured are very different in different microarray platforms. These factors imply that some lncRNAs may be overlooked and false-positive or false-negative results may be generated [7, 8]. Based on publicly available microarray expression datasets in the Gene Expression Omnibus (GEO) database, an in-depth bioinformatics analysis of lncRNAs may provide a comprehensive understanding of not only transcriptional regulation but also posttranscriptional regulation. Hence, using bioinformatics methods to analyze the comprehensive gene network, this study was performed to identify the biological processes and pathways of differentially expressed genes (DEGs) that are involved in the pathogenic mechanism of IPF. These results may be useful in elucidating the critical regulatory mechanism of IPF from a systematic perspective and providing the relevant effective interventions to attenuate or reverse the process of lung fibrosis.


## 2. Materials and Methods

### 2.1. Microarray Data Information

NCBI-GEO (http://www.ncbi.nlm.nih.gov/geo/) is a free database repository comprising microarray/gene profile, next-generation sequencing, hybridization array, and chip data. All data were derived from GEO datasets GSE2052, GSE44723, and GSE24206. The microarray data of GSE2052 were based on GPL1739 Platforms (Amersham Biosciences CodeLink Uniset Human I Bioarray, University of Pittsburgh, PA, USA) and included 15 IPF and 11 control lung tissues (submission date: 09 December 2004) [9, 10]. The GSE44723 data were based on GPL570 Platforms (Affymetrix Human Genome U133 Plus 2.0 Array, Affymetrix, Santa Clara, CA, USA) and included 10 pulmonary fibrosis and 4 normal lung tissues (submission date: 10 April 2013) [11]. The GSE24206 data were based on GPL570 Platforms (Affymetrix Human Genome U133 Plus 2.0 Array, Affymetrix, Santa Clara, CA, USA) and included 17 IPF and 6 normal lung tissues (submission date: 01 November 2011) [12]. The total RNA of the samples was extracted to analyze the genomic profile of the RNA. All data came from expression profiling with microarrays conducted for Homo sapiens.


### 2.2. Identification of Differential Gene Expression in IPF

The original data from these datasets, including SOFT-formatted family files and Series Matrix Files, were downloaded for analysis. DEGs were identified with the R package limma (http://bioconductor.org/packages/release/bioc/html/limma.html). Unsupervised hierarchical clustering was performed to center the normalized and log2-scaled expression values on the median by using Cluster 3.0 (Fig. S1–S3). After pretreatment of the genes that came from more than one probe set, the DEGs identified with cutoff criteria of logFC>1 and P<0.05 by the classic t-test were considered statistically significant.


### 2.3. Gene Ontology and KEGG Enrichment Analysis of DEGs

Functional and pathway enrichment analyses of candidate DEGs were performed with the online bioinformatics Database for Annotation, Visualization, and Integrated Discovery (DAVID,http://david.ncifcrf.gov) (version 6.7), which can integrate biological data and comprehensively annotate the biological functional information of genes. Gene Ontology (GO) analysis can provide annotation of DEGs regarding biological processes (BPs), molecular functions (MFs), and cellular components (CCs) and allows further analysis of the bioprocesses of these genes. The Kyoto Encyclopedia of Genes and Genomes (KEGG) provides high-level functions and biological system information derived from large-scale molecular datasets generated with high-throughput experimental technologies. DAVID was applied to analyze the function of DEGs, and a P value of less than 0.05 was considered statistically significant.


### 2.4. Construction of Protein-Protein Interaction (PPI) Network and Module Analysis

Functional analysis of interactions between the candidate DEG-encoded proteins can provide a new perspective on the pathogenesis and development of IPF. The protein-protein interaction network (PPI) of DEGs was constructed with the Search Tool for the Retrieval of Interacting Genes (STRING) online database (http://string-db.org) (version 11.0) considering combined scores of interaction greater than 0.4 to indicate statistical significance, and the network was visualized in the form of modules by using ClusterONE Cytoscape plug-in (version 1.0) [13]. Cytoscape (version 3.6.1) is powerful bioinformatics software that is utilized to visualize molecular interaction networks. Then, GO and KEGG enrichment analyses of genes in the module were conducted by using DAVID.


### 2.5. Selection of the Key Genes

Using the Molecular Complex Detection (MCODE) (version 1.4.2) plug-in of Cytoscape, the hub genes were selected by means of clustering the dense connection domain based on the topology of a given network [14]. The GO and pathway enrichment analyses of hub genes were performed with the ClueGO (2.5.1) plug-in of Cytoscape [15]. Subsequently, the biological pathway relationship network of these hub genes was constructed with the Biological Networks Gene Ontology tool (BiNGO) (version 3.0.3) plug-in of Cytoscape [16]. Using the hypergeometric test of the empirical Bayes approach, key genes were obtained through calculation. A P value of less than 0.05 was considered statistically significant.


### 2.6. Subjects and Blood Samples

Considering the particular and complex nature of IPF, not all patients can undertake the invasion operation including bronchial and surgical lung biopsies to obtain the lung tissue samples. Moreover, obtaining healthy control samples would be not only extremely difficult but also restricted by ethical concerns. Due to the feasibility and convenience of obtaining blood samples, we validated the expression levels of candidate genes in the peripheral blood samples of all subjects. IPF patients (n=20) were diagnosed at the Traditional Chinese Medicine Hospital Affiliated to Xinjiang Medical University. Healthy physical examinees (n=20) were selected as the control group. The cohort of 40 subjects provided written informed consent in compliance with the code of ethics of the World Medical Association. The collection and usage of the blood samples were approved by the Medical Research Ethics Committee of Traditional Chinese Medicine Hospital Affiliated to Xinjiang Medical University (the Scientific Research Project 2018XE0109-1).


### 2.7. Validation of Key Genes by RT-PCR Analysis

Purification of RNA from blood samples from 20 IPF patients and 20 normal control subjects was performed using the TRIzol™ LS Reagent (Invitrogen, USA). The RNA was reverse-transcribed using the PrimeScript™ RT reagent kit with ɡDNA Eraser (TAKARA, Japan) according to the manufacturer’s recommendations. The cDNA from each sample was used as a template with GAPDH as an internal reference. The specific primer sequences used to amplify the 4 key candidate genes are listed in TableS1. Real-time PCR (RT-PCR) was performed using the StepOnePlus™ Real-Time PCR System (Thermo Fisher Scientific, USA). The results are represented as the means of 3 repetitions and were quantified via the 2-ΔΔct method. The mRNA levels of key genes between the IPF and normal lung tissues were compared using a paired t-test (P<0.05) using GraphPad 6.0 (GraphPad Software, La Jolla, CA, USA). The data were presented as the mean±standarddeviationX¯±S (Table S2). Counting data was assessed using a χ2 test. Multiple-group comparison was assessed using one-way analysis of variance (ANOVA) followed by the Bonferroni multiple comparison test. Comparison of two groups was assessed by a two-tailed t-test.


## 2.1. Microarray Data Information

NCBI-GEO (http://www.ncbi.nlm.nih.gov/geo/) is a free database repository comprising microarray/gene profile, next-generation sequencing, hybridization array, and chip data. All data were derived from GEO datasets GSE2052, GSE44723, and GSE24206. The microarray data of GSE2052 were based on GPL1739 Platforms (Amersham Biosciences CodeLink Uniset Human I Bioarray, University of Pittsburgh, PA, USA) and included 15 IPF and 11 control lung tissues (submission date: 09 December 2004) [9, 10]. The GSE44723 data were based on GPL570 Platforms (Affymetrix Human Genome U133 Plus 2.0 Array, Affymetrix, Santa Clara, CA, USA) and included 10 pulmonary fibrosis and 4 normal lung tissues (submission date: 10 April 2013) [11]. The GSE24206 data were based on GPL570 Platforms (Affymetrix Human Genome U133 Plus 2.0 Array, Affymetrix, Santa Clara, CA, USA) and included 17 IPF and 6 normal lung tissues (submission date: 01 November 2011) [12]. The total RNA of the samples was extracted to analyze the genomic profile of the RNA. All data came from expression profiling with microarrays conducted for Homo sapiens.


## 2.2. Identification of Differential Gene Expression in IPF

The original data from these datasets, including SOFT-formatted family files and Series Matrix Files, were downloaded for analysis. DEGs were identified with the R package limma (http://bioconductor.org/packages/release/bioc/html/limma.html). Unsupervised hierarchical clustering was performed to center the normalized and log2-scaled expression values on the median by using Cluster 3.0 (Fig. S1–S3). After pretreatment of the genes that came from more than one probe set, the DEGs identified with cutoff criteria of logFC>1 and P<0.05 by the classic t-test were considered statistically significant.


## 2.3. Gene Ontology and KEGG Enrichment Analysis of DEGs

Functional and pathway enrichment analyses of candidate DEGs were performed with the online bioinformatics Database for Annotation, Visualization, and Integrated Discovery (DAVID,http://david.ncifcrf.gov) (version 6.7), which can integrate biological data and comprehensively annotate the biological functional information of genes. Gene Ontology (GO) analysis can provide annotation of DEGs regarding biological processes (BPs), molecular functions (MFs), and cellular components (CCs) and allows further analysis of the bioprocesses of these genes. The Kyoto Encyclopedia of Genes and Genomes (KEGG) provides high-level functions and biological system information derived from large-scale molecular datasets generated with high-throughput experimental technologies. DAVID was applied to analyze the function of DEGs, and a P value of less than 0.05 was considered statistically significant.


## 2.4. Construction of Protein-Protein Interaction (PPI) Network and Module Analysis

Functional analysis of interactions between the candidate DEG-encoded proteins can provide a new perspective on the pathogenesis and development of IPF. The protein-protein interaction network (PPI) of DEGs was constructed with the Search Tool for the Retrieval of Interacting Genes (STRING) online database (http://string-db.org) (version 11.0) considering combined scores of interaction greater than 0.4 to indicate statistical significance, and the network was visualized in the form of modules by using ClusterONE Cytoscape plug-in (version 1.0) [13]. Cytoscape (version 3.6.1) is powerful bioinformatics software that is utilized to visualize molecular interaction networks. Then, GO and KEGG enrichment analyses of genes in the module were conducted by using DAVID.


## 2.5. Selection of the Key Genes

Using the Molecular Complex Detection (MCODE) (version 1.4.2) plug-in of Cytoscape, the hub genes were selected by means of clustering the dense connection domain based on the topology of a given network [14]. The GO and pathway enrichment analyses of hub genes were performed with the ClueGO (2.5.1) plug-in of Cytoscape [15]. Subsequently, the biological pathway relationship network of these hub genes was constructed with the Biological Networks Gene Ontology tool (BiNGO) (version 3.0.3) plug-in of Cytoscape [16]. Using the hypergeometric test of the empirical Bayes approach, key genes were obtained through calculation. A P value of less than 0.05 was considered statistically significant.


## 2.6. Subjects and Blood Samples

Considering the particular and complex nature of IPF, not all patients can undertake the invasion operation including bronchial and surgical lung biopsies to obtain the lung tissue samples. Moreover, obtaining healthy control samples would be not only extremely difficult but also restricted by ethical concerns. Due to the feasibility and convenience of obtaining blood samples, we validated the expression levels of candidate genes in the peripheral blood samples of all subjects. IPF patients (n=20) were diagnosed at the Traditional Chinese Medicine Hospital Affiliated to Xinjiang Medical University. Healthy physical examinees (n=20) were selected as the control group. The cohort of 40 subjects provided written informed consent in compliance with the code of ethics of the World Medical Association. The collection and usage of the blood samples were approved by the Medical Research Ethics Committee of Traditional Chinese Medicine Hospital Affiliated to Xinjiang Medical University (the Scientific Research Project 2018XE0109-1).


## 2.7. Validation of Key Genes by RT-PCR Analysis

Purification of RNA from blood samples from 20 IPF patients and 20 normal control subjects was performed using the TRIzol™ LS Reagent (Invitrogen, USA). The RNA was reverse-transcribed using the PrimeScript™ RT reagent kit with ɡDNA Eraser (TAKARA, Japan) according to the manufacturer’s recommendations. The cDNA from each sample was used as a template with GAPDH as an internal reference. The specific primer sequences used to amplify the 4 key candidate genes are listed in TableS1. Real-time PCR (RT-PCR) was performed using the StepOnePlus™ Real-Time PCR System (Thermo Fisher Scientific, USA). The results are represented as the means of 3 repetitions and were quantified via the 2-ΔΔct method. The mRNA levels of key genes between the IPF and normal lung tissues were compared using a paired t-test (P<0.05) using GraphPad 6.0 (GraphPad Software, La Jolla, CA, USA). The data were presented as the mean±standarddeviationX¯±S (Table S2). Counting data was assessed using a χ2 test. Multiple-group comparison was assessed using one-way analysis of variance (ANOVA) followed by the Bonferroni multiple comparison test. Comparison of two groups was assessed by a two-tailed t-test.


## 3. Results

### 3.1. Identification of Differentially Expressed Genes in IPF

After normalization and standardization of the raw data from these three GSE2052, GSE44723, and GSE24206 datasets (Figures1(a)–1(c)), we identified a total of 8483 aberrantly expressed genes, including 988 upregulated and 7495 downregulated genes, in IPF tissues compared to normal lung tissues (Table 1). There were 29 overlapping genes between the GSE2052, GSE24206, and GSE44723 datasets according to the Venn diagram, including 29 overlapping genes between GSE2052 and GSE44723, 268 overlapping genes between GSE44723 and GSE24206, and 389 overlapping genes between GSE2052 and GSE24206 (Figure 1(d)).Figure 1
Identification of DEGs from three GEO datasets (GSE2052, GSE44723, and GSE24206). (a–c) Volcano plot of the distributions of DEGs from the three datasets, mapping upregulated genes (red dots) and downregulated genes (blue dots). No significantly changed genes are marked as gray dots. (d) Identification of 29 commonly changed DEGs using Venn diagram from these three datasets (http://www.ehbio.com/ImageGP/index.php/Home/). Differently colored areas represented different datasets. The overlap areas meant the commonly changed DEGs. DEGs were identified with a paired t-test. The P value of less than 0.05 and [logFC] value of more than 1 were considered statistically significant.
(a)(b)(c)(d)Table 1
DEGs were identified from three datasets, including 262 upregulated and 191 downregulated genes in IPF compared with normal lung tissues (upregulated genes were defined withfoldchangeFC>0 as the cutoff criterion. The opposite was defined as downregulated genes).
DEGsGene namesUpregulatedSULF1, DEAF1, SCG5, DSG2, SLC1A4, CCND2, KCNN4, ST6GAL1, SLC38A1, SEC11C, XPOT, DPY30, PFKP, DDB1, HEPH, CXCL13, ATXN10, SEL1L3, CIAO1, CCL19, STARD5, SLFN12, ROR2, VAT1L, FPGT, GMPPA, COL18A1, COL7A1, PIGF, LMO4, FAM120A, SLC29A3, TCFL5, IGFBP2, UQCRQ, CCNA2, TWSG1, TCTN3, ASPN, PAM, BPIFB1, FAIM, FBLN2, SCARA3, COMP, ABCC5, DIO2, CHEK2, MCM4, TM9SF2, NAB1, DGKA, PTGFRN, FAT1, DOK5, CNIH1, ACTN1, PLA2G12A, MAGED1, ALG1, TWIST1, TRIM5, RCN2, CXCL14, ARMC1, STMN3, HMCN1, WDR5B, CROT, LEF1, TMEM14A, PLA2G4A, FKBP10, ABCC3, SPR, ROBO1, OXR1, CRLF1, TRIAP1, KDELR3, DIRAS3, BBS2, TGFB3, LGMN, CDK2AP1, CXCL12, RRM2, STRBP, TSPAN6, DAP, COL6A3, FZD6, TDO2, GMDS, PPP2R5E, SUPT7L, ZKSCAN7, CDKN3, CNTNAP1, IGF1, GSS, LRRC8D, TMEM98, TRIM2, LTBP1, BACE2, BRD8, COLEC11, FXN, PAFAH1B3, PGAM1, COL5A2, AOC1, ANTXR1, TMEM69, IMPACT, NET1, MXRA5, RCN1, MYO1E, DUSP23, CDH3, RHOD, CYP2S1, POSTN, ICMT, PDLIM4, C1QTNF6, ACVR1, SYTL2, PLA2G7, MFAP2, ZDHHC13, TMED10, ALDH1A3, SLN, CPOX, CDKN2C, PPIC, XRCC5, CLDN1, NSG1, ITGA7, R3HDM1, ERGIC2, TRIM36, EYA2, RPL39L, CCL13, RBP5, DONSON, SERPINB5, TXNDC15, HOMER3, ARL1, UBE2E3, CRYM, MEOX1, TMEM45A, COL15A1, ATP1B1, LDLRAD4, STEAP3, NABP2, BDKRB2, DCLK1, CFH, TRO, ECM1, PFN2, IL13RA2, MYOF, FHL2, CADPS, ITGAV, PCNA, PTK7, KIF2C, MEGF8, BMP4, PDCD2L, PRMT6, TP53BP1, OSBPL6, FMO1, PDE1A, PBX3, ELOVL4, ATF7IP, SYNDIG1, TMEM158, CFI, ALDH3A1, CKAP2, MRPL2, COL14A1, EGFL6, LHX6, THBS2, RRM1, YLPM1, TM7SF3, MLEC, CFB, BCL11A, GPR87, ZNF436, CLNS1A, ATIC, LGR4, CYP24A1, SEMA3C, PDGFC, TP63, ARMCX2, NUSAP1, ASB2, SLC39A6DownregulatedMKLN1, ECHDC3, RAB32, SLC25A51, HOPX, MATN3, MAP2, ARHGAP6, EPB41L5, NRGN, HEY1, PCTP, ACVRL1, TBX5, ERMP1, NAGA, MPP1, TXNIP, LRRN3, FLOT1AATK, RCL1, CSF3R, ANXA3, TEK, GRK5, HES1, HSPA1L, GATA6, EMP2, SLCO2A1, PMM1, STARD13, SEC14L1, SPTBN1, GHRL, TSPAN7, NEBL, ZNF655, TMEM11, UIMC1, NCOA3, ZFP36, CREBBP, LDLR, RAB20, SERTAD1, PPFIBP1, REPS2, ELF1, CALCOCO2, CSRNP1, GPM6A, DLL4, FPR1, CARHSP1, ADARB1, LMO7, RCOR1, LRRC32, CTNNBIP1, CA4, PADI4, OSGIN1, CXCL2, EGFR, CHI3L2, LPIN2, ANKHD1, ARHGEF4, ARHGAP29, VAMP5, RAI2, CYTIP, PRX, IER5, DNM2, IL17RA, BHLHE40, SLC39A8, RAPGEF5, PTPRM, CNOT8, DLC1, TLK2, EPAS1, PRELP, MAFF, ABTB1, HSD17B6, NDEL1, HYAL1, HECA, HSPB8, DNAJB1, CDH13, RGS16, PTPN12, CD55, TIPARP, CRYAB, CD36, NUP153, PTPRB, ITSN2, TNNC1, MAPT, THBD, CDKN2D, AOC3, P2RY1, ZBTB16, CSF3, EDNRB, FAM167A, SRGAP2, SLCO1A2, DAPK2, AGTR1, RIMKLB, ASRGL1, ANG, CCK, BCL2L13, OSGIN2, ACSM5, KIAA0040, KDR, FUT1, DOCK9, GADD45G, CLDN5, LIFR, STXBP6, GPR4, S1PR1, SLC1A1, PLAG1, EDA, DENND3, IDI1, KHDRBS3, CLEC1A, INMT, MPP3, PLLP, MTSS1, FSTL3, CRTAC1, GTF2IRD1, F3, KLF10, KLRD1, FBLN5IZUMO4, PIR, MAOA, C1QL1, THRB, RNF182, ALDH6A1, FAM49A, ST6GALNAC3, SSFA2, SLC25A24, AMPH, ADAMTS8, PLSCR1, BCKDHA, STXBP4, FLRT3, AOX1, SYAP1, RLF, SSH2, DERA, PIM1, STARD3NL, SUN2, SEPP1, IL1R2, EIF2A, FAH, METTL7A, EIF4E3, CHRM2, 1-Mar, PDK1, TJP2, RASL11A, NKX3-1


### 3.2. GO and KEGG Enrichment Analyses of Differentially Expressed Genes

The biological processes associated with the DEGs were determined by using the DAVID online bioinformatics database. As shown in Table2, the top 6 GO results revealed that the significantly enriched BPs of IPF DEGs were mainly concentrated in the cell adhesion, biological adhesion, regulation of cell proliferation, and so on. The top 6 significantly enriched MFs were mainly concentrated in the calcium ion binding, cytokine binding, chemokine activity, and so on. The top 6 significantly enriched CCs were mainly concentrated in the extracellular region part, extracellular space, extracellular matrix, and so on. The top 6 significantly enriched KEGG pathways were mainly concentrated in the ECM-receptor interaction, cytokine-cytokine receptor interaction, and so on.Table 2
The significantly enriched analysis of differentially expressed genes in idiopathic pulmonary fibrosis.
CategoryPathway IDDescriptionCountP valueGOTERM_BPGO:0007155Cell adhesion451.46327E−09GOTERM_BPGO:0022610Biological adhesion451.50199E−09GOTERM_BPGO:0042127Regulation of cell proliferation429.51714E−07GOTERM_BPGO:0035295Tube development181.51015E−05GOTERM_BPGO:0032103Positive regulation of response to external stimulus101.56397E−05GOTERM_BPGO:0001501Skeletal system development221.77224E−05GOTERM_MFGO:0005509Calcium ion binding380.000949643GOTERM_MFGO:0019955Cytokine binding100.001081524GOTERM_MFGO:0008009Chemokine activity60.004374834GOTERM_MFGO:0042802Identical protein binding270.00485805GOTERM_MFGO:0042379Chemokine receptor binding60.005747644GOTERM_MFGO:0008017Microtubule binding70.00692901GOTERM_CCGO:0044421Extracellular region part543.94172E−10GOTERM_CCGO:0005615Extracellular space384.31806E−07GOTERM_CCGO:0031012Extracellular matrix242.55104E−06GOTERM_CCGO:0005578Proteinaceous extracellular matrix228.9433E−06GOTERM_CCGO:0044459Plasma membrane part765.14236E−05GOTERM_CCGO:0031226Intrinsic to plasma membrane480.000103242KEGG_PATHWAYhsa04512ECM-receptor interaction90.00286045KEGG_PATHWAYhsa04060Cytokine-cytokine receptor interaction170.003658523KEGG_PATHWAYhsa04510Focal adhesion130.01357663KEGG_PATHWAYhsa04610Complement and coagulation cascades70.014586106KEGG_PATHWAYhsa05414Dilated cardiomyopathy80.017081193KEGG_PATHWAYhsa00360Phenylalanine metabolism40.024804532


### 3.3. Construction and Enrichment Analysis of Modules

The IPF DEGs were used to construct the PPI network by using the STRING online database, and a total of 18 modules were obtained with the ClusterONE Cytoscape plug-in (Figure2). To obtain functional and pathway enrichment information, the genes involved in these 18 modules were further analyzed by using DAVID as shown in Figures 3(a)–3(d). The top 10 modules of significantly enriched BPs were mainly concentrated in the protein polyubiquitination, ciliary basal body-plasma membrane docking, Golgi vesicle transport, and so on. The top 10 modules of significantly enriched CCs were mainly concentrated in the ubiquitin ligase complex, microtubule organizing center part, microtubule-associated complex, and so on. The top 10 modules of significantly enriched MFs were mainly concentrated in the ubiquitin-protein transferase activity, structural constituents of the cytoskeleton, microtubule motor activity, and so on (Table 3). The KEGG pathways of these 18 modules were concentrated in the ubiquitin-mediated proteolysis, spliceosome, purine metabolism, Fanconi anemia pathway, and so on (Table 4).Figure 2
The module genes were filtered into the DEG protein-protein interaction (PPI) network complex.Figure 3
Gene Ontology and KEGG analysis of module gene function in idiopathic pulmonary fibrosis. (a–c) Gene Oncology (GO) analysis was conducted to identify overrepresented GO terms for the significantly enriched analysis of module genes. (d) KEGG analysis was conducted to identify the biological information of module genes.
(a)(b)(c)(d)Table 3
The significant GO enrichment analysis of module gene function in idiopathic pulmonary fibrosis.
CategoryPathway IDPathway descriptionCountP valueGO_BP_m1GO:0000209Protein polyubiquitination471.77E−64GO_BP_m10GO:0097711Ciliary basal body-plasma membrane docking326.55E−71GO_BP_m11GO:0048193Golgi vesicle transport309.70E−41GO_BP_m12GO:0006898Receptor-mediated endocytosis196.36E−21GO_BP_m13GO:0060337Type I interferon signaling pathway201.37E−38GO_BP_m14GO:0007059Chromosome segregation292.58E−38GO_BP_m15GO:0009165Nucleotide biosynthetic process202.61E−27GO_BP_m16GO:0036297Interstrand cross-link repair161.73E−32GO_BP_m17GO:0006289Nucleotide-excision repair267.46E−53GO_BP_m18GO:0000280Nuclear division247.63E−35GO_CC_m1GO:0000151Ubiquitin ligase complex421.34E−54GO_CC_m10GO:0044450Microtubule organizing center part191.93E−29GO_CC_m11GO:0005875Microtubule-associated complex271.16E−46GO_CC_m12GO:0030136Clathrin-coated vesicle228.55E−33GO_CC_m13GO:0042611MHC protein complex56.87E−10GO_CC_m14GO:0000775Chromosome, centromeric region354.74E−59GO_CC_m15GO:0008074Guanylate cyclase complex, soluble22.57E−06GO_CC_m16GO:0043240Fanconi anemia nuclear complex79.55E−17GO_CC_m17GO:1990391DNA repair complex93.70E−17GO_CC_m18GO:0005819Spindle175.08E−23GO_MF_m1GO:0004842Ubiquitin-protein transferase activity695.10E−99GO_MF_m10GO:0005200Structural constituent of cytoskeleton64.58E−08GO_MF_m11GO:0003777Microtubule motor activity213.95E−38GO_MF_m12GO:0030276Clathrin binding92.00E−14GO_MF_m13GO:0042605Peptide antigen binding52.31E−09GO_MF_m14GO:0043515Kinetochore binding49.27E−10GO_MF_m15GO:0016776Phosphotransferase activity104.97E−20GO_MF_m16GO:0140097Catalytic activity, acting on DNA135.26E−17GO_MF_m17GO:0003684Damaged DNA binding169.67E−31GO_MF_m18GO:0004674Protein serine/threonine kinase activity86.57E−07Table 4
The significant KEGG enrichment analysis of module gene function in idiopathic pulmonary fibrosis.
PathwayPathway descriptionCountP valueGenesKEGG_Pathway_m1Ubiquitin mediated proteolysis301.21E−4610273, 10477, 11065, 22954, 23327, 25898, 51433, 51434, 54926, 55070, 57154, 65264, 7318, 7320, 7321, 7323, 7326, 7328, 7332, 7428, 83737, 8454, 8697, 9021, 90293, 9246, 92912, 9820, 991, 6921KEGG_Pathway_m2Spliceosome301.23E−4410084, 10262, 10285, 10594, 10907, 10946, 10992, 1665, 22827, 51340, 51690, 57187, 57819, 6428, 6429, 6430, 6432, 6434, 6625, 6626, 6628, 6635, 7307, 8175, 84321, 8449, 84991, 9775, 988, 9984KEGG_Pathway_m15Purine metabolism268.05E−3910201, 11128, 124583, 171568, 205, 2987, 29922, 3704, 4831, 4832, 4833, 4881, 50484, 50808, 51251, 5138, 5315, 5422, 6240, 6241, 84284, 953, 955, 956, 2982, 2983KEGG_Pathway_m16Fanconi anemia pathway202.23E−372067, 2072, 2176, 2178, 2188, 2189, 22909, 29089, 51426, 5429, 55120, 5889, 6118, 6119, 672, 80010, 80198, 84464, 80233, 91442KEGG_Pathway_m4Protein digestion and absorption211.16E−361277, 1278, 1281, 1285, 1287, 1288, 1289, 1290, 1292, 1293, 1294, 1299, 1300, 1303, 1306, 1308, 255631, 7373, 80781, 81578, 85301KEGG_Pathway_m5Glutathione metabolism201.44E−36124975, 2678, 2687, 27306, 2878, 2879, 2880, 2882, 2937, 2938, 2941, 2948, 2949, 2950, 2952, 373156, 4257, 4258, 493869, 51060KEGG_Pathway_m6Ribosome204.77E−3511222, 29088, 51069, 51073, 51116, 51263, 51264, 51318, 54460, 54948, 6183, 63875, 63931, 64928, 64963, 64965, 64983, 65003, 65008, 79590KEGG_Pathway_m17Nucleotide excision repair164.27E−301022, 1069, 1642, 1643, 2067, 2072, 2968, 3978, 5111, 5425, 6118, 6119, 7507, 7508, 8451, 902KEGG_Pathway_m9Neuroactive ligand-receptor interaction257.26E−2810161, 10800, 1131, 1241, 148, 154, 185, 1902, 1910, 2149, 2151, 2925, 3061, 3357, 3358, 4829, 4923, 5021, 5028, 5031, 624, 680, 6915, 7201, 9002KEGG_Pathway_m3Neuroactive ligand-receptor interaction253.41E−211129, 1268, 150, 152, 1813, 187, 1901, 1902, 1903, 2357, 2358, 2359, 2587, 2913, 2918, 4543, 4887, 59340, 624, 6752, 6755, 719, 728, 9294, 9568KEGG_Pathway_m8Ribosome biogenesis in eukaryotes125.45E−2010171, 10799, 134430, 23560, 27341, 3692, 51096, 51119, 51602, 55341, 55651, 84916KEGG_Pathway_m7Ribosome177.99E−1911224, 25873, 51065, 6128, 6141, 6156, 6157, 6160, 6168, 6169, 6181, 6201, 6205, 6229, 6231, 6232, 6234KEGG_Pathway_m12Endocytosis163.91E−1610109, 1173, 1759, 1785, 1956, 22905, 26119, 27131, 273, 3949, 408, 5868, 6456, 6643, 867, 8976KEGG_Pathway_m13Herpes simplex infection118.93E−1310379, 3105, 3113, 3115, 3122, 3134, 3434, 3661, 4938, 4940, 6041KEGG_Pathway_m18Cell cycle79.18E−114085, 701, 7272, 8379, 890, 891, 9133KEGG_Pathway_m14Oocyte meiosis88.35E−104085, 5516, 5525, 5528, 5529, 8379, 891, 9133KEGG_Pathway_m11Vasopressin-regulated water reabsorption61.07E−0910540, 1639, 51164, 79659, 84516, 8655KEGG_Pathway_m10Pathogenic Escherichia coli infection35.99E−05203068, 7277, 7846


### 3.4. Selection and Analysis of Key Genes

The biological network of differentially expressed IPF genes was constructed by using the BiNGO plug-in of Cytoscape, and the results revealed that most of the DEGs were significantly enriched in mitochondrial translation, cellular macromolecule metabolic process, cellular process, and so on (Figure4(a)). ClueGO, another plug-in of Cytoscape, can annotate and visualize the pathway networks of DEGs integrating GO terms as well as KEGG pathways. The results from ClueGO revealed that most of the DEGs were significantly enriched in the glutathione metabolism, Fanconi anemia pathway, etc. (Figure 4(b)).Figure 4
Identification and biological analysis of key genes in idiopathic pulmonary fibrosis. (a) Using BiNGO plug-in of Cytoscape, the biological network of the significant DEGs was conducted. (b) Using ClueGO plug-in of Cytoscape, GO and KEGG analyses were conducted to identify the significant DEGs. (c) The ncRNAs in module genes were identified using the hypergeometric test. (d) The TFs in module genes were identified using the hypergeometric test.
(a)(b)(c)(d)Subsequently, the key genes were obtained through calculation of the hypergeometric test. The 30 miRNAs and 4 lncRNAs enriched in 13 modules and the 44 transcription factors (TFs) enriched in 10 modules are presented in (Figures4(c) and 4(d)). According to the enrichment scores, the corresponding relevant noncoding RNAs (ncRNAs) were closely associated with ubiquitin-mediated proteolysis module m1, spliceosome module m2, cell cycle modules m14 and m18, and endocytosis module m12, which included long noncoding RNAs (lncRNAs) MALAT1 (modulelinks=14, P=7.6∗10−3), FENDRR (modulelinks=18, P=2.5∗10−3), RNU1-1 (modulelinks=23, P=0), and TUG1 (modulelinks=17, P=4.03∗10−7). The transcription factors (TFs) identified based on the enrichment scores were closely associated with GPR signaling pathway module m3, ECM-receptor interaction module m4, glutathione metabolism module m5, neuroactive ligand-receptor interaction module m9, endocytosis module m12, cell adhesion module m13, nucleotide excision repair module m17, homologous recombination module m16, and cell cycle modules m14 and m18, which included E2F1 (modulelinks=5, P=3∗10−4), TP53 (modulelinks=6, P=2∗10−4), YBX1 (modulelinks=4, P=1.24∗10−5), E2F4 (modulelinks=3, P=2∗10−4), SP1 (modulelinks=7, P=4.3∗10−3), BRCA1 (modulelinks=3, P=2.8∗10−3), CREB1 (modulelinks=5, P=4.74∗10−5), and CIITA (modulelinks=5, P=4.19∗10−7). As shown in Table 5, these key genes, such as MALAT1, RNU1-1, FENDRR, TUG1, E2F1, TP53, SP1, YBX1, BRCA1, E2F4, CREB1, and CIITA, play significant functional roles in their associated modules, suggesting that these genes may play roles in cell cycle regulation, methylation, acetyltransferase activity, and the splicing cycle. According to the integrated analysis results, these key genes of lncRNAs and TFs might play pathogenic roles in the occurrence and progression of IPF.Table 5
Functional roles of the strongest correlations of lncRNAs and TFs in IPF.
No.Gene symbolFull nameFunction1MALAT1Metastasis-associated lung adenocarcinoma transcript 1Form molecular scaffolds for ribonucleoprotein complexes, acting as a transcriptional regulator for numerous genes, and involved in cell cycle regulation2RNU1-1RNA, U1 small nuclear 1Its related pathways are spliceosomal splicing cycle3FENDRRFOXF1 adjacent noncoding developmental regulatory RNABind to polycomb repressive complex 2 and/or TrxG/MLL complexes to promote the methylation of the promoters of target genes4TUG1Taurine upregulated 1Interacts with the polycomb repressor complex and functions in the epigenetic regulation of transcription, acting as a sponge for microRNAs5E2F1E2F transcription factor 1Bind preferentially to retinoblastoma protein pRB in a cell cycle-dependent manner and mediate both cell proliferation and p53-dependent/independent apoptosis6TP53Tumor protein P53Regulate expression of target genes, inducing cell cycle arrest, apoptosis, senescence, DNA repair, or changes in metabolism7SP1Sp1 transcription factorBe involved in many cellular processes, including cell differentiation, cell growth, apoptosis, immune responses, response to DNA damage, and chromatin remodeling8YBX1Y-box binding protein 1Be implicated in numerous cellular processes including regulation of transcription and translation, pre-mRNA splicing, DNA reparation, and mRNA packaging9BRCA1Breast cancer type 1 susceptibility proteinPlay a role in transcription, DNA repair of double-stranded breaks, and recombination which forms a large multisubunit protein complex known as the BRCA1-associated genome surveillance complex and interacts with histone deacetylase complexes10E2F4E2F transcription factor 4Act as proliferation-associated suppression genes and bind to all three of the tumor suppressor proteins pRB, p107, and p13011CREB1CAMP responsive element binding protein 1Its related pathways are development of HGF signaling pathway and circadian entrainment12CIITAClass II major histocompatibility complex transactivatorOnce it does not bind DNA but rather uses an intrinsic acetyltransferase activity to act in a coactivator-like fashion


### 3.5. Validation of the lncRNAs and TFs in IPF with qRT-PCR

Demographic and clinical features of the IPF patients and the healthy control group are listed in Table6. Patients with IPF smoked fewer cigarettes than in the control group. Moreover, the clinical features of the IPF group were decreased in pulmonary function. To validate the results obtained through integrated analysis of the three datasets related to IPF, the relative expression of the key genes was analyzed by RT-PCR (Fig. S4). We found that 3 of the 4 candidate genes have statistically significant differences between the IPF and normal groups (MALAT1, E2F1, and YBX1 with P<0.01, FENDRR with P>0.05).Table 6
Demographic and clinical characteristics of the subjects.
Healthy controlsIPFP valuen=20n=20Age, mean (SD)68.4 (4.7)69.0 (6.6)0.74Gender,n (%)0.07Male11 (55)12 (60)Female9 (45)8 (40)Smoker,n (%)0.001Current4 (20)2 (10)Former7 (35)10 (50)Never9 (45)8 (40)Smoking dose (pack-year)37 (14.5-58)35.7 (12-56)0.14FVC (% predicted), mean (SD)96.2 (11.7)63.88 (17.0)<0.01DLCO (% predicted), mean (SD)80.8 (17.6)42.04 (16.5)<0.01FEV1 (% predicted), mean (SD)106.0 (18.6)71.7 (11.0)<0.001FEV1/FVC (%), mean (SD)79.9 (5.0)57.8 (8.3)<0.001FVC: forced vital capacity; DLCO: diffusing capacity of carbon monoxide; FEV1: forced expiratory volume in one second.


## 3.1. Identification of Differentially Expressed Genes in IPF

After normalization and standardization of the raw data from these three GSE2052, GSE44723, and GSE24206 datasets (Figures1(a)–1(c)), we identified a total of 8483 aberrantly expressed genes, including 988 upregulated and 7495 downregulated genes, in IPF tissues compared to normal lung tissues (Table 1). There were 29 overlapping genes between the GSE2052, GSE24206, and GSE44723 datasets according to the Venn diagram, including 29 overlapping genes between GSE2052 and GSE44723, 268 overlapping genes between GSE44723 and GSE24206, and 389 overlapping genes between GSE2052 and GSE24206 (Figure 1(d)).Figure 1
Identification of DEGs from three GEO datasets (GSE2052, GSE44723, and GSE24206). (a–c) Volcano plot of the distributions of DEGs from the three datasets, mapping upregulated genes (red dots) and downregulated genes (blue dots). No significantly changed genes are marked as gray dots. (d) Identification of 29 commonly changed DEGs using Venn diagram from these three datasets (http://www.ehbio.com/ImageGP/index.php/Home/). Differently colored areas represented different datasets. The overlap areas meant the commonly changed DEGs. DEGs were identified with a paired t-test. The P value of less than 0.05 and [logFC] value of more than 1 were considered statistically significant.
(a)(b)(c)(d)Table 1
DEGs were identified from three datasets, including 262 upregulated and 191 downregulated genes in IPF compared with normal lung tissues (upregulated genes were defined withfoldchangeFC>0 as the cutoff criterion. The opposite was defined as downregulated genes).
DEGsGene namesUpregulatedSULF1, DEAF1, SCG5, DSG2, SLC1A4, CCND2, KCNN4, ST6GAL1, SLC38A1, SEC11C, XPOT, DPY30, PFKP, DDB1, HEPH, CXCL13, ATXN10, SEL1L3, CIAO1, CCL19, STARD5, SLFN12, ROR2, VAT1L, FPGT, GMPPA, COL18A1, COL7A1, PIGF, LMO4, FAM120A, SLC29A3, TCFL5, IGFBP2, UQCRQ, CCNA2, TWSG1, TCTN3, ASPN, PAM, BPIFB1, FAIM, FBLN2, SCARA3, COMP, ABCC5, DIO2, CHEK2, MCM4, TM9SF2, NAB1, DGKA, PTGFRN, FAT1, DOK5, CNIH1, ACTN1, PLA2G12A, MAGED1, ALG1, TWIST1, TRIM5, RCN2, CXCL14, ARMC1, STMN3, HMCN1, WDR5B, CROT, LEF1, TMEM14A, PLA2G4A, FKBP10, ABCC3, SPR, ROBO1, OXR1, CRLF1, TRIAP1, KDELR3, DIRAS3, BBS2, TGFB3, LGMN, CDK2AP1, CXCL12, RRM2, STRBP, TSPAN6, DAP, COL6A3, FZD6, TDO2, GMDS, PPP2R5E, SUPT7L, ZKSCAN7, CDKN3, CNTNAP1, IGF1, GSS, LRRC8D, TMEM98, TRIM2, LTBP1, BACE2, BRD8, COLEC11, FXN, PAFAH1B3, PGAM1, COL5A2, AOC1, ANTXR1, TMEM69, IMPACT, NET1, MXRA5, RCN1, MYO1E, DUSP23, CDH3, RHOD, CYP2S1, POSTN, ICMT, PDLIM4, C1QTNF6, ACVR1, SYTL2, PLA2G7, MFAP2, ZDHHC13, TMED10, ALDH1A3, SLN, CPOX, CDKN2C, PPIC, XRCC5, CLDN1, NSG1, ITGA7, R3HDM1, ERGIC2, TRIM36, EYA2, RPL39L, CCL13, RBP5, DONSON, SERPINB5, TXNDC15, HOMER3, ARL1, UBE2E3, CRYM, MEOX1, TMEM45A, COL15A1, ATP1B1, LDLRAD4, STEAP3, NABP2, BDKRB2, DCLK1, CFH, TRO, ECM1, PFN2, IL13RA2, MYOF, FHL2, CADPS, ITGAV, PCNA, PTK7, KIF2C, MEGF8, BMP4, PDCD2L, PRMT6, TP53BP1, OSBPL6, FMO1, PDE1A, PBX3, ELOVL4, ATF7IP, SYNDIG1, TMEM158, CFI, ALDH3A1, CKAP2, MRPL2, COL14A1, EGFL6, LHX6, THBS2, RRM1, YLPM1, TM7SF3, MLEC, CFB, BCL11A, GPR87, ZNF436, CLNS1A, ATIC, LGR4, CYP24A1, SEMA3C, PDGFC, TP63, ARMCX2, NUSAP1, ASB2, SLC39A6DownregulatedMKLN1, ECHDC3, RAB32, SLC25A51, HOPX, MATN3, MAP2, ARHGAP6, EPB41L5, NRGN, HEY1, PCTP, ACVRL1, TBX5, ERMP1, NAGA, MPP1, TXNIP, LRRN3, FLOT1AATK, RCL1, CSF3R, ANXA3, TEK, GRK5, HES1, HSPA1L, GATA6, EMP2, SLCO2A1, PMM1, STARD13, SEC14L1, SPTBN1, GHRL, TSPAN7, NEBL, ZNF655, TMEM11, UIMC1, NCOA3, ZFP36, CREBBP, LDLR, RAB20, SERTAD1, PPFIBP1, REPS2, ELF1, CALCOCO2, CSRNP1, GPM6A, DLL4, FPR1, CARHSP1, ADARB1, LMO7, RCOR1, LRRC32, CTNNBIP1, CA4, PADI4, OSGIN1, CXCL2, EGFR, CHI3L2, LPIN2, ANKHD1, ARHGEF4, ARHGAP29, VAMP5, RAI2, CYTIP, PRX, IER5, DNM2, IL17RA, BHLHE40, SLC39A8, RAPGEF5, PTPRM, CNOT8, DLC1, TLK2, EPAS1, PRELP, MAFF, ABTB1, HSD17B6, NDEL1, HYAL1, HECA, HSPB8, DNAJB1, CDH13, RGS16, PTPN12, CD55, TIPARP, CRYAB, CD36, NUP153, PTPRB, ITSN2, TNNC1, MAPT, THBD, CDKN2D, AOC3, P2RY1, ZBTB16, CSF3, EDNRB, FAM167A, SRGAP2, SLCO1A2, DAPK2, AGTR1, RIMKLB, ASRGL1, ANG, CCK, BCL2L13, OSGIN2, ACSM5, KIAA0040, KDR, FUT1, DOCK9, GADD45G, CLDN5, LIFR, STXBP6, GPR4, S1PR1, SLC1A1, PLAG1, EDA, DENND3, IDI1, KHDRBS3, CLEC1A, INMT, MPP3, PLLP, MTSS1, FSTL3, CRTAC1, GTF2IRD1, F3, KLF10, KLRD1, FBLN5IZUMO4, PIR, MAOA, C1QL1, THRB, RNF182, ALDH6A1, FAM49A, ST6GALNAC3, SSFA2, SLC25A24, AMPH, ADAMTS8, PLSCR1, BCKDHA, STXBP4, FLRT3, AOX1, SYAP1, RLF, SSH2, DERA, PIM1, STARD3NL, SUN2, SEPP1, IL1R2, EIF2A, FAH, METTL7A, EIF4E3, CHRM2, 1-Mar, PDK1, TJP2, RASL11A, NKX3-1


## 3.2. GO and KEGG Enrichment Analyses of Differentially Expressed Genes

The biological processes associated with the DEGs were determined by using the DAVID online bioinformatics database. As shown in Table2, the top 6 GO results revealed that the significantly enriched BPs of IPF DEGs were mainly concentrated in the cell adhesion, biological adhesion, regulation of cell proliferation, and so on. The top 6 significantly enriched MFs were mainly concentrated in the calcium ion binding, cytokine binding, chemokine activity, and so on. The top 6 significantly enriched CCs were mainly concentrated in the extracellular region part, extracellular space, extracellular matrix, and so on. The top 6 significantly enriched KEGG pathways were mainly concentrated in the ECM-receptor interaction, cytokine-cytokine receptor interaction, and so on.Table 2
The significantly enriched analysis of differentially expressed genes in idiopathic pulmonary fibrosis.
CategoryPathway IDDescriptionCountP valueGOTERM_BPGO:0007155Cell adhesion451.46327E−09GOTERM_BPGO:0022610Biological adhesion451.50199E−09GOTERM_BPGO:0042127Regulation of cell proliferation429.51714E−07GOTERM_BPGO:0035295Tube development181.51015E−05GOTERM_BPGO:0032103Positive regulation of response to external stimulus101.56397E−05GOTERM_BPGO:0001501Skeletal system development221.77224E−05GOTERM_MFGO:0005509Calcium ion binding380.000949643GOTERM_MFGO:0019955Cytokine binding100.001081524GOTERM_MFGO:0008009Chemokine activity60.004374834GOTERM_MFGO:0042802Identical protein binding270.00485805GOTERM_MFGO:0042379Chemokine receptor binding60.005747644GOTERM_MFGO:0008017Microtubule binding70.00692901GOTERM_CCGO:0044421Extracellular region part543.94172E−10GOTERM_CCGO:0005615Extracellular space384.31806E−07GOTERM_CCGO:0031012Extracellular matrix242.55104E−06GOTERM_CCGO:0005578Proteinaceous extracellular matrix228.9433E−06GOTERM_CCGO:0044459Plasma membrane part765.14236E−05GOTERM_CCGO:0031226Intrinsic to plasma membrane480.000103242KEGG_PATHWAYhsa04512ECM-receptor interaction90.00286045KEGG_PATHWAYhsa04060Cytokine-cytokine receptor interaction170.003658523KEGG_PATHWAYhsa04510Focal adhesion130.01357663KEGG_PATHWAYhsa04610Complement and coagulation cascades70.014586106KEGG_PATHWAYhsa05414Dilated cardiomyopathy80.017081193KEGG_PATHWAYhsa00360Phenylalanine metabolism40.024804532


## 3.3. Construction and Enrichment Analysis of Modules

The IPF DEGs were used to construct the PPI network by using the STRING online database, and a total of 18 modules were obtained with the ClusterONE Cytoscape plug-in (Figure2). To obtain functional and pathway enrichment information, the genes involved in these 18 modules were further analyzed by using DAVID as shown in Figures 3(a)–3(d). The top 10 modules of significantly enriched BPs were mainly concentrated in the protein polyubiquitination, ciliary basal body-plasma membrane docking, Golgi vesicle transport, and so on. The top 10 modules of significantly enriched CCs were mainly concentrated in the ubiquitin ligase complex, microtubule organizing center part, microtubule-associated complex, and so on. The top 10 modules of significantly enriched MFs were mainly concentrated in the ubiquitin-protein transferase activity, structural constituents of the cytoskeleton, microtubule motor activity, and so on (Table 3). The KEGG pathways of these 18 modules were concentrated in the ubiquitin-mediated proteolysis, spliceosome, purine metabolism, Fanconi anemia pathway, and so on (Table 4).Figure 2
The module genes were filtered into the DEG protein-protein interaction (PPI) network complex.Figure 3
Gene Ontology and KEGG analysis of module gene function in idiopathic pulmonary fibrosis. (a–c) Gene Oncology (GO) analysis was conducted to identify overrepresented GO terms for the significantly enriched analysis of module genes. (d) KEGG analysis was conducted to identify the biological information of module genes.
(a)(b)(c)(d)Table 3
The significant GO enrichment analysis of module gene function in idiopathic pulmonary fibrosis.
CategoryPathway IDPathway descriptionCountP valueGO_BP_m1GO:0000209Protein polyubiquitination471.77E−64GO_BP_m10GO:0097711Ciliary basal body-plasma membrane docking326.55E−71GO_BP_m11GO:0048193Golgi vesicle transport309.70E−41GO_BP_m12GO:0006898Receptor-mediated endocytosis196.36E−21GO_BP_m13GO:0060337Type I interferon signaling pathway201.37E−38GO_BP_m14GO:0007059Chromosome segregation292.58E−38GO_BP_m15GO:0009165Nucleotide biosynthetic process202.61E−27GO_BP_m16GO:0036297Interstrand cross-link repair161.73E−32GO_BP_m17GO:0006289Nucleotide-excision repair267.46E−53GO_BP_m18GO:0000280Nuclear division247.63E−35GO_CC_m1GO:0000151Ubiquitin ligase complex421.34E−54GO_CC_m10GO:0044450Microtubule organizing center part191.93E−29GO_CC_m11GO:0005875Microtubule-associated complex271.16E−46GO_CC_m12GO:0030136Clathrin-coated vesicle228.55E−33GO_CC_m13GO:0042611MHC protein complex56.87E−10GO_CC_m14GO:0000775Chromosome, centromeric region354.74E−59GO_CC_m15GO:0008074Guanylate cyclase complex, soluble22.57E−06GO_CC_m16GO:0043240Fanconi anemia nuclear complex79.55E−17GO_CC_m17GO:1990391DNA repair complex93.70E−17GO_CC_m18GO:0005819Spindle175.08E−23GO_MF_m1GO:0004842Ubiquitin-protein transferase activity695.10E−99GO_MF_m10GO:0005200Structural constituent of cytoskeleton64.58E−08GO_MF_m11GO:0003777Microtubule motor activity213.95E−38GO_MF_m12GO:0030276Clathrin binding92.00E−14GO_MF_m13GO:0042605Peptide antigen binding52.31E−09GO_MF_m14GO:0043515Kinetochore binding49.27E−10GO_MF_m15GO:0016776Phosphotransferase activity104.97E−20GO_MF_m16GO:0140097Catalytic activity, acting on DNA135.26E−17GO_MF_m17GO:0003684Damaged DNA binding169.67E−31GO_MF_m18GO:0004674Protein serine/threonine kinase activity86.57E−07Table 4
The significant KEGG enrichment analysis of module gene function in idiopathic pulmonary fibrosis.
PathwayPathway descriptionCountP valueGenesKEGG_Pathway_m1Ubiquitin mediated proteolysis301.21E−4610273, 10477, 11065, 22954, 23327, 25898, 51433, 51434, 54926, 55070, 57154, 65264, 7318, 7320, 7321, 7323, 7326, 7328, 7332, 7428, 83737, 8454, 8697, 9021, 90293, 9246, 92912, 9820, 991, 6921KEGG_Pathway_m2Spliceosome301.23E−4410084, 10262, 10285, 10594, 10907, 10946, 10992, 1665, 22827, 51340, 51690, 57187, 57819, 6428, 6429, 6430, 6432, 6434, 6625, 6626, 6628, 6635, 7307, 8175, 84321, 8449, 84991, 9775, 988, 9984KEGG_Pathway_m15Purine metabolism268.05E−3910201, 11128, 124583, 171568, 205, 2987, 29922, 3704, 4831, 4832, 4833, 4881, 50484, 50808, 51251, 5138, 5315, 5422, 6240, 6241, 84284, 953, 955, 956, 2982, 2983KEGG_Pathway_m16Fanconi anemia pathway202.23E−372067, 2072, 2176, 2178, 2188, 2189, 22909, 29089, 51426, 5429, 55120, 5889, 6118, 6119, 672, 80010, 80198, 84464, 80233, 91442KEGG_Pathway_m4Protein digestion and absorption211.16E−361277, 1278, 1281, 1285, 1287, 1288, 1289, 1290, 1292, 1293, 1294, 1299, 1300, 1303, 1306, 1308, 255631, 7373, 80781, 81578, 85301KEGG_Pathway_m5Glutathione metabolism201.44E−36124975, 2678, 2687, 27306, 2878, 2879, 2880, 2882, 2937, 2938, 2941, 2948, 2949, 2950, 2952, 373156, 4257, 4258, 493869, 51060KEGG_Pathway_m6Ribosome204.77E−3511222, 29088, 51069, 51073, 51116, 51263, 51264, 51318, 54460, 54948, 6183, 63875, 63931, 64928, 64963, 64965, 64983, 65003, 65008, 79590KEGG_Pathway_m17Nucleotide excision repair164.27E−301022, 1069, 1642, 1643, 2067, 2072, 2968, 3978, 5111, 5425, 6118, 6119, 7507, 7508, 8451, 902KEGG_Pathway_m9Neuroactive ligand-receptor interaction257.26E−2810161, 10800, 1131, 1241, 148, 154, 185, 1902, 1910, 2149, 2151, 2925, 3061, 3357, 3358, 4829, 4923, 5021, 5028, 5031, 624, 680, 6915, 7201, 9002KEGG_Pathway_m3Neuroactive ligand-receptor interaction253.41E−211129, 1268, 150, 152, 1813, 187, 1901, 1902, 1903, 2357, 2358, 2359, 2587, 2913, 2918, 4543, 4887, 59340, 624, 6752, 6755, 719, 728, 9294, 9568KEGG_Pathway_m8Ribosome biogenesis in eukaryotes125.45E−2010171, 10799, 134430, 23560, 27341, 3692, 51096, 51119, 51602, 55341, 55651, 84916KEGG_Pathway_m7Ribosome177.99E−1911224, 25873, 51065, 6128, 6141, 6156, 6157, 6160, 6168, 6169, 6181, 6201, 6205, 6229, 6231, 6232, 6234KEGG_Pathway_m12Endocytosis163.91E−1610109, 1173, 1759, 1785, 1956, 22905, 26119, 27131, 273, 3949, 408, 5868, 6456, 6643, 867, 8976KEGG_Pathway_m13Herpes simplex infection118.93E−1310379, 3105, 3113, 3115, 3122, 3134, 3434, 3661, 4938, 4940, 6041KEGG_Pathway_m18Cell cycle79.18E−114085, 701, 7272, 8379, 890, 891, 9133KEGG_Pathway_m14Oocyte meiosis88.35E−104085, 5516, 5525, 5528, 5529, 8379, 891, 9133KEGG_Pathway_m11Vasopressin-regulated water reabsorption61.07E−0910540, 1639, 51164, 79659, 84516, 8655KEGG_Pathway_m10Pathogenic Escherichia coli infection35.99E−05203068, 7277, 7846


## 3.4. Selection and Analysis of Key Genes

The biological network of differentially expressed IPF genes was constructed by using the BiNGO plug-in of Cytoscape, and the results revealed that most of the DEGs were significantly enriched in mitochondrial translation, cellular macromolecule metabolic process, cellular process, and so on (Figure4(a)). ClueGO, another plug-in of Cytoscape, can annotate and visualize the pathway networks of DEGs integrating GO terms as well as KEGG pathways. The results from ClueGO revealed that most of the DEGs were significantly enriched in the glutathione metabolism, Fanconi anemia pathway, etc. (Figure 4(b)).Figure 4
Identification and biological analysis of key genes in idiopathic pulmonary fibrosis. (a) Using BiNGO plug-in of Cytoscape, the biological network of the significant DEGs was conducted. (b) Using ClueGO plug-in of Cytoscape, GO and KEGG analyses were conducted to identify the significant DEGs. (c) The ncRNAs in module genes were identified using the hypergeometric test. (d) The TFs in module genes were identified using the hypergeometric test.
(a)(b)(c)(d)Subsequently, the key genes were obtained through calculation of the hypergeometric test. The 30 miRNAs and 4 lncRNAs enriched in 13 modules and the 44 transcription factors (TFs) enriched in 10 modules are presented in (Figures4(c) and 4(d)). According to the enrichment scores, the corresponding relevant noncoding RNAs (ncRNAs) were closely associated with ubiquitin-mediated proteolysis module m1, spliceosome module m2, cell cycle modules m14 and m18, and endocytosis module m12, which included long noncoding RNAs (lncRNAs) MALAT1 (modulelinks=14, P=7.6∗10−3), FENDRR (modulelinks=18, P=2.5∗10−3), RNU1-1 (modulelinks=23, P=0), and TUG1 (modulelinks=17, P=4.03∗10−7). The transcription factors (TFs) identified based on the enrichment scores were closely associated with GPR signaling pathway module m3, ECM-receptor interaction module m4, glutathione metabolism module m5, neuroactive ligand-receptor interaction module m9, endocytosis module m12, cell adhesion module m13, nucleotide excision repair module m17, homologous recombination module m16, and cell cycle modules m14 and m18, which included E2F1 (modulelinks=5, P=3∗10−4), TP53 (modulelinks=6, P=2∗10−4), YBX1 (modulelinks=4, P=1.24∗10−5), E2F4 (modulelinks=3, P=2∗10−4), SP1 (modulelinks=7, P=4.3∗10−3), BRCA1 (modulelinks=3, P=2.8∗10−3), CREB1 (modulelinks=5, P=4.74∗10−5), and CIITA (modulelinks=5, P=4.19∗10−7). As shown in Table 5, these key genes, such as MALAT1, RNU1-1, FENDRR, TUG1, E2F1, TP53, SP1, YBX1, BRCA1, E2F4, CREB1, and CIITA, play significant functional roles in their associated modules, suggesting that these genes may play roles in cell cycle regulation, methylation, acetyltransferase activity, and the splicing cycle. According to the integrated analysis results, these key genes of lncRNAs and TFs might play pathogenic roles in the occurrence and progression of IPF.Table 5
Functional roles of the strongest correlations of lncRNAs and TFs in IPF.
No.Gene symbolFull nameFunction1MALAT1Metastasis-associated lung adenocarcinoma transcript 1Form molecular scaffolds for ribonucleoprotein complexes, acting as a transcriptional regulator for numerous genes, and involved in cell cycle regulation2RNU1-1RNA, U1 small nuclear 1Its related pathways are spliceosomal splicing cycle3FENDRRFOXF1 adjacent noncoding developmental regulatory RNABind to polycomb repressive complex 2 and/or TrxG/MLL complexes to promote the methylation of the promoters of target genes4TUG1Taurine upregulated 1Interacts with the polycomb repressor complex and functions in the epigenetic regulation of transcription, acting as a sponge for microRNAs5E2F1E2F transcription factor 1Bind preferentially to retinoblastoma protein pRB in a cell cycle-dependent manner and mediate both cell proliferation and p53-dependent/independent apoptosis6TP53Tumor protein P53Regulate expression of target genes, inducing cell cycle arrest, apoptosis, senescence, DNA repair, or changes in metabolism7SP1Sp1 transcription factorBe involved in many cellular processes, including cell differentiation, cell growth, apoptosis, immune responses, response to DNA damage, and chromatin remodeling8YBX1Y-box binding protein 1Be implicated in numerous cellular processes including regulation of transcription and translation, pre-mRNA splicing, DNA reparation, and mRNA packaging9BRCA1Breast cancer type 1 susceptibility proteinPlay a role in transcription, DNA repair of double-stranded breaks, and recombination which forms a large multisubunit protein complex known as the BRCA1-associated genome surveillance complex and interacts with histone deacetylase complexes10E2F4E2F transcription factor 4Act as proliferation-associated suppression genes and bind to all three of the tumor suppressor proteins pRB, p107, and p13011CREB1CAMP responsive element binding protein 1Its related pathways are development of HGF signaling pathway and circadian entrainment12CIITAClass II major histocompatibility complex transactivatorOnce it does not bind DNA but rather uses an intrinsic acetyltransferase activity to act in a coactivator-like fashion


## 3.5. Validation of the lncRNAs and TFs in IPF with qRT-PCR

Demographic and clinical features of the IPF patients and the healthy control group are listed in Table6. Patients with IPF smoked fewer cigarettes than in the control group. Moreover, the clinical features of the IPF group were decreased in pulmonary function. To validate the results obtained through integrated analysis of the three datasets related to IPF, the relative expression of the key genes was analyzed by RT-PCR (Fig. S4). We found that 3 of the 4 candidate genes have statistically significant differences between the IPF and normal groups (MALAT1, E2F1, and YBX1 with P<0.01, FENDRR with P>0.05).Table 6
Demographic and clinical characteristics of the subjects.
Healthy controlsIPFP valuen=20n=20Age, mean (SD)68.4 (4.7)69.0 (6.6)0.74Gender,n (%)0.07Male11 (55)12 (60)Female9 (45)8 (40)Smoker,n (%)0.001Current4 (20)2 (10)Former7 (35)10 (50)Never9 (45)8 (40)Smoking dose (pack-year)37 (14.5-58)35.7 (12-56)0.14FVC (% predicted), mean (SD)96.2 (11.7)63.88 (17.0)<0.01DLCO (% predicted), mean (SD)80.8 (17.6)42.04 (16.5)<0.01FEV1 (% predicted), mean (SD)106.0 (18.6)71.7 (11.0)<0.001FEV1/FVC (%), mean (SD)79.9 (5.0)57.8 (8.3)<0.001FVC: forced vital capacity; DLCO: diffusing capacity of carbon monoxide; FEV1: forced expiratory volume in one second.


## 4. Discussion

Chronic and progressive airway remodeling is a major characteristic of IPF with unknown etiology. Although accumulating evidence reveals that activated fibroblasts have important effects on the pathogenesis and progression of IPF, the underlying molecular mechanisms involved in the regulation of IPF remain unclear. Previous findings of gene regulation on IPF have mainly focused on protein-coding genes which can delay but do not inhabit the development of fibrosis. Recently, with the development of high-throughput sequencing technology, epigenetic researches provide new insights into the underlying molecular and etiological mechanisms of IPF. Epigenetics, such as functional ncRNAs, refers to heritable changes in DNA and chromatin that influence gene expression other than changes in DNA sequence and has gradually become the research hotspot. Multiple studies have indicated that lncRNAs can influence the pathological process involving the structural remodeling of pulmonary architecture and eventually lead to respiratory failure. As multifunctional adaptor molecules, lncRNAs play multifunctional roles in the regulation of gene expression by regulating mRNA decay, splicing, and gene looping by binding to DNA, proteins, and certain other RNAs [17, 18]. In this study, we integrated three publicly available microarray datasets (GSE2052, GSE44723, and GSE24206) and found that differential expression of 8483 genes comprised 988 upregulated and 7495 downregulated genes. Consistent with the results of previous studies on the molecular mechanism of IPF, we found that DEGs were mainly concentrated on the extracellular matrix and these biological functions were mainly related to cell adhesion, proliferation, cytoskeleton development, and cytokine interaction. After a series of bioinformatics analysis, the regulatory network consisting of key lncRNAs and transcription factors (TFs), which may contribute to the pathogenesis of IPF, was ultimately obtained. We found that the biological functions of these key genes, which were related to epithelial-mesenchymal transition (EMT), mainly focused on mitochondrial translation, RNA processing, and ubiquitin-mediated proteolysis. We performed a comprehensive literature search and judged by integrating degrees, closeness, and betweenness centrality of the regulatory network ultimately identifying 2 lncRNAs and 2 TFs (MALAT1, FENDRR and E2F1, YBX1, respectively). Subsequently, we further validate the expression levels of these key genes related to the regulation of pulmonary fibrosis in blood samples between the IPF and control groups using real-time polymerase chain reaction (RT-PCR). As a result, differential expression of these genes including downregulated YBX1 and upregulated MALAT1 and E2F1 reached statistical significance except for FENDRR between two groups.The research on epithelial-to-mesenchymal transition (EMT) related to a fibrotic process has received an increased attention in recent years. Considering the possible false-negative or false-positive results between the individual studies and the research sample size, we integrated and analyzed the potential lncRNAs and TFs related to the pathogenesis and progression of IPF from the public microarray data. Metastasis-associated lung adenocarcinoma transcript 1 (MALAT1) located on chromosome 11q13.1, also known as nuclear-enriched abundant transcript 2 (NEAT2), is involved in various biological functions including molecular scaffolds for ribonucleoprotein complexes, transcriptional regulator for genes, and regulation of cell cycle [19]. Substantial research has confirmed that MALAT1 has many important physiological and pathological function in a wide range of diseases such as various solid cancers, septic lung injury, myocardial or renal ischemia-reperfusion injury, cardiac fibrosis, liver fibrosis, and silica-induced pulmonary fibrosis [20–23]. Furthermore, MALAT1 also play an important role in EMT related to the pulmonary fibrosis [24]. Although MALAT1 was reported to be mainly localized in the nucleus, it could transfer from the nucleus to the cytoplasm during the G2/M-phase cell cycle [25]. E2F1 located on chromosome 20q11.2, which was screened out and validated in this study, belongs to TFs of the nuclear factor of the E2F family and participates in the cell cycle G1/S phase regulation mediating both cell proliferation and apoptosis [26]. Many studies have confirmed that E2F1 could activate the expression of stromal markers related to EMT such as vimentin and fibronectin and facilitate the pathogenesis processes such as fibrosis and tumor progression [27]. YBX1 located on chromosome 1p34.2, which is another screened and validated transcription factor in this study, belongs to a member of cold-shock protein family and acts as an important regulator related to cell proliferation and cell cycle [28]. Some studies have reported a correlation between abnormal expression of YBX1 and EMT markers such as vimentin and N-cadherin [29]. The above results suggest that genes screened and validated in this study might act as key regulators of the pathogenesis and progression of IPF.However, the limitations of this study are as follows. First, although key genes related to the pathogenesis of IPF have been screened and validated through integrating three datasets and performing a series of bioinformatics approaches, expression levels of these key genes need further validated experiments such as western blotting (WB) and immunohistochemistry analysis (IHC). Second, differential gene analysis is one of the crucial data analysis strategies for expression profiling of IPF in GEO datasets. However, the three datasets combined and analyzed in this study from the GEO database microarray and platform were not unified. Meanwhile, the sample sizes of these three datasets were relatively small and imbalanced. The potential selection bias and information bias were inevitable. Therefore, the accuracy and reliability of candidate genes could be improved greatly by integrating more various types of datasets. Third, verification of the expression levels of candidate genes in clinical samples is far from enough. Further functional verification of these candidate genes was necessary to perform by loss-of-function and gain-of-function experiments in vivo and in vitro. Lastly, the verification and discussion of the underlying molecular mechanisms of these candidate genes involved in the pathogenesis and progression of IPF will be necessary to confirm through a chromatin immunoprecipitation assay (CHIP) or dual-luciferase reporter gene assay and so on. Regardless of the limitations mentioned above, this study provided preliminary evidence for the candidate genes related to the pathogenesis of IPF. As a time-saving and cost-saving method for analysis of biomedical data, we took an extensive bioinformatics data-mining approach from different microarray platforms to obtain candidate lncRNAs and TFs in IPF. This study provided a framework and broad application prospects for exploring pathological molecular networks related to IPF.


## 5. Conclusion

This study provides reliable and comprehensive perspectives on the pathogenesis and progression of IPF; potential lncRNAs and TFs related to the pathogenesis of IPF were obtained through bioinformatics analysis. Ultimately, the 3 key genes that were found to show abnormal expression in IPF compared to normal lung tissues may be considered as biomarkers for the diagnosis and treatment of IPF, which should be verified in subsequent studies.



---

*Source: 6537462-2020-06-03.xml*