Identification of key pathways and genes in carotid atherosclerosis through bioinformatics analysis of RNA-seq data

Basic analysis of sequencing data

After filtering the original data, checking the sequencing error rate, and verifying the distribution of GC content, the clean reads for follow-up analysis were obtained. Carotid atherosclerosis samples ASA_1 and ASA_2 generated 57.14, and 72.76 million clean reads, respectively. Normal control samples NA_1, NA_2, and NA_3 generated 59.84, 70.13, and 59.98 million clean reads, respectively (Table 1). Quality-controlled clean reads were compared to the human reference genome. Results of the comparison of each sample in this project with the human reference genome are presented in Table 2. Total mapped reads for ASA_1 and ASA_2 were 96.08 and 94.70%, whereas those for NA_1, NA_2, and NA_3 were 94.95, 96.28, and 95.99%, respectively.

Gene expression analysis

Box and density diagrams depicted uniform distribution of gene expression levels in different samples (Figure 2A–2B). Based on correlation analysis of gene expression levels among samples, the experiment was reliable with accurate sample selection (Figure 2C). Principal component analysis (PCA) revealed that samples were scattered among groups and clustered within groups (Figure 2D).

Figure 2. Analysis of DEGs. Boxplot (A) and density plot (B) showing that gene expression levels of different samples are evenly distributed. (C) Correlation diagram of gene expression level between samples indicates that the correlation coefficient between samples of the same group is high and the sample selection is reliable. (D) Principal component analysis diagram (PCA) showing that samples of different groups were significantly different, and samples within the same group were relatively uniform. (E) Heat map of DEGs, red represents high expression, blue represents low expression. (F) Volcano map of DEGs, screening parameters were |logFC|>3, adj.P.Val <0.05, 710 up-regulated genes (red dots), and 717 down-regulated genes (green dots) were identified, genes with no significant changes were labeled as black dots. ASA: Carotid atherosclerosis group. NA: Normal control group. FC: Fold change. Adj.p.Val: corrected p-value.

Identification of differentially expressed genes

After normalizing the sequencing data, we constructed the heat map (Figure 2E). A total of 1427 DEGs were identified after screening at |logFC|>3 and adj.p.Val <0.05 (Figure 2F, Supplementary Table 1); 710 DEGs and, 717 DEGs were up-regulated and down-regulated, respectively.

GO and KEGG enrichment analysis of DEGs

GO analysis revealed that biological processes (BP) of these targets were mainly enriched in leukocyte migration, regulation of activation, and inflammatory response. Enriched cell components (CC) included extracellular matrix, plasma membrane protein complex, and synaptic membrane. Molecular functions (MF) were mainly enriched in Metal ion Transmembrane Transporter activity, Substrate-specific channel activity, and amide binding (Figure 3A). KEGG pathway analysis demonstrated that enriched pathways for these targets are implicated in neuroactive ligand-receptor interaction, cytokine-cytokine receptor interaction, and cell inflammatory molecules (Figure 3B).

Figure 3. GO and KEGG enrichment analysis of DEGs. (A) GO enrichment analysis of DEGs, BP represents the biological process, CC represents the cellular component, and MF represents the molecular function. (B) KEGG pathway analysis of DEGs, adj.p.Val is corrected p-value, the count represents the number of DEGs. (C) Interaction between DEGs and top 5 terms of GO entries. (D) Interaction between DEGs and top 5 disease pathways in KEGG. DEGs: Differentially expressed genes. KEGG: Kyoto Encyclopedia of Genes and Genomes. FC: Fold change. Adj.p.Val: corrected p-value.

Cluster analysis for the top five GO terms explored the associations of genes such as MMP12, MMP8, and MMP9 with phagocytosis, leukocyte migration, acute response, regulation of immune activation, and regulation of immune effector process (Figure 3C). We performed cluster analysis for the top 5 KEGG pathways to explore the association of genes such as LBP, CSF3, and FCGR2A with leishmaniasis, hematopoietic cell, strain, tuberculosis, and osteoclast differentiation (Figure 3D).

Enrichment analysis of DEGs in MCODE

MCODE analysis demonstrated that DEGs were grouped into 18 modules. Briefly, Module 1 was comprised of 29 genes, including GAL, GNAl1, and OPRK1, etc. Module 2 comprised 11 genes, including TACR1, NTSR1, and TAC1, etc. Module 3 comprised 11 genes including CD74, CTSS, and CTSV, etc. (Figure 4). We did KEGG enrichment analysis on the 18 modules. Pathways in module 1 were mainly associated with neuroactive ligand-receptor interaction, whereas those in module 2 and module 3 were mainly enriched in cAMP signaling and chemokine signaling, respectively, and so on (Table 3).

