Transcriptome Proling of Peripheral Blood Mononuclear Cells Reveals COVID-19 Patients Are Not Recovered to Normal After Discharge for 5 Months

Background: Highly pathogenic coronavirus disease-2019 (COVID-19) initiated by severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) infection has swiftly expanded throughout the world, and the fatality rate is still expanding due to the second wave in 2020 winters. This ongoing epidemic threatens public health with its new strain that emerged in some countries and might cause devastating deaths. Therefore, the host transcriptomic profile from patients during recovery is important for understanding this disease. MethodsWe performed transcriptome profiling of the RNAs isolated from the peripheral blood mononuclear cells (PBMCs) of recovered COVID-19 patients at hospital discharge of three months and five months respectively.ResultsOur results exposed diverse inflammatory genes and cytokine profiles to infection in recovered patients, and emphasize the highly expressed genes in COVID-19 patients like CCL4, CCL3, CXCL9, CXCL16, IL10, CSF2, VEGFA showed a decreasing trend in recovered patients. Furthermore, the integrated analysis predicted that JUN, CTSL, DDIT4, RRAS, BIRC5, CTSZ, CCNB2, CDK1, OAS1/2, IFIT3, RSAD2, and TP53I3 genes may be valuable for the recovery of COVID-19 patients. ConclusionsOur analysis confirms the presence of some inflammatory genes in recovered patients, suggesting COVID-19 patients did not return to their normal expression even after 5-months of discharge. Identification of transcriptome profiling of recovered patients provides useful information regarding its pathogenesis and might help for the development of better treatment for COVID-19.

Background SARS-CoV-2 is an enveloped, easily transmittable fatal single-stranded RNA virus consisted of 29,903 nucleotides, which caused an outbreak along worldwide since 2019. It belongs to the order Nidovirales, the family Coronaviridae, and subfamily Coronavirinae of viruses that become a major reason for minor to destructive respiratory tract disorders in humans [1][2][3]. The clinical manifestations of the COVID-19 caused by this virus depend on the host immune response and viral pathogenesis lead from asymptomatic to minor infections (fever, sneezing, anosmia, etc.), to severe disorders (ARDS, organ failure, etc.), and leading to death [4][5][6][7][8]. Patients with older age and other comorbidities (diabetes, COPD, hypertension, etc.) and patients in immuno-compromised states were faced complications and at utmost risk with severe COVID-19 [9][10][11]. While worldwide scientists and clinicians have made countless efforts to ght against the disease, several vaccines, anti-viral and immunomodulatory drugs are in different phases of clinical trials [12][13][14]. In November, the pharma giant P zer and BioNTech (got temporary authorization), Moderna, and a Russian institute have all suggested initial evidence that spikebased vaccines can attain >90% e cacy [15,16].
Recently, several new variants of SARS-CoV-2 have been reported by the UK, South Africa, and Nigeria, including a new variant with N501Y mutation on S protein, VUI-202012/01 from the UK, which has about 71% higher transmissibility than the existing strains [17]. The development of the COVID-19 is quite similar to the 1918 "Spain Flu", especially for the second wave with new variants and increasing spread rate globally. In the second wave of Spain u, the mutated variants become more deadly than the original strain and mortality were very high [18,19]. Predictably more and more new variants will emerge with the rapidly expanding epidemic. Thus, a deep understanding of anti-viral innate and adaptive immune responses is a prerequisite for accurate diagnosis and adopting a therapeutic strategy. Dysregulation of the immune system is thought to be correlated with the acuteness of COVID-19 infections, and enhanced expression of proin ammatory chemokines and cytokines is exhibited in mild and severe patients [4,20,21]. Different studies have revealed that PBMCs are indulgent to SARS-CoV disease and help in the replication of the virus, and are used for transcriptomics evaluations, including changes in the expression of genes in response to infection or treatment [22,23].
Though, how signi cant immune cells population change and basic molecular procedures of the in ammatory responses and pathological processes of the COVID-19 have remained unclear. However, to investigate the biological process and host transcriptional uctuations in response to disease and after the infection is necessary, especially with the emergence of the novel variants from different regions and countries.
Here, we performed transcriptome pro ling of PMBCs from recovered COVID-19 individuals, to map the major expression pro le of pro-in ammatory genes and cytokines induced by infection. We observed lower expressions of CCL4, CCL3, CXCL9, CXCL16, IL10, CSF2, and VEGFA. CXCL3 was down-regulated in all recovered patients except X13-X14. Our study represented a high-resolution transcriptome landscape after the COVID-19 infection in different periods. It helps the researchers to get a better understanding of the host immune response, abnormal expression levels, and cytokine storm uctuations throughout the recovery stage, which might be helpful for the production of vaccines and antibodies for therapeutics.

