T cell receptor signaling pathway and cytokine-cytokine receptor interaction affect the rehabilitation process after respiratory syncytial virus infection

Data source and data preprocessing

The E-MTAB-5195 dataset was acquired from EBI ArrayExpress database (https://www.ebi.ac.uk/arrayexpress/) (Brazma et al., 2003), which was detected on the platform of Affymetrix GeneChip Human Genome U133 Plus 2.0. E-MTAB-5195 contained 60 blood samples from infants, including 39 samples in the acute phase of infection (14 females and 25 males; mean age = 146.49 days) and 21 samples in the recovery period (11 females and 10 males; mean age = 137.95 days). This study analyzed the expression profile dataset obtained from public database, and no animal or human experiments were involved. Therefore, no ethical review or informed consents were needed.

The raw expression profile data in .CEL format were conducted with format transformation, filling of missing data (median method), background correction (MicroArray Suite method), and data standardization (quantile method) using the R package oligo (version 1.41, http://www.bioconductor.org/packages/release/bioc/html/oligo.html) (Parrish & Spenceriii, 2004).

Differential expression analysis

The platform information of E-MTAB-5195 was downloaded. Combined with the Transcript ID, RefSeq ID, and location information provided in platform information, the mRNAs and lncRNAs in the expression profile E-MTAB-5195 were annotated based on the 19,198 protein coding genes and 3909 lncRNAs included in HUGO Gene Nomenclature Committee (HGNC, http://www.genenames.org/) database (Bruford et al., 2008).

According to disease stage information, the samples were divided into samples in the acute phase of infection and samples in the recovery period. Using the R package limma (version 3.34.0, https://bioconductor.org/packages/release/bioc/html/limma.html) (Ritchie et al., 2015), differential expression analysis for the two groups were carried out. The false discovery rate (FDR) < 0.05 and —log fold change (FC)— > 0.263 were defined as the criteria for screening differentially expressed RNAs (DE-RNAs). Based on the expression levels of the screened RNAs, the R package pheatmap (version 1.0.8, Kolde, 2015) (Wang et al., 2014) was used to perform hierarchical clustering analysis (Szekely & Rizzo, 2005) and present the result in clustering heatmap.

Weighed gene co-expression network analysis (WGCNA) and enrichment analysis

The systems biology method WGCNA can be applied for integrating gene expression and identify disease-associated modules (Li et al., 2018). Using the R package WGCNA (version 1.61, Langfelder & Horvath, 2008) (Langfelder & Horvath, 2008), the DE-RNAs were analyzed to identify the modules significantly related to disease states and clinical factors. The analysis procedures included assumption of the scale-free network, definition of co-expression correlation matrix, definition of adjacency function, calculation of dissimilarity coefficients, and identification of disease-associated modules. The modules having significant correlation with disease symptoms (p-value < 0.05) and higher correlation coefficients compared with the control (grey module) were taken as the significant modules.

Then, DAVID online tool (version 6.8, https://david.ncifcrf.gov/) (Lempicki, 2008) was utilized to perform Gene Ontology (GO) (The Gene Ontology Consortium, 2015) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2016) enrichment analyses for the mRNAs in the significant modules. The p-value < 0.05 was selected as the threshold of enrichment significance.

Co-expression network analysis

For the lncRNAs and mRNAs in the significant modules, pearson correlation coefficients (PCCs) (Hauke & Kossowski, 2011) were calculated using the cor function in R (https://www.rdocumentation.org/packages/stats/versions/3.6.0/topics/cor). Afterwards, visualization of co-expression network was conducted using Cytoscape software (http://www.cytoscape.org/) (Kohl, Wiese & Warscheid, 2011). Moreover, KEGG pathways were enriched for the mRNAs in the co-expression network using DAVID tool (Lempicki, 2008).

Network analysis for the RSV-correlated pathways

Using “Respiratory Syncytial Virus” as keyword, the KEGG pathways directly related to RSV were searched form Comparative Toxicogenomics Database (2017 update, http://ctd.mdibl.org/) (Davis et al., 2013). By comparing the searched pathways and the pathways enriched for the mRNAs in the co-expression network, the overlapped pathways and the mRNAs involved in them were obtained. Furthermore, the clusters (including both lncRNAs and mRNAs) of the RNAs involved in the RSV-correlated pathways were extracted from the co-expression network to construct the network of RSV-correlated pathways.

Differential expression analysis

The raw data of expression profile dataset E-MTAB-5195 were preprocessed (Fig. 1). After annotation was performed based on platform information, a total of 16,984 protein coding RNAs and 939 lncRNAs were obtained. Then, the samples were grouped according to disease stage information. There were a total of 2,489 DE-RNAs between the two groups, including 2,386 DE-mRNAs (1,393 up-regulated and 993 down-regulated) and 103 DE-lncRNAs (62 up-regulated and 41 down-regulated) (Fig. 2). The clustering heatmap of the DE-RNAs indicated that the DE-RNAs could separate the two groups of samples and thus had sample characteristics (Fig. 3).

WGCNA and enrichment analysis

The 2,489 DE-RNAs were further analyzed and screened using WGCNA. In order to satisfy the premise of scale-free network distribution as far as possible, the value of the weight parameter “power” of adjacency matrix needed to be explored. The selection range of network construction parameters was set, and the scale-free topological matrix was calculated. The value of “power” when the square value of the correlation coefficient reached 0.9 for the first time was selected, namely, power = 18. At this time, the average node degree of the constructed co-expression network was 1, which fully conformed to the properties of small world network (Fig. 4A). Then, the dissimilarity coefficients between the gene points were calculated, and the system cluster tree was obtained. The minimum number of genes was set as 70 for each module, and the pruning height was set as cutHeight = 0.99. A total of 13 modules (except of grey module) (Fig. 4B), and the number and type of RNAs involved in each module were listed in Table 1. Finally, the correlation between each module and different disease symptoms were calculated. The black, salmon, blue, tan and turquoise modules were significantly correlated with stage (acute/recovery), meanwhile, brown, blue, magenta and pink modules were significantly related to disease severity (Mild/Moderate/Severe) (Fig. 4C). The RNAs in black, salmon, blue, tan and turquoise modules were taken as RNA set1 correlated with stage. Similarly, the RNAs in brown, blue, magenta and pink modules were considered as RNA set2 associated with disease severity.

The mRNAs in RNA set1 and RNA set2 separately were conducted with enrichment analysis. For the mRNAs in RNA set1, 17 GO terms (such as regulation of transcription and lymphocyte activation) and seven KEGG pathways (such as cell adhesion molecules (CAMs) and primary immunodeficiency) were enriched. For the mRNAs in RNA set2, 20 GO terms (including defense response and response to wounding) and eight KEGG pathways (including cytokine-cytokine receptor interaction and lysosome) were predicted (Table 2).

Co-expression network analysis

For the RNA set1 and RNA set2, the PCCs between the lncRNAs and mRNAs separately were calculated and the lncRNA-mRNA pairs with PCC >0.6 were remained. Subsequently, the co-expression networks for RNA set1 (Fig. 5A) and RNA set2 (Fig. 5B) separately were constructed. In the co-expression network for RNA set1, there were 813 nodes (including 777 mRNAs and 36 lncRNAs) and 3773 edges. In the co-expression network for RNA set2, there were 529 nodes (including 521 mRNAs and eight lncRNAs) and 967 edges. Moreover, 10 (such as hematopoietic cell lineage and primary immunodeficiency) and 11 (such as glutathione metabolism and cytokine-cytokine receptor interaction) KEGG pathways separately were enriched for the mRNAs in the co-expression networks for RNA set1 and RNA set2 (Table 3).

Network analysis for the RSV-correlated pathways

Based on Comparative Toxicogenomics Database, 44 KEGG pathways directly related to RSV were obtained. After comparing the searched pathways and the pathways enriched for the mRNAs in the co-expression networks, four (including hsa04660:T cell receptor signaling pathway, involving CD40 ligand (CD40LG) and RAS guanyl releasing protein 1 (RASGRP1); hsa04514:Cell adhesion molecules (CAMs); hsa04650:Natural killer cell mediated cytotoxicity; and hsa04350:TGF-beta signaling pathway) and two (including hsa04060:Cytokine-cytokine receptor interaction, involving interleukin 1 beta (IL1B), interleukin 1 receptor type II (IL1R2), interleukin 18 (IL18) and interleukin 18 receptor 1 (IL18R1); and hsa04621:NOD-like receptor signaling pathway) overlapped pathways separately were obtained for the mRNAs in the co-expression networks for RNA set1 and RNA set2.

Furthermore, the networks of RSV-correlated pathways were built (Fig. 6). In the pathway network for RNA set1, long intergenic non-protein coding RNA 891 (LINC00891), long intergenic non-protein coding RNA 526 (LINC00526) and long intergenic non-protein coding RNA 1215 (LINC01215) were co-expressed with CD40LG, and LINC01215 was co-expressed with RASGRP1. Moreover, long intergenic non-protein coding RNA 1503 (LINC01503) and SMIM25 co-expressed with IL1B, colorectal neoplasia differentially expressed (CRNDE) co-expressed with IL1R2, BAIAP2 antisense RNA 1 (BAIAP2-AS1) and CRNDE co-expressed with IL18, and SMIM2 5 co-expressed with IL18R1 were involved in the pathway network for RNA set2.