Assessment of the targeted effect of Sijunzi decoction on the colorectal cancer microenvironment via the ESTIMATE algorithm.

OBJECTIVE: Sijunzi decoction (SJZD) was used to treat patients with colorectal cancer (CRC) as an adjuvant method. The aim of the study was to investigate the therapeutic targets and pathways of SJZD towards the tumor microenvironment of CRC via network pharmacology and the ESTIMATE algorithm.
METHODS: The ESTIMATE algorithm was used to calculate immune and stromal scores to predict the level of infiltrating immune and stromal cells. The active targets of SJZD were searched in the Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (TCMSP) and UniProt database. The core targets were obtained by matching the differentially expressed genes in CRC tissues and the targets of SJZD. Then, GO, KEGG and validation in TCGA were carried out.
RESULTS: According to the ESTIMATE algorithm and survival analysis, the median survival time of the low stromal score group was significantly higher than that of the high stromal score group (P = 0.018), while the patients showed no significant difference of OS between different immune groups (P = 0.19). A total of 929 genes were upregulated and 115 genes were downregulated between the stromal score groups (|logFC| > 2, adjusted P < 0.05); 357 genes were upregulated and 472 genes were downregulated between the immune score groups. The component-target network included 139 active components and 52 related targets. The core targets were HSPB1, SPP1, IGFBP3, and TGFB1, which were significantly associated with poor prognosis in TCGA validation. GO terms included the response to hypoxia, the extracellular space, protein binding and the TNF signaling pathway. Immunoreaction was the main enriched pathway identified by KEGG analysis.
CONCLUSION: The core genes (HSPB1, SPP1, IGFBP3 and TGFB1) affected CRC development and prognosis by regulating hypoxia, protein binding and epithelial-mesenchymal transition in the extracellular matrix.

Introduction

Cancer development and progression are complex processes caused by the accumulation of genetic modifications in cancer cells and influenced by the surrounding microenvironment. Cancer cells recruit vascular and stromal components (including immune cells, fibroblasts, cytokines, and the extracellular matrix that surrounds them) to the tumor microenvironment (TME), and the activated TME in turn modifies the malignant behaviors of cancer cells [1]. Numerous studies have indicated that infiltrating immune cells in the TME fail to execute antitumor functions but interact intimately with tumor cells to promote oncogenesis and progression. Both innate and adaptive immune cells are involved in this cross-talk [2].

Colorectal cancer (CRC) is a kind of digestive tract tumor, mostly affecting the colon, that is influenced by multiple factors, such as genetics, environment, and diet. With changes in dietary structure and the deterioration of the natural environment, the incidence and mortality of CRC continue to rise [3-5]. In 2018, there were 18.1 million new cancer cases and 9.6 million cancer-related deaths worldwide. China has the highest cancer incidence and mortality rates in the world, and CRC ranks fourth in the incidence and second in the mortality rate among malignant tumors in China [6]. At present, surgical resection is still the only treatment that can completely cure CRC. However, since the initial symptoms of patients with CRC are not obvious, most patients are diagnosed at the late stage when surgery is not an option [7]. Moreover, up to 15% of patients develop major complications, and postoperative recurrence and metastasis are common [8]. Although the chemotherapy combination of 5-fluorouracil with oxaliplatin has been shown to reduce the rates of metastasis and recurrence in advanced CRC [9, 10], long-term chemoradiotherapy, like many chemotherapy-related strategies, has adverse effects, such as organ damage, myelosuppression, cumulative neuropathy, severe reduction in the quality of life of patients and drug resistance [11, 12].

Traditional Chinese medicine (TCM) plays an indispensable role in the integrative treatment of advanced CRC and has been widely used in China and other parts of Asia. It has been widely reported that TCM could prolong the survival time and improve the quality of life of patients with CRC by preventing tumorigenesis, suppressing tumor growth, and reducing metastasis and recurrence [13-16]. Moreover, as an adjunct therapy, TCM could increase the sensitivity and alleviate the side effects of chemotherapy or radiotherapy, improve immunity, lessen the damage induced by conventional treatments, and ameliorate bone marrow suppression [11, 17–20]. Sijunzi decoction (SJZD), which consists of Panax Ginseng C. A. Mey, Atractylodes Macrocephala Koidz, Poria Cocos Wolf and licorice, is a classical prescription for curing temper deficiency syndrome in TCM. It has long been widely used in the clinic to treat CRC since it combats vomiting, diarrhea and nausea to restore homeostasis in the digestive tract [21-25]. Emerging evidence has shown that SJZD is beneficial for advanced CRC treatment, as it improves the quality of life of patients with Qi deficiency syndrome [26]. However, previous studies mainly focused on the differentially expressed genes (DEGs) between tumor tissues and nontumor tissues, but these studies were often insufficient and limited because they ignored the influence of TCM on the TME. Here, we predicted the levels of infiltrating stromal and immune cells by algorithmically calculating the stromal and immune score. The therapeutic targets of SJZD and the DEGs between different infiltration levels were intersected, and the core genes affected the TME of CRC were found. The result will provide the TME of CRC with a novel auxiliary therapeutic perspective.