Methodology
Sample Isolation, Sequencing, and Filtering The purple EDTA anticoagulant tubes were used to collect the whole blood from patients recovering from SARS-COV-2 in the First A liated Hospital of USTC, Anhui, China during their checkups after discharge.
Blood samples were collected for each patient three months and ve months after recovery. Ficoll-Hypaque density gradient centrifugation was applied to separate lymphocytes in peripheral blood. In brief, under the action of centrifugation, Ficoll-Hypaque solution with a density of 1.077 can separate PBMCs from Red blood cells, granulocytes, plasma, and other components. The collected PBMCs were washed twice with DMEM medium and counted under a microscope. One million cells from each sample were sent to BGI (Shenzhen, China) for transcriptome sequencing. A large number of sequencing data was generated in a single run and quality control was then performed to determine whether the sequencing data is suitable for subsequent analysis. The original data were puri ed by SOAPnuke-v1.5.2 (https://github.com/BGI-exlab/SOAPnuke) with parameters -q 0.2 -l 15 -n 0.05. It was performed to remove adaptor pollution, low-quality reads and unknown "N" bases greater than 5%. In the end, clean reads were saved in FASTQ format [24].

Expression Quanti cation
Then output was fed to RSEM-v1.2.8 (http://deweylab.biostat.wisc.edu/rsem/rsem-calculateexpression.html) to calculate the gene expression level of each sample with default parameters. After that, it is converted by FPKM or TPM to obtain the standardized gene or transcript expression level. After obtaining the Read Counts, gene or transcript expression differential analysis was performed. Here, gene quanti cation analysis and other analyses based on gene expression like principal component analysis (PCA), sample correlation heat map, expression quanti cation distribution (boxplot, density map, stacked histogram), and differential gene screening were conducted.

Differential Gene Detection
PossionDis method was used for the detection of differential genes with parameters: Fold Change >= 2 and FDR <= 0.001. The sequencing-based differential gene detection method used a rigorous process to screen DEGs between samples [26]. R-software pheat-map (https://cran.rproject.org/web/packages/pheatmap/) was utilized to perform hierarchical-clustering analysis on the union set differential genes with default parameters. Additionally, we showed the differential expression of genes between different comparison groups and the intersection/union via VENN diagrams.

Functional Enrichment and Variation Analysis
Functional enrichment analysis was executed on all the transcripts and DEGs found from this RNA-Seq, which were compared with KEGG (http://www.genome.jp/kegg/) and GO (http://www.geneontology.org/) databases. According to the annotated classi cation of KEGG Pathway and GO analysis, the phyper-function in R was utilized to execute the enrichment study, compute the P-value, and then carry-out FDR-correction. Generally, the function with Qvalue <= 0.05 was considered as signi cant enrichment. The rMATS-v3.2.5 (http://rnaseqmats.sourceforge.net) statistical model was used to quantify the expression of alternative splicing events with parameters: -analysis U -t paired -a 8. Additionally, gene fusion events were detected by aligning the sequences at the paired-end relationship in the genome and transcript via Ericscript-v0.5.5 (http://ericscript.sourceforge.net/) with default parameters.

Cytokine Storm Detection
Cytokines were screened using Perl and Python scripts and performed GO enrichment, KEGG pathway, and disease enrichment analysis, in which Qvalue <= 0.05 was viewed as a signi cant enrichment.

Results
To explore the in uence and mechanism of SARS-CoV-2 disease on the host immune system after the recovery, we utilized RNA-sequence to nd transcriptome modi cations in PMBC samples from recovered individuals in two different periods.  Table 1.
The sequencing dataset passed stringent high-quality ltering and plotted the statistics of reads mapped to the genome and genes. The mapping rate for the mapped reads of each sample with the selected reference genome and genes total ranges from 87.8 % to 88.62 % and 52.05 % to 60.79 %, respectively. The average alignment ratio of the observed sample comparison genome was 88.13%, and the gene set was 57.03%. Whereas the proportion of reads that map to the unique position of the reference genome and genes ranges from 78.2% to 83.11 % and 46.86% to 56.1%, respectively (Table. S1). We conducted a Randomness analysis, most of the transcripts were entirely covered, and reads were evenly scattered in various regions of the transcript. Additionally, the coverage of the transcripts of each sample suggested that most of the transcripts were fully covered, and reads were regularly distributed in each area of the transcript. Finally, the sequencing saturation analysis for each sample suggested that the amount of sequencing data meets the desired requirements ( Fig. S1A/B/C/D/E/F).

Transcriptome Pro ling of PBMCs
To obtain the gene signature, we found a total of 17686 expressed genes and 87324 expressed transcripts in the PBMCs samples of all recovered individuals. Gene expression overview in the PBMCs comparisons group was examined in the PCA (Fig. 1A). To signify the correlation of gene expression between samples, the Pearson correlation coe cients between every two samples were calculated and re ected in the form of a heatmap (Fig. 1B). To get an overview of gene expression quanti cation in each sample, we formed a boxplot from which the dispersed degree of the data distribution can be observed ( Fig. 1C). The density map was constructed to show the trend of gene abundance as the expression level changes (Fig. 1D). Also, the statistical analysis was used to visualize the number of genes in the diverse FPKM range of each sample (Fig. 1E).

Gene and Transcripts Analysis in Comparison Groups
To get a comparison overview of all genes and transcripts in control and treated groups, we performed the intersection/union analysis of each sample and group. In the rst group of a total of 16309 genes, in group-2, 16387, group-3, 16399, group-4, 16320, and group-5, consisted of 16272 total expressed genes. Overall, a little more gene was present in control individuals (17,186)

Differential Expression Pattern in Comparison Groups
To characterize the genomic changes in SARS-CoV-2 recovered individuals, we evaluated transcriptomic data from PBMCs and found a total of 5317 DEGs, of which 3170 genes were up-regulated and 2147 genes were down-regulated ( Fig. 3A/B). We also identi ed 85,207 differentially expressed transcripts (41,865 down-regulated and 43,342 up-regulated) (Fig. S3). The up-regulated and down-regulated genes of each comparison group were listed up in Table S2. An overview of statistical signi cance among the gene expression in each sample group was observed, which assisted rapid visual recognition of genes with large fold changes (Fig. 4A). Cluster the expression of the DEGs gave a general overview of gene expression variation among each group (Fig. 4B). The number of DEGs was illustrated graphically which exhibited overlapped genes/transcripts in the comparison group ( Fig. 3C/D).

Characterization of Longitudinal Analysis
We did a longitudinal analysis (3 rd to 5 th months) of individuals one by one, which clari ed the number of expressed genes in each sample and those genes that were speci c to each patient. Overall, in controlgroup (3 rd -month) 5254 genes (3107-up, 2147-down), and treated-group (5 th -month) 5273 (3170-up, 2103down) genes were differentially expressed. While in the comparison of X11 to X12, we identi ed a total of 885 and 889 DEGs, respectively, from which 9 genes were speci cally down-regulated in X11, while 13 genes were speci cally up-regulated in X12. In X13, DEGs were 15626 from which 10 genes were particularly down-regulated, whereas X14 contained 15855 genes and only 12 were uniquely showed higher expression compared to X13. From the third group, X5 had 794 and X6 had 755 DEGs, where 5 down-regulated genes were speci c to X5 and 11 up-regulated genes, were speci c to X6. The 3 rd month's sample of X7 differentially expressed 1081 and during 5 th months (X8) 1085 different expressions of genes, from which X7 speci c down-regulated genes were 10 uniquely up-regulated genes. DEGs in X9 were 1408 and X10 was 1411, where 10 genes showed down and 13 genes showed a high expression that speci c to X9 and X10, respectively (Table. S3).

Functional Enrichment and Gene Annotation
We then performed the functional enrichment and classi cation analysis to identify the changes in ve different groups' DEGs. GO annotation of DEGs was classi ed and mapped. Among the up-regulated genes of group-1, 536 genes were involved in biological process, 569 in a cellular component, and 537 in molecular function. While from the up-regulated genes of the other groups (2,3,4,5), 569, 525, 624, 681 (biological process), 583,548, 656, 699 (cellular component) and 563,530,618, 677 (molecular function) genes were involved respectively. Based on GO enrichment analysis, up-regulated genes were enriched in a cellular component like "mitochondrial inner membrane, main axon, phagocytic vesicle, plasma membrane, lysosome, cytoskeleton, etc.". While in-case of molecular function genes were signi cantly enriched in "ATPase activity, peroxisome proliferator-activated receptor binding, cytokine/chemokine activity, heme and IgG and nucleotide-binding, etc." According to the Go_P annotation classi cation, they were signi cantly enriched in the biological process including "blood coagulation, oxygen transport, neutrophil chemotaxis, cellular response to DNA damage stimulus, positive regulation of interferon, immune system process, neutrophil degranulation and so on". In group-2, GO annotation analysis revealed that up-regulated genes were annotated and involved in biological processes, cellular components, and molecular functions, but not signi cantly enriched (

Pathway and Disease Enrichment Analysis
KEGG Pathway-based evaluation was performed to further gure out the biological function of DEGs in each group. KEGG Pathway classi cation involved in genes was divided into several branches. A classi ed statistic was further made for the metabolic pathways under each branch. Interestingly, the KEGG enrichment study discovered that up-regulated genes were enriched in many pathways (Leishmaniasis, Williams-Beuren/Stickler syndrome, and PPAR signaling, etc.) and also showed enrichment for several diseases (NK cell defects, Neutropenic disorders, Agammaglobulinemias, TNDM, and polydactyly disorders, etc.) Functional pro ling of the PBMCs up-regulated DEGs were recognized important terms like an in ammatory response, cytokine activity, chemokine activity, and type-1-interferon signaling pathway signifying a cellular-mediated response to the immune system ( Fig.  S5A/B/C/D). While down-regulated genes were also enriched in many pathways (TNF/IL-17 and Chemokine signaling, and involved in malaria and cancer pathways, etc.) and signi cantly enriched in many diseases like Hypertrophic, cardiomyopathy, etc. While group-2 down-regulated genes were not signi cantly enriched in any pathway (Fig. S6A/B/C/D).

Variation Analysis
We also performed statistics of alternative splicing events for each sample and each comparison group (Fig. S7). Gene fusion events were detected by aligning the sequences at the paired-end relationship in the genome and transcript and also displayed their location on chromosomes (Fig. S8).

Regulated Cytokines Analysis
To evaluate the immunological response in patients after recovery, we also fetched genes involved in cytokine storm and identi ed a total of 108 cytokines (59-up and 49-down). While in the case of comparison groups (1, 2, 3, 4, 5), a total of 14, 11, 38, 7, 7 up-regulated and 17, 12, 2, 29,27 downregulated cytokines were observed, respectively (Table. S4). Moreover, JUN, FOS, IRF1, DUSP1, MKI67, IFI6, TNFSF13, ISG15, CSF1, and S100A8 genes were up-regulated in some individuals. CCL4, CCL3, CXCL9, CXCL16, IL10, CSF2, VEGFA, IL10, CXCL3, FOSB, FOSL2, IFNG, and many other genes showed less expression. We also observed some unusual behaviors of most cytokines in all groups, because they show different expression levels in different groups such as IL2R, IL1B, and IL6 were up-regulated in one and were down-regulated in another. KEGG pathway and GO enrichment, and disease enrichment analysis for the regulated cytokines exposed certain important pathways that were involved in severe immune degenerative disorders (Fig. 5A/B/C/D/E).  [27,28]. However, the RNA virus tends to evolve with different variants after generations [29], so, SARS-CoV-2 can generate various mutations, which will increase the transmissibility and might increase mortality rate or affect speci c people like youngers. Therefore, host transcriptomic pro le from different patients is important for accurate diagnosis and therapeutics. Our results exhibited a diverse immune response of recovered COVID-19 patients, as shown by the differential expression pattern of genes and cytokines. We identi ed gene signatures from each recovered individual and each comparison group and found 5317 genes were differentially expressed (3170-up and 2147-down). In the case of each comparison group, we found 898, 1143, 3760, 41095, and 1421 DEGs in group-1,2,3,4, and 5 respectively. The overaggressive response of the immune system leading to immunopathology was observed in previous studies [30]. We also did a We also observed in ammatory cytokines and chemokines, like DUSP1/6, IRF1, JUN, IFITM3, and FOS and interferon-stimulated genes (IFRD1, IRF6, and IFI1), were expressed at high levels in recovered patients, which were consistent with the previous study, indicating that these genes were associated with CD14++ in ammatory monocytes and CD4+ T cells. The high expression of some genes IL1B, IL6, CSF1, and CSF2 may be involved with cytokine storm in COVID-19 patients [31]. IL6 was highly expressed in acute COVID-19 patients compared to healthy people, and no signi cant difference was observed between the acute and recovering patients [32]. In our analysis, IL6 expression was down in recovered patients except for X13-X14, and no expression was observed in X5-X6. In contrast, we observed low expression of JUNB, CCL3, CCL4, CSF2, and no differential expression for CXCR4, and KLF6, which showed that reduced expression of in ammatory genes after the recovery. The dysregulation in the monocyte population balance in patients showed that classical-monocytes CD14++ enhanced circulation to stimulate in ammation in SARS-CoV-2 infection [31]. In this study of recovered patients, we demonstrated that IFITM3 expression was elevated in the X5-X6 group, compared to X13-X14, and X9-X10. IFITM3 was involved in anti-viral response and provides immunity against infection by interrupting the intracellular cholesterol homeostasis and hinders the virus entry into the cytoplasm [33,34]. Former study at the whole transcriptome level showed some sets of genes that were upregulated during COVID-19 recovery also showed high expression in some patients of our study even after 5-months of recovery like OSM, IL1B, JUN, NR4A3, AREG, DUSP8, and CD69, whereas NR4A3, ZEB1, IL1B, OSM, MAP3K8, SOSC3, DUSP8 genes presented decreasing trend in some patients after 3 and 5 months of recovery. It disclosed that up-regulation of JUN, AREG, and CD69, and down-regulation of ZEB, MAP3K8, and SOSC3 during the recovery stage exhibited the same trend even after 5-months of recovery. Additionally, according to the previous investigation, some antiviral genes like OAS1, IFIT3, RSAD2 were downregulated at the rehabilitation stage, but current analysis showed these genes were up-regulated after recovery, except group-2 and group-4 where IFIT3 exhibited decreasing trend even after 5-months of recovery [35].
Single-cell RNA-seq exhibited up-regulation of the TNF/IL1B-driven in ammatory response in COVID-19 [36]. Furthermore, previously stated that in COVID-19 highly expressed IL1B gene may be a novel candidate target, however, its expression was decreased in observed recovered patients, except in the X5-X6 group where it had a high expression. There was no differential expression of IL1B in the X13-X14 group, might be its expression uniquely depends on other factors like age, because the X13-X14 individual was 8 years old and X5-X6 was 58 years old. IL1B and CXCL8 were increased during aging and upregulate in COVID-19 patients [37] still showed high expression in aged people after recovery of this infection. Many genes behave differently in both of these groups as compared to others. We observed high expression of CCL2, CXCL10, S100A8, and other B cell activation-related genes TNFSF13, and TNFSF13B in one group (X5-X6), but lessen the expression of TNFSF12-TNFSF13 and TNFSF13B in another group (X13-X14). Based on the previous study, IL18, TNFSF13, TNFSF13B, IL2, and IL4 to elevate the cells proliferation and then produced antibodies into the blood may be helpful for the recovery [31] but we did not observe differential expression of IL2, IL4, and IL18 in our recovered patients. While CCL13, CXCL8 (related with chemoattraction of macrophages or neutrophils (innate-immunity)) [32], CXCL2, and CXCL10, were elevated in COVID-19 patients compared to recovering stage, exhibited altered expression in recovered patients like to up-regulate in X13-X14 and downregulate in others, while no expression in X5-X6. Might be CCL2, CCL13, CXCL2, and CXCL10 were also related to aging. While CCL13 was normal in recovered patients. Meanwhile, a previous study showed that OAS2 and IL16 linked with T cells (adaptiveimmunity) were highly expressed in COVID-19 patients with recovery stage and healthy people. In this current analysis, we also observed high expression of OAS2 in some groups of recovered patients even after 5 months of recovery. Some analysis described that the enhanced expression of CXCL10, TNFSF14, S100A8, IL6, and OSM genes in COVID-19 infection is strongly associated with clinical severity [38,39]. We observed down the expression of TNFSF14 in group-4 and group-5 of recovered patients. Many studies had stated notably higher levels of in ammatory cytokines which were linked with an acuteness of disease in MERS, SARS, and SARS-CoV-2 patients suggested that in ammation was an important component of the immune response during COVID-19 infection [40,41].
Different in ammatory cells and monocyte populations might play a key role as they were recognized to fuel in ammation [42][43][44]. We also observed less expression of IL10, IL1A, IL1RN, IL4I1, IL21R, IL31RA, ILDR2, and IL2RA in recovered COVID-19 individuals. Several studies demonstrated that upregulation of IL10 and IL2RA was involved in anti-in ammatory signaling during the infection [45,46]. The increased level of CXCL10 was correlated with disease severity [47,48], thus most of the recovered patients showed decreased expression of CXCL10, CCL2, CCL3, and CCL4. This altered expression proposed that these cytokines might be a helpful candidate target for COVID-19 screening and therapeutics. A previous study showed a correlation between COVID-19 pathogenesis and excessive release of these cytokines [23].
We found some decreasing trends of genes such as CTSL, DDIT4, RRAS, BIRC5, CTSZ, CCNB2, CDK1, and TP53I3 in some recovered patients, which were enriched to P53 signaling pathways and apoptosis according to the previous study, while TNFSF10 showed an increasing trend in recovered patients. The former analysis presented those lymphocytes had been signi cantly reduced in COVID-19 patients. While, TP53, was a signi cant gene in the procedure of apoptosis, showed no altered expression in recovered patients [23,49]. Moreover, we identi ed decreased expression of VEGFA genes (growth factor family) in few patients. FCGBP encodes IgG-binding protein's Fc fragment, which has been stated with the virus disease and viral vector design. It was highly expressed in some recovered patients (X9, X10, X13, X14) in this study, while, FCGBP was also showed high expression in infected COVID-19 patients [50]. Overall, signi cantly expressed genes that were enriched in COVID-19 patients (CCL4, CCL3, CXCL9, CXCL16, IL10, CSF2, VEGFA) showed a decreasing trend in recovered patients and these genes might be the hallmark of COVID-19. The presence of in ammatory gene signatures even after 5-months of discharge, suggesting these COVID-19 patients did not return to their normal life yet, was consistent with a previous study of single-cell sequencing of COVID-19 patients in the recovery phase [31]. Furthermore, we performed gene annotation, functional enrichment, and classi cation analysis in PBMCs to identify the involved disease pathways and biological function of DEGs. The down-regulated cytokines were mostly involved in in ammatory bowel disease, T B+severe combined immunode ciency syndrome, allergic rhinitis and type 1 diabetes mellitus, and so on. Additionally, we also carried out alternative splicing and gene fusion analysis, which might have a pivotal role in SARS-COV-2 illness precision medicine.

Conclusions
In conclusion, by studying the differential expression of genes after SARS-CoV-2 infection, we provided deep perceptions on the possible usage of transcriptomics data and pointed out the prominent in ammatory genes and cytokines in recovered patients, which might explain immune response alterations after infection, and why some patients feel unwell after being discharged. Moreover, it is bene cial to gure out the disease consequences. Additionally, the cytokine expression pro le suggested altered expression of pro-in ammatory cytokine might be a hallmark of recovered patients. Furthermore, a deep investigation is required to verify either these genes could be used as a potential biomarker or have potential medical e cacy in immunotherapies of COVID-19.

Limitations
In this article the limited number of patient samples assessed and uninfected controls are also not included not only due to time limits, but also limited patients. MA and SD conceived the presented idea, prepared the samples, analyzed the data, wrote the original draft, and formatted the manuscript for submission. DZ, ZK, HH and HM participated the research, reviewed and edited the original version of the manuscript. TJ conceptualized the main idea, provided funding during the whole study, and supervised the whole paper. All the authors read and approved the nal version of the manuscript for publication.  Genes-Transcripts Overview. It is exhibited that number of genes (A, B) and transcripts (C, D) intersection and union in each 3rd to 5th-months sample. Each circle represented group of gene sets, and superimposed areas by circles represented the intersection. Non-overlapping part indicated the uniquely expressed genes, and the numbers on the gure indicated the number of genes in the corresponding area.   Pathway Enrichment of Cytokines. In Go-CPF analysis the X-axis represented the number of genes annotated to the GO entry, and the Y-axis showed the GO functional classi cation. While in pathway enrichment X-axis is the enrichment ratio (the ratio of the number of genes annotated to an entry in the selected gene set to the total number of genes annotated to the entry in the species), Y-axis is KEGG