Basic information of study groups
We collected two study groups (propofol group and sevoflurane group) demographic data, including age, gender, body mass index, ASA Physical Status, smoke history, diabetes type Ⅱ, arterial hypertension, coronary artery disease, UICC tumor stage, duration of surgery and duration of one-lung ventilation (Table 1). All of these data were not statistically significant between groups (P value > 0.05). In order to explore the effects of propofol (P) and sevoflurane (S) during lung cancer resection, we further matched two groups and used next-generation sequencing (NGS) to reveal the molecular mechanism (Fig. 1).
Quality of NGS and mapping rate
We sequenced a total number of 16 samples, 4 patients in each anesthesia group (P or S group) with tumor and paired normal samples. However, 2 patient data (one from P group, one from S group) did not pass the quality control step, which resulted in 12 samples (sequence libraries) left. After sequencing, we obtained a total number of 170.05 GB clean reads with an average number of 14.17 GB data for each library. The sequence quality score Q30 ranged from 92.58–93.92% and the average GC content was 47.24% (Table 2). Then we mapped clean reads to human reference genome. The total mapping rate ranged from 96.87 to 97.59% and the uniquely mapped rate was 89.85–92.56% (Table 2). After mapping step, we obtained the expression level for each gene, the distributions of Fragments per Kilobase per Million (FPKM) in all 12 samples were consistent, indicating a good sequencing quality (Fig. 2).
Differentially expressed genes (DEG) identified by single factor model
We first did differential gene analysis with paired comparisons, (1) P group: compare P group cancer samples vs P group normal samples, aiming to discover the effect of propofol on lung cancer; (2) S group: compare S group cancer samples vs S group normal samples, aiming to discover the effect of sevoflurane on lung cancer. Thus, we obtained two DEG lists. The threshold was defined as FDR < = 0.01 and log2FC >=±2. In total, we obtained 810 DEG in P group and 508 DEG in S group.
In order to provide a straightforward and detailed perspective on DEG expression, we used hierarchical cluster method to cluster DEG. From the heatmap, there were clearly two clusters in each group (Fig. 3A, 3B); all cancer samples were in one cluster and normal samples in another cluster, which suggested propofol and sevoflurane had effects during lung cancer resection.
We further compared these two DEG lists and found that 191 DEG were identified in both groups. Then we are wondering whether these 191 DEG can reveal the effects of the propofol and sevoflurane were different on lung cancer or not. However, from statistical perspective, simply look at the overlap of DEG in two comparisons, is not sufficient because too many false positives will be produced when we add up data. And if we directly compare P cancer samples vs S cancer samples, it did not a good comparison due to two reasons, 1) the control group was not the same; 2) the heterogeneity of patients. Therefore, an advanced analysis was needed.
DEG identified by multi-factor model
We developed an integrated statistical model to identify DEG, which takes the information of anesthesia drugs (propofol and sevoflurane) and conditions (cancer and normal) into account (Details in Materials and Methods). The threshold was set as FDR < = 0.05 and log2FC >=±1. As a result, 124 DEG were identified, which responded to not only the condition (normal or cancer), but also the effect of propofol and sevoflurane. A number of 53 DEG were down-regulated and 71 DEG were up-regulated.
After obtaining this comprehensive list of DEG, we performed a principal component analysis (PCA), a statistical procedure to transform data from high dimension into low dimension. This type of plot is useful for visualizing the overall effect of experimental conditions. In PCA figure, four groups (P cancer group, P normal group, S cancer group, S normal group) were separated adequately (Fig. 4A). More importantly, there was no overlap between groups. For example, P group and S group separated very well from principal component 1 (PC1), cancer group and normal group separated very well from principal component 2 (PC2). Overall, PC1 can explain 50% variances and PC2 can explain 23% variances.
We also performed a hierarchical cluster analysis, where the data in the same cluster displayed high similarity (Fig. 4B). From the heatmap, the P group and S group were clearly divided into two clusters; what’s more, the cancer and normal conditions were also separated well in each group, indicating the accuracy of the integrated statistical model and the identified DEG. In conclusion, from these results, it suggested the effects of propofol and sevoflurane during lung cancer resection were different.
Gene ontology (GO) and KEGG pathway analysis
GO enrichment analysis and KEGG pathways analysis were performed on 124 DEG. The threshold was taken as P-value < 0.01. The GO analysis consists three parts: biological process, cellular component and molecular function. Terms that were significant in each part and the number of genes in each process were plotted as a bar graph (Fig. 4C). The biological process was the most interesting and informative one, the top10 significant terms in biological process were 1) positive regulation of transcription, 2) inflammatory response, 3) positive regulation of cell proliferation, 4) immune response, 5) cellular response to lipopolysaccharide, 6) calcium ion-regulated exocytosis of neurotransmitter, 7) chemokine-mediated signaling pathway, 8) cellular response to interleukin-1, 9) cellular response to tumor necrosis factor, 10) chemotaxis.
In KEGG pathway analysis, the top 10 significant pathways with their P value and the number of genes involved in each pathway were plotted in the bubble diagram (Fig. 4D), including 1) Cytokine-cytokine receptor interaction, 2) TNF signaling pathway, 3) Malaria, 4) Salmonella infection, 5) Chemokine signaling pathway, 6) Hepatitis B, 7) Jak-STAT signaling pathway, 8) African trypanosomiasis, 9) NOD-like receptor signaling pathway, 10) Pertussis.
From GO and KEGG pathway analyses, we compared and integrated the significant processes, then chose to focus on seven pathways (Table 3) for downstream analysis. Among these pathways (Cellular response to interleukin-1, Chemokine-mediated signaling pathway, Chemokine signaling pathway, Cytokine-cytokine receptor interaction, Inflammatory response, Immune response, TNF signaling pathway), genes enrolled in these pathways were also extracted (Table 3). Several genes were involved in more than one pathway, indicating important roles in network regulation. Therefore, we chose these genes for qRT-PCR validation.
Candidate genes confirmed by qRT-PCR
We chose four genes (CCL20, CXCR1, CXCL8, TNFAIP3) for validation. In order to better reflect the differences between P normal (Pn), P cancer (Pc), S normal (Sn) and S cancer (Sc), we calculated the ratio of normal samples to cancer samples (Pn/Pc and Sn/Sc). If this ratio is bigger than 1, it means this gene was down-regulated in cancer samples. If this ratio is smaller than 1, it means this gene was up-regulated in cancer samples. From qRT-PCR, three genes (CXCR1, CXCL8 and TNFAIP3) were consistent with RNA-seq results and statistically significant between P and S (Fig. 5), whereas CCL20 showed no significance between P and S. Both CXCR1 and TNFAIP3 showed decreased expression in cancer samples (ratio > 1), however, when compared between P and S, CXCR1 depressed more under S, TNFAIP3 depressed more under P, indicating different level of responses to anesthetic drugs. CXCL8 showed increased expression in S cancer samples (ratio < 1), whereas decreased expression in P cancer samples (ratio > 1), which may suggest cancer cells responded differently to P and S during lung cancer resection.