Materials and methods

Data mining and processing of the CRC expression profile chip

The Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database was generated and maintained through the National Center for Biotechnology Information (NCBI) Gene Expression database. The GEO includes high-throughput gene expression data submitted by research institutions around the world [27]. The datasets were considered to be appropriate for further analysis according to two criteria: (1) studies with CRC tissue samples and (2) studies with clinical prognostic data available for the samples. Based on these criteria, GSE31595 [28] and GSE17536 [29] were evaluated.

Potential variables and heterogeneity are commonly considered to be major sources of variability and bias in high-throughput experiments [30]. The obtained sample data were prone to error due to the small sample size and differences in sample handling and study design. Therefore, all samples from both datasets were integrated to significantly increase the sample size. Then, the scale function was used in R 4.0.2 software to normalize and standardize the samples.

Evaluation of the infiltration degree of stromal cells and immune cells based on the Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) algorithm

The ESTIMATE algorithm [31], which utilizes the unique properties of the transcriptional profiles of cancer samples to infer tumor cellularity and different degrees of infiltration of tumor and normal cells, was applied. To identify specific signatures related to the infiltration of stromal and immune cells, immune and stromal scores were calculated with the ESTIMATE algorithm, and these scores reflected the predicted level of infiltrating stromal and immune cells. In this study, the ESTIMATE package (https://CRAN.R-project.org/package=EstimateGroupNetwork) in R 4.0.2 was used to perform ESTIMATE scoring of the CRC samples. Besides, in order to compare the variation trends of immune and stromal score at different pathological stages, box plot was used to visualized the distribution characteristics of immune and stromal score at different American Joint Committee on Cancer (AJCC) stages.

Prognostic analysis of the ESTIMATE score

Based on the results of the ESTIMATE score (stromal score or immune score), the 214 samples were divided into two groups (high: ESTIMATE score ≥ median, and low: ESTIMATE score < median). To determine the validity of the ESTIMATE score, a prognosis analysis was utilized to analyze the relationship between the ESTIMATE score and the survival time of the patients with CRC, and survival curves were evaluated by Kaplan–Meier method. In the survival analysis, the clinical data of all patients were obtained from the previously downloaded data sets. Overall survival (OS) was defined as the time from the beginning of random assignment to death (for patients still alive at the end of the study, the time of death was the end of follow-up).

Gene set enrichment analysis

To explore the significant biological pathway of different infiltration levels in the TME of CRC, gene set enrichment analysis (GSEA, http://software.broadinstitute.org/gsea/index.jsp) was performed in the immune and stromal group respectively. The expression matrix and grouping information (stromal high vs stromal low, immune high vs immune low) of 214 samples were used as input data. GSEA v4.1.0 software was used to evaluate the enrichment of biological pathway with “weighted” enrichment statistics. Enriched gene sets with false discovery rate (FDR) < 0.25, |normalized enrichment score (NES)| > 1, and nominal P < 0.05 were regarded as being statistically significant [32].

Collection of active ingredients

The Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (TCMSP, http://tcmspw.com/tcmsp.php/) [33] was used to retrieve active compounds of SJZD based on absorption, distribution, metabolism, and excretion (ADME) [34]. This is a unique pharmacology platform devised for TCM to capture the disease-herb-target relationship.

As oral administration was the chief route of drug delivery of TCM, oral bioavailability (OB) was set as the major parameter in active constituent retrieval [35]. In this study, OB≥30% and drug likeness (DL)≥0.18 were used as the standards to screen the effective active compound of SJZD [36].

Drug target screening and core target gene acquisition of SJZD

The effective active ingredients acquired from TCMSP in the previous step were further searched to obtain the corresponding target proteins. The standardized genes were retrieved through the UniProt database (https://www.Uniprot.org/) [37]. The active genes of Panax Ginseng C. A. Mey, Atractylodes Macrocephala Koidz, Poria Cocos Wolf and licorice were sorted out and introduced into Cytoscape V3.7.2 to establish the "compound-target" network diagram to visualize the overall situation of drug compounds and targets. The network topology analysis plug-in CytoNCA [38] in Cytoscape V3.7.2 was used to obtain the core compound and target gene of SJZD with a degree of freedom > 34 (twice the median). "Degree" reflects the number of links with each node in the network and the interaction relationship between nodes.

Screening of core genes related to SJZD in CRC

The stromal scores and immune scores of the samples were used as the grouping parameters for differential expression analysis. Stromal scores were classified into two groups: high stromal score group (HSS) and low stromal score group (LSS). Similarly, high immune score group (HIS) and low immune score group (LIS) were obtained. The DEGs of stromal groups or immune groups based on the TME were obtained using differential expression analysis (|Log2FC| > 2, adjusted P ≤ 0.05) according to the score groups (HSS vs LSS, HIS vs LIS) in the R package limma [39]. The therapeutic targets of SJZD for CRC obtained previously were intersected with these DEGs to obtain the core target genes of SJZD in CRC.

Analysis of core target genes related to SJZD in CRC

To explore the specific pathways and biological processes of the core target genes involved, gene ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed based on the clusterProfiler package [40]. Protein-protein interaction (PPI) were analyzed and constructed using the Search Tool for the Retrieval of Interacting Genes (STRING) database (https://string-db.org/cgi/input.pl) [41], which could retrieve further protein interactions of the core targets. The PPI network was visualized in Cytoscape V3.7.2. Topological analysis was carried out, and the degree value was used as a reference to determine the important core targets. To investigate the effect of these core target genes on the prognosis of CRC patients, a survival analysis was performed using the clinical data from GEO database (GSE31595 and GSE17536). In addition, considering the limited number of patient samples, we performed further survival analysis validation by using 530 colon adenocarcinoma (COAD) samples from the TCGA database.

Survival curves were drawn to screen core target genes with prognostic value (P ≤ 0.05). The core target genes with prognostic value were called microenvironment-related prognostic genes.

Results

Immune scores and stromal scores were associated with CRC stage

First, the gene expression profiles and clinical characteristics of all 214 patients pathologically diagnosed with CRC were obtained from the GEO database. Among them, 103 patients were female (48.1%), and 111 patients were male (51.9%). The pathological classification of AJCC stage revealed 24 (11.2%) cases of stage I disease, 77 (36.0%) cases of stage II disease, 74 (34.6%) cases of stage III disease and 39 (18.2%) cases of stage IV disease.

According to the ESTIMATE algorithm, the ranges of stromal scores and immune scores were −323.92~5403.97 and -859.48~6701.10, respectively. According to the overall trend of the box plot, the stromal and immune score generally increased with the progression of AJCC stage (Fig 1A and 1B). Just in terms of the median value, the median of stage IV is lower than that of stage III, but still higher than that of stage II and stage I (Fig 1A). Similarly, the rank order of immune score for CRC stages was stage III > stage IV > stage II > stage I (Fig 1B), indicating that immune scores and stromal scores based on the ESTIMATE algorithm were positively correlated with the pathological diagnosis and clinical stage to a certain extent.

Fig 1

Results of correlation analysis.

(A) Association between stromal score and CRC AJCC stage. (B) Association between the immune score and CRC AJCC stage. (C) Correlation between the stromal score group and the survival time of patients with CRC. (D) Correlation between the immune score group and the survival time of patients with CRC.

To assess the potential relationships between overall survival (OS) and immune scores and stromal scores, all the patients were divided into high score groups (≥ median ESTIMATE score) and low score groups (< median ESTIMATE score). Kaplan‐Meier survival curves indicated that the median survival time of the low stromal score group was significantly higher than that of the high stromal score group (Fig 1D, P = 0.018), while the patients showed no significant difference of OS between different immune groups (Fig 1C, P = 0.19)

Results of gene set enrichment analysis (GSEA)

In the C5 enrichment analysis results, 658 gene sets were significantly enriched in the high score group; there were 10 enriched gene sets in the low score group. According to the immunization score-based grouping, 678 and 1 gene sets were significantly enriched for the high group and low group, respectively (S1 Appendix).

The results of C5 enrichment analysis in the high stromal score group showed that the infiltration of stromal cells in the TME of CRC might be related to the following biological functions: cytokine mediated signaling pathway, G protein coupled receptor signaling pathway coupled to cyclic nucleotide second messenger, B cell activation, side of membrane, and adaptive immune response (Fig 2B). According to the GSEA results, cytokine mediated signaling pathway was the most enriched biological function in high stromal score group (Fig 2A, S1 Appendix), which could modulate the ability of cancer cell to grow, proliferate, invade, and metastasize [42].

Fig 2

Estimate score C5 enrichment analysis results based on the ESTIMATE algorithm.

(A) GSEA of top one pathway enriched for high stromal group (NES = 2.54, nominal P < 0.0001). (B) The top5 significant biological pathways in high stromal score group (sort by NES). (C) The enrichment plot of amide binding (NES = 2.61, nominal P < 0.0001). (D)The top5 significant biological pathways in high immune score group (sort by NES).

The results of C5 enrichment analysis in the high immune score group suggested that the infiltration of immune cells in the TME of CRC might be related to the following biological functions: amide binding, response to wounding, endosomal part, calcium ion transmembrane transport, and endocytic vesicle (Fig 2D). The biological function of amide binding was of the highest NES in the high immune score (Fig 2C, S1 Appendix), which indicated that this gene set might regulate the immune response in patients with CRC by inducing related amide binding, that was consistent with the previous research [43].

DEGs

The DEGs of HSS vs LSS and HIS vs LIS were integrated to obtain the union of 1557 DEGs (S2 Appendix). Heat map showed the different expression levels of the top 50 genes between the stromal score groups and immune score groups (Fig 3A and 3B). The results showed that 357 genes were upregulated and 472 genes were downregulated between the immune score groups (|logFC| > 2, adjusted P < 0.05, Fig 3C); 929 genes were upregulated and 115 genes were downregulated between the stromal score groups (|logFC| > 2, adjusted P < 0.05, Fig 3D).

Fig 3

The expression profiles and differential expression analysis results of the top 50 characteristic genes between the high and low stromal and immune score groups.

(A) The expression heat map for the top 50 DEGs based on immune score grouping. (B) The expression heatmap for the top 50 DEGs based on stromal score grouping. (C) Volcano maps based on differential expression analysis of immune score groups. (D) Volcano maps based on differential expression analysis of stromal score groups.

Network pharmacological analysis indicated the importance of each gene in SJZD treatment

Through the TCMSP database, a total of 139 active compounds that met the screening conditions were obtained. Some active ingredients of the four herbs of SJZD were shown in Table 1. Complete information on the active ingredients of SJZD can be found in S1 Table. The official human genes corresponding to the 230 target proteins were retrieved by using UniProtKB. Compound-target networks were constructed by using Cytoscape (Fig 4A). The degree value represents the degree of connection between the targets, and a high degree of connection indicates great importance. The top 12 targets are shown in detail in Table 2. In the network diagram, the square nodes are drug targets, and the red nodes are compounds of SJZD. The larger the node is, the greater the importance of the drug component or target is.

Table 1

Basic information for some active components of SJZD.

Mol ID	Chemical component	OB (%)	DL	Herb
MOL000273	(2R)-2-[(3S,5R,10S,13R,14R,16R,17R)-3,16-dihydroxy-4,4,10,13,14-pentamethyl-2,3,5,6,12,15,16,17-Octahydro-1H-cyclopenta[a]phenanthren-17-yl]-6-methylhept-5-enoic acid	30.93	0.81	Poria Cocos Wolf.
MOL000275	Trametenolic acid	38.71	0.80
MOL000276	7,9(11)-Dehydropachymic acid	35.11	0.81
MOL000279	Cerevisterol	37.96	0.77
MOL001484	Inermine	75.18	0.54	licorice
MOL001792	DFV	32.76	0.18
MOL000211	Mairin	55.38	0.78
MOL002879	Diop	43.59	0.39	Panax Ginseng C. A. Mey.
MOL000449	Stigmasterol	43.83	0.76
MOL000358	β-Sitosterol	36.91	0.75
MOL003648	Inermin	65.83	0.54
MOL000022	14-Acetyl-12-senecioyl-2E,8Z,10E-atractylentriol	63.37	0.3	Atractylodes Macrocephala Koidz.
MOL000033	(3S,8S,9S,10R,13R,14S,17R)-10,13-dimethyl-17-[(2R,5S)-5-propan-2-yloctan-2-yl]-2,3,4,7,8,9,11,12,14,15,16,17-dodecahydro-1H-cyclopenta[a]phenanthren-3-ol	36.23	0.78
MOL000049	3β-Acetoxyatractylone	54.07	0.22
MOL000072	8β-Ethoxy atractylenolide Ⅲ	35.95	0.21

Fig 4

Network pharmacology analysis of SJZD.

(A) Drug-active ingredient-target network diagram. The blue nodes are drug targets, and the red nodes are compounds of SJZD. (B) Intersection of 230 therapeutic targets and 1557 DEGs related to SJZD. (C) PPI network diagram of the 52 core genes. The larger nodes in the inner ring represented the most important hub nodes. The smaller nodes in the outer ring represent the other nodes (several nodes with degrees of 1 were deleted).

Table 2

A partial information table for the core targets.

Targets	Degree	Targets	Degree
PTGS2	86.0	PRSS1	57.0
ESR1	77.0	GSK3B	56.0
NOS2	67.0	SCN5A	53.0
PPARG	67.0	ESR2	53.0
AR	67.0	CCNA2	53.0
HSP90AB1	64.0	PTGS1	52.0

52 core target genes were obtained after the intersection of 1557 DEGs and 230 protein-coding genes (Fig 4B). These 52 genes were utilized to obtain a PPI network diagram (Fig 4C). Genes with larger nodes and deeper colors have higher degree values.

Biological pathways related to the target DEGs

The functional enrichment analysis of the 52 core target genes showed that they were mainly enriched in the signal transduction and expression regulation in TME of CRC. Among the enriched GO terms, there were 153 significantly enriched biological processes (BP), 17 significantly enriched cellular components (CC) and 25 significantly enriched molecular functions (MF) (P < 0.05). The main enriched GO terms included response to hypoxia and inflammatory response (Fig 5A), extracellular space, extracellular region and extracellular exosome (Fig 5B), and protein binding, receptor binding, enzyme binding and serine−type endopeptidase activity (Fig 5C). In addition, all pathways identified as enriched by the KEGG analysis were associated with TNF signaling pathway, transcriptional misregulation in cancer and pathways in cancer (Fig 5D).

Fig 5

GO (BP, MF, CC) enrichment analysis and KEGG pathway enrichment analysis of 52 nodes.

Survival analysis and validation

Survival analysis of the 52 core target genes showed that HSPB1, COL3A1, SPP1, OLR1, IGFBP3 and TGFB1 were associated with adverse prognosis of the 214 GEO CRC patients. Validation of 530 TCGA samples further confirmed that HSPB1, SPP1, IGFBP3 and TGFB1 were indeed associated with poor prognosis. These 4 genes were identified as the main targets of SJZD in the CRC TME, called the core prognostic targets of SJZD for patients with CRC. The survival curves were shown in Figs 6 and 7.

Fig 6

Survival analysis results for the 52 core target genes (performed in 214 CRC samples from GEO).

Fig 7

Verification results of survival analysis (performed in 530 COAD samples from TCGA).

Discussion

Network pharmacology analysis clarified a potential mechanism underlying the role of SJZD in the CRC TME and indicated that HSPB1, SPP1, IGFBP2 and TGFB1 could be potential therapeutic targets in the CRC TME. GO and KEGG analysis of 52 core genes showed that they were mainly enriched in the response to hypoxia, the extracellular space, protein binding and the TNF signaling pathway. After a comprehensive analysis of these results, we supposed that the effect of SJZD on the microenvironment of CRC might be realized via the following three mechanisms: SJZD remodels the TME and prolongs the survival of patients with CRC by altering the expression of SPP1, TGFB1A and IGFBP3, attenuating tumor hypoxia, decreasing the reactive oxygen species (ROS) level and enhancing the sensitivity to drugs. Hypoxia is a hallmark of the TME and plays roles in metastasis, therapy resistance and adverse clinical outcomes [44]. The prognosis of CRC is significantly related to the hypoxia response, immunity, and combined treatment response, and CRC commonly exhibits an immunosuppressive phenotype [45]. Epithelial-mesenchymal transition (EMT) is promoted by hypoxia to facilitate invasion, migration and extravasation of CRC cells [44]. ① SPP1 is considered to be related to hypoxia, whereas the differential expression of SPP1 might indicate a hypoxic response. Osteopontin, encoded by SPP1, promotes the carcinogenesis, disease progression and recurrence of CRC [46, 47]. ② The hypoxia-associated gene TGFB1A is upregulation in the extracellular space, and this upregulation is related to immune evasion and immunological rejection [48, 49]. ③ The IGFBP3 gene is induced under hypoxic conditions; it regulates cell proliferation, senescence, apoptosis, and EMT and promotes the upregulation of a major cell surface receptor of hyaluronic acid, CD44 (CD44H). IGFBP3 inhibits ROS cytotoxicity to promote CD44 expression in the hypoxic microenvironment [49-52].

SJZD might inhibit EMT by altering the expression of TGFB1 to suppress the progression and influence the prognosis of CRC. The cellular EMT program is crucial for malignant progression; EMT is a coordinated process involving recombination of cell-extracellular matrix, protein binding, and regulation of drug resistance and tumor cell metastasis [53-55]. The hub gene TGFB1 has been shown to efficiently activate EMT. When EMT is activated, the intracellular and extracellular structures and interactions are altered; these changes include the disruption of cell-cell junctions and reorganization of the extracellular matrix [56]. EMT endows cancer cells with tumor invasion and metastasis capabilities [57]. Research has supported that the EMT program, which is mediated by EMT-related transcription factors, could transform nonmetastatic primary carcinoma cells into invasive and metastatic cells [58].

SJZD affects protein binding in the CRC TME by altering the expression of HSPB1, IGFBP-3 and SPP1. CRC tumor cells have high metastatic ability due to their low adhesion properties [59]. Protein binding mediates the response to intestinal injury to fuel regeneration [60]. ① HSPB1, also called HSP27, is a marker of metastatic tumors that partitions damaged proteins and prevents protein aggregation [61]. Overexpression of HSPB1 is associated with a poor prognosis in CRC [62]. ② Research has shown that IGFBP-3 cannot enter cells; thus, it interacts with surface proteins to accelerate tumor survival and invasion [63]. ③ SPP1 alters the normal expression pattern of intestinal epithelial cells by binding the IRF8 protein, and the expression of SPP1 is related to intestinal diseases such as ulcerative colitis [47, 64].

Nevertheless, there are still limitations and deficiencies in this study, such as the small sample size. We will explore the specific mechanism of these core genes in further research to provide a more accurate analysis of how the CRC TME affects prognosis.

Conclusions

In conclusion, our study explored core genes involved in the therapeutic mechanism of SJZD in the CRC TME. The results confirmed that the core genes HSPB1, SPP1, IGFBP3 and TGFB1 affected CRC development and prognosis by regulating hypoxia, protein binding and EMT in the extracellular matrix, which presents a new idea for the treatment of CRC.

Basic information of some active components of SJZD.

(DOCX)

Click here for additional data file.

The details of the GSEA results.

(XLSX)

Click here for additional data file.

The list of the DEGs.

(XLSX)

Click here for additional data file.

6 Dec 2021

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). You should upload this letter as a separate file labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. You should upload this as a separate file labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. You should upload this as a separate file labeled 'Manuscript'.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter. Guidelines for resubmitting your figure files are available below the reviewer comments at the end of this letter.

If applicable, we recommend that you deposit your laboratory protocols in protocols.io to enhance the reproducibility of your results. Protocols.io assigns your protocol its own identifier (DOI) so that it can be cited independently in the future. For instructions see: https://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols. Additionally, PLOS ONE offers an option for publishing peer-reviewed Lab Protocol articles, which describe protocols hosted on protocols.io. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols.

We look forward to receiving your revised manuscript.

Kind regards,

Ilya Ulasov, Ph.D

Academic Editor

PLOS ONE

Journal Requirements:

Additional Editor Comments (if provided):

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

Reviewer #2: Partly

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: No

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: No

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: The work is written in normal English. Possesses undoubted novelty and relevance. The work is based on an adequate bioinformatics analysis using standard algorithms in the R programming language and publicly available databases. However, the use of these algorithms and databases has allowed new data on compounds and mechanisms of their potential action in CRC to be obtained. This has important fundamental and applied significance.

Reviewer #2: The manuscript “Assessment of the targeted effect of Sijunzi decoction on the colorectal cancer microenvironment via the ESTIMATE algorithm “ by Jiaxin Du et al. describes the bioinformatic identification of potential targets of Sijunzi decoction (SJZD) that might play a role in colorectal cancer (CRC) via the regulation of the tumor microenvironment (TME).

The authors utilize publicly available gene expression and clinical data from two CRC patient cohorts to identify genes that are differentially regulated between tumors with high or low levels of immune and stromal cell infiltration, and then integrate these expression data with SJZD compound-target networks to identify TME-related genes, or their protein products, respectively, that are potentially targeted by SJZD active compounds, thus potentially explaining the beneficial effects of adjuvant treatment of CRC patients with SJZD.

However, there are numorous issues regarding methodology, data presentation and their interpretation that make the manuscript not acceptable for publication.

Major points :

1.

One major weakness is the rather limited number of patient samples analyzed in the present study. Could the authors validate their findings in an independent patient cohort (e.g. TCGA-COAD) ?

2.

In the methods section, the description regarding the correlation analysis of the ESTIMATE score are described inadequately. What exactly is the “correlation analysis model”, or a “correlation curve” and how was the survival analysis performed ?

3.

The relationship between immune cell and stromal infiltration and patient prognosis has been described previously. Hence, the analysis does not provide much additional insight.

Figure 1 A and B: These do not depict “correlation analyses”. It is unclear by which statistical method the p-values were obtained.

Lines 201 ff : It is not convincing to conclude any effects on survival, if the difference is not statistically significant.

4.

A description regarding the purpose of the GSEA and conclusions regarding its results are lacking. The results are not put into context with regard to the authors´research aim. Currently, the authors name a number of enriched gene sets without any interpretation of the results.

In the text, the authors refer to 613/48 and 730/41 gene sets in Figs 2 and 3, but only show a select few gene sets in the actual figure. Are these the statistically most significant gene sets ? Based on which criteria have they benn selected ?

Figures 2 and 3 do not have any figure legends. Normalized enrichment scores and statistical significance should be provided for each depicted gene set. It should be described in more detail which two groups were actually compared for the GSEA and how.

5.

Actually, the identification of DEGs is only described for normal vs tumor samples in the methods section. The identification of immune and stromal high / low DEGs should also be described, and lists with the identified DEGs should be provided.

Figures 4 A and B display identical heatmaps.

Figures 4 C,B : the labelling requires extensive editing

6.

What is the rationale to determine the overlap between potential SJZD targets and genes differentially expressed between normal and tumor tissues, since the authors identified immune and stromal high / low DEGs before ? This is confusing and should be explained more clearly.

Line 263 : Which DEGs were analyzed ? How were the 1607 genes in Fig 5B obtained ? This should be described explicitly and mentioned in the figure legend.

Line 264 : should be rephrased : “The active genes of SJZD”.

Furthermore, the statement “The active genes of SJZD also had differential expression in the tumor purity groups.” is not supported by any figure. Actually, this is also the more crucial analysis, since the authors aim to identify SJZD targets related to the TME, and not targets upregulated in tumor vs normal tissue.

Line 277: “The functional enrichment analysis of the DEGs” is actually not providing any information if the genes in the intersection (the 52 core genes) are related to “immune response”. Furthermore, a conclusion is lacking.

7.

Fig. 7: X-axis not properly labelled. A figure legend is missing.

Line 290 : “cytokine microarray”. Where is this shown ?

8.

Line 89 : The claim “we provided novel evidence that clarifies the influence of SJZD on the CRC TME” is an overstatement.

From this manuscript, it is actually not clear if SJZD has a defined (the authors only speak of "altered" expression in the discussion section) and measurable effect on the activity of any of the identified genes, or on the TME as such, in CRC.

Minor points :

1. Line 20 : please rephrase “Sijunzi decoction (SJZD) was used to treat patients with colorectal cancer (CRC) as an adjuvant method . The aim of the study was to investigate the effects of SJZD on the immune microenvironment of CRC”. The authors actually do not study the effects of SJZD on the tumor microenvironment, but rather aim to determine potential targets of SJZD related to the microenvironment using bioinformatic methods

2. Line 41 : delete “SJZD”

3. Lines 73f : please rephrase : “ could improve the survival and quality of life of patients with CRC by preventing tumorigenesis”

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Denis S Kutilin

Reviewer #2: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email PLOS at figures@plos.org. Please note that Supporting Information files do not need this step.

13 Jan 2022

Dear editor,

Thank you! We are grateful to the reviewers for the valuable comments and suggestions. We also thank the editor for an opportunity to resubmit the manuscript (PONE-D-21-21615) after modifications. We have studied comments carefully and have made correction which we hope meet with approval. We listed the changes and marked in red in revised paper. Our point-by-point answers and explanations to the reviewers’ comments are as follows:

Question #1: One major weakness is the rather limited number of patient samples analyzed in the present study. Could the authors validate their findings in an independent patient cohort (e.g. TCGA-COAD) ?

Answer 1: Thank you for your suggestions! We performed further survival analysis validation by using 530 colon adenocarcinoma (COAD) samples from the TCGA database. The TCGA samples further confirmed the poor prognosis of HSPB1, SPP1, IGFBP3 and TGFB1. COL3A1 and OLR1 showed no significance in this progress, so the previous 6 core prognostic targets (HSPB1, COL3A1, SPP1, OLR1, IGFBP3 and TGFB1) were revised to 4 core prognostic targets (HSPB1, SPP1, IGFBP3 and TGFB1). Please see the Method (Line 188 - 193) and Result (Line 313 - 316) section.

Questions #2-3:

2.In the methods section, the description regarding the correlation analysis of the ESTIMATE score are described inadequately. What exactly is the “correlation analysis model”, or a “correlation curve” and how was the survival analysis performed?

3.The relationship between immune cell and stromal infiltration and patient prognosis has been described previously. Hence, the analysis does not provide much additional insight.

Figure 1 A and B: These do not depict “correlation analyses”. It is unclear by which statistical method the p-values were obtained.

Lines 201 ff : It is not convincing to conclude any effects on survival, if the difference is not statistically significant.

Answer 2-3: Sorry for our incorrect writing! We inappropriately described Figure 1 A and B as "correlation analysis models".We amended these sentences. Please see the part of methods (Line 119-122). For survival analysis, a more detailed introduction has also been added in the methods section (Line 124-133). In line 222 the sentence has been corrected as “while the patients showed no significant difference of OS between different immune groups”.

Question #4:

A description regarding the purpose of the GSEA and conclusions regarding its results are lacking. The results are not put into context with regard to the authors´research aim. Currently, the authors name a number of enriched gene sets without any interpretation of the results.

In the text, the authors refer to 613/48 and 730/41 gene sets in Figs 2 and 3, but only show a select few gene sets in the actual figure. Are these the statistically most significant gene sets ? Based on which criteria have they benn selected ?

Figures 2 and 3 do not have any figure legends. Normalized enrichment scores and statistical significance should be provided for each depicted gene set. It should be described in more detail which two groups were actually compared for the GSEA and how.

Answer 4: Thank you for your reminders. We added the method and purpose of GSEA in Line 135~143. Corresponding interpretation of the results can be seen in the results section "Results of Gene Set Enrichment Analysis". Enriched gene sets with FDR < 0.25, |normalized enrichment score (NES)| > 1, nominal P < 0.05 were regarded statistically significant [32]. The results of GSEA have been re-visualized and the top 5 gene sets (sort by NES) were shown in Figure 2. New figure legend has updated in figure 2, and reloaded online. More details of GSEA result were profiled in the S1 Appendix.

[32] Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proceedings of the National Academy of Sciences of the United States of America. 2005;102(43):15545-50. Epub 2005/10/04. doi: 10.1073/pnas.0506580102. PubMed PMID: 16199517; PubMed Central PMCID: PMCPMC1239896.

Question #5:

5.Actually, the identification of DEGs is only described for normal vs tumor samples in the methods section. The identification of immune and stromal high / low DEGs should also be described, and lists with the identified DEGs should be provided.

Figures 4 A and B display identical heatmaps.

Figures 4 C,B : the labelling requires extensive editing

Answer 5: Sorry for our incorrect writing! In fact, we didn't want to identify the DEGs between the normal and tumor samples. The incorrect description of "normal vs tumor" in the methods has been modified to "stromal high vs stromal low or immune high vs immune low" and more details have been added (Line169-176). Lists with the identified DEGs can be seen in S2 Appendix. The modified heatmaps and labels were shown in Figure 3 now.

Question #6:

6.What is the rationale to determine the overlap between potential SJZD targets and genes differentially expressed between normal and tumor tissues, since the authors identified immune and stromal high / low DEGs before ? This is confusing and should be explained more clearly.

Line 263 : Which DEGs were analyzed ? How were the 1607 genes in Fig 5B obtained ? This should be described explicitly and mentioned in the figure legend.

Line 264 : should be rephrased : “The active genes of SJZD”.

Furthermore, the statement “The active genes of SJZD also had differential expression in the tumor purity groups.” is not supported by any figure. Actually, this is also the more crucial analysis, since the authors aim to identify SJZD targets related to the TME, and not targets upregulated in tumor vs normal tissue.

Line 277: “The functional enrichment analysis of the DEGs” is actually not providing any information if the genes in the intersection (the 52 core genes) are related to “immune response”. Furthermore, a conclusion is lacking.

Answer 6: Thank you for your suggestions! The overlap between potential SJZD targets and stromal high / low DEGs was what we actually did. But we incorrectly described it as "overlap between potential SJZD targets and genes differentially expressed between normal and tumor tissues". This error has been corrected in answer 5. The origin of 1557 genes in Fig 5B was given in line 254-260. The incorrect sentences and figures have been rectified in the manuscript. The correct figures and the details of the DEGs have uploaded in the platform (Fig 3 & S2 Appendix).

Question #7:

Fig. 7: X-axis not properly labelled. A figure legend is missing.

Line 290 : “cytokine microarray”. Where is this shown ?

Answer 7: Sorry for our mistake. The statements of “cytokine microarray” was a clerical error and has been deleted. The correct figure has re-uploaded as Fig 6.

Question #8:

Line 89 : The claim “we provided novel evidence that clarifies the influence of SJZD on the CRC TME” is an overstatement. From this manuscript, it is actually not clear if SJZD has a defined (the authors only speak of "altered" expression in the discussion section) and measurable effect on the activity of any of the identified genes, or on the TME as such, in CRC.

Answer 8: Thanks for your kind suggestion. In Line 86, the statement has been adjusted that “The therapeutic targets of SJZD and the DEGs between different infiltration levels were intersected, and the core genes affected the TME of CRC were found. The result will provide the TME of CRC with a novel auxiliary therapeutic perspective.”

Question #9:

1. Line 20 : please rephrase “Sijunzi decoction (SJZD) was used to treat patients with colorectal cancer (CRC) as an adjuvant method. The aim of the study was to investigate the effects of SJZD on the immune microenvironment of CRC”. The authors actually do not study the effects of SJZD on the tumor microenvironment, but rather aim to determine potential targets of SJZD related to the microenvironment using bioinformatic methods.

2. Line 41 : delete “SJZD”

3. Lines 73f : please rephrase : “ could improve the survival and quality of life of patients with CRC by preventing tumorigenesis”

Answer 9: Thank you for your reminding. The objective of this study has been corrected in Line 21. In Line 71, the statement has been rephrased by “It has been widely reported that TCM could prolong the survival time and improve the quality of life of patients with CRC by preventing tumorigenesis, suppressing tumor growth, and reducing metastasis and recurrence”.

Kind regards

Xin-lin Chen

Submitted filename: Response to Reviewers.docx

Click here for additional data file.

16 Feb 2022

Assessment of the targeted effect of Sijunzi decoction on the colorectal cancer microenvironment via the ESTIMATE algorithm

PONE-D-21-21615R1

Dear Dr. XIN-lin Chen,

We’re pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it meets all outstanding technical requirements.

Within one week, you’ll receive an e-mail detailing the required amendments. When these have been addressed, you’ll receive a formal acceptance letter and your manuscript will be scheduled for publication.

An invoice for payment will follow shortly after the formal acceptance. To ensure an efficient process, please log into Editorial Manager at http://www.editorialmanager.com/pone/, click the 'Update My Information' link at the top of the page, and double check that your user information is up-to-date. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to help maximize its impact. If they’ll be preparing press materials, please inform our press team as soon as possible -- no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

Kind regards,

Ilya Ulasov, Ph.D

Academic Editor

PLOS ONE

10 Mar 2022

PONE-D-21-21615R1

Assessment of the targeted effect of Sijunzi decoction on the colorectal cancer microenvironment via the ESTIMATE algorithm

Dear Dr. Chen:

I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

If we can help with anything else, please email us at plosone@plos.org.

Thank you for submitting your work to PLOS ONE and supporting open access.

Kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Ilya Ulasov

Academic Editor

PLOS ONE