Figure 4. MCODE analysis of DEGs. DEGs were divided into 18 modules based on functions, and protein-protein interaction networks were constructed for DEGs of each module. DEGs: Differentially expressed genes.

Constructing PPI network and screening for hub genes

The STRING online Database was employed to construct the PPI network of DEGs (Figure 5A), which was then visualized and optimized using Cytoscape software (Figure 5B). Analysis of hub genes was achieved using the degree method in cytoHubba; the top 20 genes were identified as hub genes (Figure 5C).

Figure 5. PPI network constructing and hub gene screening of DEGs. (A) PPI network constructed using STRING Online Database contained 1055 nodes and 1883 edges. (B) PPI network was optimized through Cytoscape software. (C) 20 hub genes identified through degree method in cytoHubba. PPI network: protein-protein interaction network.

Bioinformatics analysis of hub genes

Enrichment analysis of top 20 hub genes in the Metascape database revealed that these genes are mainly associated with Galpha (i) signaling events, GPCR ligand binding, and neuroactive ligand-receptor interaction (Figure 6A). TRRUST database analysis uncovered RELA and NF-κB1 to be the main transcription factors regulating the hub genes (Figure 6B). DisGeNET database disease enrichment analysis demonstrated that the hub genes were associated with depressive disorder, mental depression, and alcoholic intoxication (Figure 6C). Besides, PaGenBase database tissue characteristic enrichment analysis showed that hub genes were mainly enriched in the lung (Figure 6D). Analysis of the CMAP database showed that CeForanide, Chenodeoxycholic acid, and 0317956-0000 compounds (Table 4) potentially exerted a negative regulatory effect on hub genes. We deduced that these three drugs are potential candidates for the management of carotid atherosclerosis.

Figure 6. Bioinformatics analysis of hub genes. (A) Biological functions of hub genes analyzed through Metascape database. (B) Enrichment of transcriptional regulators of hub genes using the TRRUST database. (C) DisGeNET database enrichment analysis of diseases involving hub genes. (D) Tissue characteristics of hub genes as analyzed using PaGenBase database.

Patients and samples

We analyzed two CAS vessels of CAS patients and 3 carotid artery vessels of healthy subjects. All samples were collected from subjects who visited Liaocheng People's Hospital between February 2018 and August 2018. Each subject gave consent to undergo an electrocardiogram (ecg), echocardiogram, head, and chest CT, CTA, and other examinations. In the experimental group, blood vessels were acquired from CAS patients CAS with stroke and presenting with dizziness, verbal slurred speech, limb dysfunction, and other manifestations of cerebral ischemia. After admission, patients underwent craniocerebral magnetic resonance imaging, head and neck CTA, and whole-brain angiography, which conformed to the latest guide on diagnosis and treatment of atherosclerosis. CAS plaques and intima of vessels were obtained from patients diagnosed with CAS, particularly those who accepted treatment via carotid endarterectomy. These patients were categorized in the observation group. In the control group, we obtained blood vessels from patients who succumbed to acute craniocerebral trauma but had been in good health after active treatment. Part of the artery of family members who volunteered to donate organs was retained as healthy control, though the organ was donated. Detailed clinical information of all participants is outlined in Table 5. Approval for this study was issued by the Ethics Committee of the hospital. All patients or their family members signed informed consent.

Tissue sample acquisition

The anterior edge of the sternocleidomastoid muscle was cut longitudinally under general anesthesia. The skin, subcutaneous, and platysma muscles were then cut layer by layer to expose and open the carotid sheath. The common carotid artery, external cervical artery, and superior thyroid artery were then blocked. Blood vessels were cut longitudinally or obliquely along the beginning of the internal carotid artery after blocking the internal carotid artery under the microscope. Intima and plaques were separated from the bifurcation to the distal end of the internal carotid artery with nerve stripping ions. Plaques of the common carotid artery and external carotid artery were stripped down in turn. The arterial lumen was washed with heparin saline repeatedly. The internal carotid artery was sutured continuously, then opened to exhaust thoroughly before knotting. After inhibiting gas occlusion, it was blocked. Blood flow of external carotid artery, superior thyroid artery, common carotid artery, and internal carotid artery were restored sequentially. We further verified that there is no blood osmosis in the anastomosis and then closed the incision layer by layer. After removing the intima and plaque of the carotid artery, plaque specimens were placed in an EP tube and frozen in liquid nitrogen. In vitro cryopreservation of specimen was conducted in 5 minutes.

