RNA sequencing and bioinformatics analysis of differentially expressed genes in the peripheral serum of ankylosing spondylitis patients

Patients and samples

A total of 18 patients with AS patients and 18 healthy age-matched controls were selected from the Affiliated Hospital of Qingdao University (Qingdao, China) between December 2020 and September 2021. All patients met the modified New York 1984 criteria [19] and were initially diagnosed with AS, drug-naive patients with short disease durations. None of the patients or controls had any previous history of cardiovascular disease, diabetes, hepatitis, malignancy, or other autoimmune and inflammatory illnesses.

Serum samples were collected using standard phlebotomy procedures and centrifuged at 3000 g for 10 min. The separated sera were stored in RNase-free centrifuge tubes at 80 °C until further processing. This research was approved and reviewed by the Medical Ethics Review Committee of the Affiliated Hospital of Qingdao University (approval number: QYFY WZLL 27251). All participants provided written informed consent in accordance with policies of the hospital ethics committee.

RNA extraction and sequencing

We randomly selected three patients with AS and three normal control (NC) for high-throughput sequencing [20]. Total RNA was extracted from serum using TRIzol reagent (Invitrogen; Thermo Fisher Scientific, Inc., Waltham, MA, USA) according to the manufacturer's instructions. Subsequently, the concentration and integrity of the isolated RNA were determined using the Qubit 3.0 Fluorometer (Invitrogen, Carlsbad, CA, USA) and an Agilent 2100 Bioanalyzer (Applied Biosystems, Carlsbad, CA, USA), respectively. RNA-seq libraries were prepared using the SMARTer Stranded Total RNA-Seq kit v.2 (Takara Bio USA, Mountain View, CA, USA) as previously [21]. The RNA samples were fragmented and reversely transcribed into first-strand cDNA, followed by second-strand synthesis. After cDNA synthesis, a tailing and adapter ligation was performed, and then the cDNA was amplified by PCR. Subsequently, the cDNA library quality and concentration were evaluated using the Agilent 2100 Bioanalyzer (Applied Biosystems, Carlsbad, CA, USA). The qPCR-based KAPA Biosystems Library Quantification kit (Kapa Biosystems, Inc.) was used for the quantification of the cDNA library. Ribosomal RNA depletion was performed during library construction according to the manufacturer’s protocol. Sequencing was carried out in a 150-bp paired-end run (PE150) using the NovaSeq 6000 system (Illumina, San Diego, CA, USA).

Analysis of DEGs

Reads were first mapped to the latest UCSC transcript set using Bowtie2 version 2.1.0 [22] and gene expression levels were estimated using RSEM v1.2.15 [23]. The trimmed mean of M-values was used to standardize gene expression. DEGs were then identified using edgeR software [24]. Genes showing altered expression with p < 0.05 and more than 1.5-fold changes were considered to be differentially expressed.

Gene ontology (GO), kyoto encyclopedia of genes and genomes (KEGG), and reactome analyses

To better understand the function and pathways of DEGs in AS, we performed GO functional annotation, KEGG enrichment [25], and Reactome analyses using the R package clusterProfiler [26]. GO analysis was used to investigate biological functions based on differentially expressed coding genes. This analysis classifies functions according to the three following aspects: biological process (BP), cellular component (CC), and molecular function (MF). A lower p-value indicated a higher significance of a GO term (p-value < 0.05). The KEGG and Reactome enrichment analyses were used to predict the related pathways of each DEG. A p-value of < 0.05 reflected significant enrichment.

Gene set enrichment analysis (GSEA)

GSEA was done as described in Subramanian et al. [27]. We used the R package fgsea to analyze the expression of filtered genes against MSigDB, a well-known molecular feature database. In addition, only two typical gene sets from MSigDB, H (hallmark gene sets) and C2:: CP (curated gene sets, canonical pathways), were analyzed by GSEA here. Finally, we retained results with statistical significance p-values < 0.05.

Construction of a protein–protein interaction (PPI) network and identification of hub genes and key modules

To gain insights into the correlation among DEGs at the protein level, the Search Tool for the Retrieval of Interacting Genes (STRING, https://stringdb.org) database was used to construct a PPI network of these DEGs [28]. The minimum required interaction score used to construct this PPI network was 0.4, and the isolated nodes were abandoned. Cytoscape [29] was used to visualize this PPI network and we used the plug-in Cytohubba to explore the hub genes of this PPI network [30]. Based on the centrality score, the key nodes in this PPI network were determined, and then the hub genes were deduced. Simultaneously, we used the molecular complex detection (MCODE) plugin for clustering analysis of gene networks to select the key subnetwork modules.

Validation of DEGs

Reverse transcription-quantitative polymerase chain reaction (qRT-PCR) experiments were conducted to validate the DEGs identified using high-throughput sequencing. cDNA was synthesized using the Prime Script RT reagent kit with genomic DNA eraser (TaKaRa, Tokyo, Japan), and qRT–PCR was performed on a LightCycler 480 (Roche, Indianapolis, IN, USA) using SYBR Green Master Mix (TaKaRa, Tokyo, Japan). The primer sequences are listed in Table 1. Validation experiments were performed using serum from 15 AS samples and 15 NC samples. mRNA levels of the selected new target genes were quantified by the 2−ΔΔCt method after normalization to the housekeeping gene GAPDH.

Table 1 Primers used in the present studyReceiver operating characteristic (ROC) analyses

To assess the diagnostic value of DEGs in AS, we performed ROC analyses using the pROC R package. We calculated the area under the curve (AUC) under the binomial exact confidence interval and Ggplot2 was applied for further visualization.

Statistical analysis

Statistical analysis was performed using SPSS 26.0 software (SPSS Inc., Chicago, IL, USA). All data are expressed as the mean ± SD. Statistical significance was determined by Student’s t test, and P values < 0.05 were considered to indicate a statistically significant difference.

留言 (0)

沒有登入
gif