RNA extraction

Total RNA was extracted from samples using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's protocol. RNA concentration and purity were assessed using a Nanodrop ND-2000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA). RNA integrity was evaluated using 2% agarose gel. Agilent 2100 bioanalyzer (Agilent, Palo Alto, CA, USA) was used to obtain RNA integrity values (RIN).

Transcriptome deep sequencing

After total RNA treatment with DNase I, mRNA was enriched with Oligo D (T) magnetic beads (New England Biolabs, Inc., Ipswich, MA, USA), and digested into short fragments. mRNA fragment was used as a template, whereas random oligonucleotide was used as a primer. The first strand of cDNA was synthesized using the M-MuLV reverse transcriptase system. Then, the RNA strand was degraded using RNaseH. The second strand was synthesized from dNTPs under the DNA polymerase I system. The double-stranded cDNA was purified, end-repaired, and “a-tailed” for adapter ligation. AMPure XP beads were used to screen out 250-300bp of cDNA for PCR amplification, followed by purification of the PCR products to obtain the library. Agilent 2100 Bioanalyzer (Agilent Technologies, Inc.) and ABI StepOnePlus Real-time PCR System (Applied Biosystems; Thermo Fisher Scientific, Inc.) were applied to quantify and identify the sample library for quality control. Sequencing of mRNA libraries was achieved with Illumina HiSeq TM 2000 (Illumina, Inc., San Diego, CA, USA).

Data quality control

Sequences of the raw data containing a small number of reads with sequencing adaptors or low-quality sequences were filtered out. Since the sequencing process is prone to machine error, we established the quality of sequencing data based on sequencing error rate distribution examination. GC content distribution was also determined. HISAT2 [33] was employed to map clean reads to the reference genome (Ensemble GRCh38.p7).

Differential gene screening

We evaluated the distribution of gene expression levels in different samples. A correlation diagram of gene expression levels between samples was applied to verify the reliability of the experiment and the accuracy of sample selection. Inter-group differences and intra-group sample duplications were evaluated via Principal component analysis (PCA). DESeq2 package in R software (R version 3.6.1: https://www.r-project.org/) was adopted to explore the differential gene expression [34]. Original read count was standardized (normalization), then the BH method was used for multiple hypotheses testing corrections. Differentially expressed genes (DEGs) were screened using the conditions, |logFC|>3, adj.p.Val <0.05. A heatmap of the selected DEGs was constructed using the pheatmap package [35].

GO and KEGG pathway analysis of differential genes

Biological processes (BP), cellular component (CC), and molecular function (MF) categories were explored through GO analysis. KEGG (http://www.genome.jp/) analysis is a systematic approach to assess gene function for the discovery of bioregulatory pathways. Herein, we first converted the gene ID from the Official Symbol of the differential gene obtained by org.Hs.eg. Db package [36] and clusterProfiler package [37] were adopted for GO and KEGG pathway analysis, respectively. Besides, the GOplot package [38] was used for cluster analysis. Adj. p. Val <0.05 denoted statistical significance.

MCODE analysis

Differentially expressed genes (DEGs) were analyzed using Metascape online database (http://metascape.org/gp/index.html#/main/step1) [39]. KEGG pathway enrichment analysis was performed for each MCODE. P < 0.01 denoted statistical significance.

PPI network construction and hub gene screening

STRING [40] database (https://string-db.org/) was used to construct a PPI network. An interaction with a composite score greater than 0.9 was considered statistically significant. Cytoscape [41] is an open-source bioinformatics software program for exploring molecular interaction networks. CytoHubba [42], a Cytoscape Plugin, was used to evaluate PPI network hub genes. The top 20 hub genes were screened out through CytoHubba Plugin by degree. Red colors denoted that the gene had a higher degree.

Bioinformatics analysis of hub gene

The top 20 hub Genes were analyzed via Metascape online database (http://metascape.org/gp/index.html#/main/step1). Gene list enrichments were identified in the following ontology categories: TRRUST, DisGeNET, and PaGenBase. All genes in the genome were used as the enrichment background. Terms with a p-value < 0.01, a minimum count of 3, and an enrichment factor > 1.5 were screened out. Interaction of small-molecule drugs with hub genes was predicted using CMAP online database (https://portals.broadinstitute.org/cmap).