# High-Throughput Screen of Natural Compounds and Biomarkers for NSCLC Treatment by Differential Expression and Weighted Gene Coexpression Network Analysis (WGCNA)

Authors: Ling Kui; Min Li; Xiaonan Yang; Ling Yang; Qinghua Kong; Yunbing Pan; Zetan Xu; Shouling Wang; Dandan Mo; Yang Dong; Yao Liu; Jianhua Miao
Journal: BioMed Research International
 (2021)

---

## Abstract

Lung cancer is known as the leading cause which presents the highest fatality rate worldwide; non-small-cell lung cancer (NSCLC) is the most prevalent type of lung carcinoma with high severity and affects 80% of patients with lung malignancies. Up to now, the general treatment for NSCLC includes surgery, chemotherapy, and radiotherapy; however, some therapeutic drugs and approaches could cause side effects and weaken the immune system. The combination of conventional therapies and traditional Chinese medicine (TCM) significantly improves treatment efficacy in lung cancer. Therefore, it is necessary to investigate the chemical composition and underlying antitumor mechanisms of TCM, so as to get a better understanding of the potential natural ingredient for lung cancer treatment. In this study, we selected 78 TCM to treat NSCLC cell line (A549) and obtained 92 transcriptome data; differential expression and WGCNA were applied to screen the potential natural ingredient and target genes. The sample which was treated withA. pierreana generated the most significant DEG set, including 6130 DEGs, 2479 upregulated, and 3651 downregulated. KEGG pathway analyses found that four pathways (MAPK, NF-kappa B, p53, and TGF-beta signaling pathway) were significantly enriched; 16 genes were significantly regulated in these four pathways. Interestingly, some of them such as EGFR, DUSP4, IL1R1, IL1B, MDM2, CDKNIA, and IDs have been used as the target biomarkers for cancer diagnosis and therapy. In addition, classified samples into 14 groups based on their pharmaceutical effects, WGCNA was used to identify 27 modules. Among them, green and darkgrey were the most relevant modules. Eight genes in the green module and four in darkgrey were identified as hub genes. In conclusion, we screened out three new TCM (B. fruticose, A. pierreana, and S. scandens) that have the potential to develop natural anticancer drugs and obtained the therapeutic targets for NSCLC therapy. Our study provides unique insights to screen the natural components for NSCLC therapy using high-throughput transcriptome analysis.


---

## Body

## 1. Introduction

Lung cancer remained the leading cause which presents the highest fatality rate worldwide, with an estimated 1.8 million deaths (18%) in each year reported by the International Agency for Research on Cancer/World Health Organization [1]. Particularly, lung cancer is the highest incidence of cancer and the high fatality rate in men, in 93 countries [2]. Approximately 85% of patients have a group of histological subtypes collectively known as non-small-cell lung cancer (NSCLC), which is the most prevalent type of lung carcinoma with high severity and affects 80% of patients with lung malignancies [3–5]. Up to now, the general treatment for NSCLC includes surgery, chemotherapy, and radiotherapy [6], but most chemotherapy drugs kill both cancer and normal cells [7–10]. The side effects of chemotherapies somewhat impact the quality of patients’ life. Therefore, it is urgent necessity to investigate the novel therapeutic strategies, to identify the prognostic markers based on better understanding of molecular mechanism of NSCLC. It helps to improve the precision rate for diagnosis and therapy and contributes to take the further steps for NSCLC research [3, 4].Traditional Chinese medicine (TCM) presented significant advantages in disease treatment special in cancer therapy; it was reported that the function of relieving adverse effects and enhancing the efficacy of drugs is presented [11–13]. Importantly, it has cytoprotective properties without hindering the anticancer activity of conventional drugs during combinational chemotherapy [14]. Finding active anticancer ingredients from TCM and reported their anticancer mechanisms became more and more popular for medical practitioners [15, 16]. A new trend in research of the new drug development for cancer treatment is obtaining natural ingredients with low toxicity and high efficiency from TCM [17, 18]. However, the anticancer compounds and molecular mechanisms of TCM are still unclear; the systemic methods need to be established. Undoubtedly, high-throughput screening is an effective and systematic approach to analyze TCM-mediated gene expression modifications in cells. Based on gene expression and pathway enrichment variations in different treatments, the underlying anticancer mechanisms of both conventional drugs and traditional medicines can be analyzed. Finally, it would indicate the potential TCM for lung cancer therapy according to the pharmaceutical effect and effective component comprehensively. Moreover, it may help to explain the molecular mechanism and screen the anticancer compounds, providing the new sight for the natural anticancer drug development on cancer treatment.The transcriptome is all the genes expressed in a certain cell in a certain functional status, and it can be used to compare differences in gene expression levels among various tissues or under various physiological conditions [19, 20]. The regulatory mechanisms for cancers are closely tied with differential genes and changes in signal pathways within the transcriptome. Comparative transcriptome analysis of cell lines with drug intervention allows the molecular mechanisms of cell growth conditions and physical sign changes following drug treatment to be holistically understood, and the interaction between drugs and cellular activities to be probed; the functional genes and their expression patterns can be analyzed, regulatory mechanisms mastered, and cancer-related signal pathways and key genes identified to find the target genes on which the drugs act [21–25]. These data can be used for targeted screening of medicinal ingredients, which can be verified and analyzed to find potential drugs and provide scientific, accurate, and comprehensive data for drug research and development [23].Weighted gene coexpression network analysis (WGCNA) is a bioinformatics data mining method that has been used to explore relationships between different gene modules of various cancer cell lines [26, 27]. The modules of coexpressing genes are found to maintain a consistent phenotype-independent expression relationship, and they may coregulate and share common biological functions [27]. Using WGCNA, expression alterations in gene sets, intrinsic properties of gene sets, correlation between gene modules, phenotype-correlated modules, candidate biomarker genes, and the targets for therapeutic drugs can be analyzed [28]. In previous studies, WGCNA was used to analyze biomarkers and targets of various diseases like schizophrenia, Alzheimer’s disease, sickle cell disease, and cancers [29–32]. As a systematical biological method, WGCNA has been used in various cancer studies including non-small-cell lung cancer (NSCLC), bladder cancer, clear cell renal cell carcinoma (ccRCC), acute myeloid leukemia (AML), and pancreatic ductal adenocarcinoma (PDAC) [28, 33–36].In this study, high-throughput transcriptome sequencing is used to provide the transcriptome data of NSCLC cell line distilled with 78 TCM and 10 chemical compounds (positive control) aforesaid. High-quality data sets are obtained from transcriptome assembly and comparative analysis. Then, identified and validated key genes were significantly associated with proinflammatory effects and metastasis process in NSCLC cancer by differential expression and WGCNA. We propose to establish a high-throughput method to screen potential natural compounds and identify the target genes, providing new sight on natural anticancer drug development for NSCLC cancer therapy.


## 2. Materials and Methods

### 2.1. Preparation of Medicinal Plant Extracts

All the 78 medicinal plants were provided by Guangxi Botanical Garden of Medicinal Plants. Plant materials were collected to dry and ground into fine powder for preparation of extracts. Some of the plant materials were boiled by hot water, then freeze dried; the obtained extracts were named as “W.” The extraction method named by “GX” refers to the medicinal plants extracted by 60% ethanol for 2 h. The extracted samples were vacuum concentrated and loaded on macroreticular resin column, and the product of interest was eluted by water. The eluted fractions were concentrated by vacuum evaporation, and the completely dried samples were used for drug preparation. The dried fruit ofMyristica fragrans was extracted using the CO2 supercritical extraction method, and the obtained extract was named as “C.” The remaining plant materials were firstly extracted by petroleum ether for 2 h and filtered, the remained residue was then extracted by ethyl acetate for another 2 h, and the extracts were vacuum concentrated to dry powder for drug preparation and named as “S.” More description of samples is shown in Supplemental Table 1. Ten conventional anticancer drugs purchased from various companies were set as positive control (Table 1). All the plant samples and drugs were dissolved in DMSO for further experiments, marked as fractions and positive control, respectively.Table 1
The list of positive control.
Drug
Source (company)
Sorafenib
Beijing Solarbio Science & Technology Co., Ltd
cis-Platinum
Beijing Solarbio Science & Technology Co., Ltd
Lenvatinib
Shanghai Topscience Bio-Technology Co., Ltd
Cabozantinib
Shanghai Yuanye Bio-Technology Co., Ltd
Doxorubicin hydrochloride
Shanghai Yuanye Bio-Technology Co., Ltd
Gefitinib
Shanghai Macklin Biological Co., Ltd
Paclitaxel
Shanghai Macklin Biological Co., Ltd
Docetaxel
Shanghai Macklin Biological Co., Ltd
Vinorelbine
Shanghai Aladdin Reagent Co., Ltd
Gemcitabine
Beijing OKA Bio-Technology Co., Ltd


### 2.2. Cell Culture and Drug Preparation

Purchased from the Shanghai Cell Bank of the Chinese Academy of Sciences, the A549 cell line was cultured in DMEM containing 10% fetal bovine serum and placed in a cell incubator with 5% CO2 at 37°C. Then, the cells with complete culture media when their growth density reached 80%-90% trypsin (0.125%) were subcultured. We selected logarithmic growth phase to test and screen the anticancer activity of A549 cell line, which was treated with 78 TCM and 10 conventional anticancer drugs (Supplemental Table 1, Supplemental Table 6, Table 1). The samples, which were treated with 10 conventional anticancer drugs, were considered as positive control, while 4 samples without any treatment were regarded as negative control. The fractions were extracted from the TCM species mainly including Leguminosae, Compositae, Araliaceae, Euphorbiaceae, and Rutaceae, all of which have been reported to have antitumor effects through previous studies. The positive control was conventional anticancer drugs purchased from the companies (Table 1).After the treatment experiments of fractions and positive control, a concentration of 100μg/ml and a final volume of 200 μl per well were prepared for first screening; then, the cell survival rate after 24 h of drug treatment was calculated. A 6-well plate with the same concentration was applied to screen the cells whose survival rate was more than 80%. If the survival rate was less than 80%, the 96-well plate screening was carried out again with a reduced concentration. A gradient of 2 times was set to screen the fractions; the screen concentration was set base on the cell survival rate during first screening. Finally, we selected the concentration at which the cell survival rate reached about 80% to carry out the 6-well plate screening. Subsequently, the cell growth conditions under a microscope were observed and this concentration as the final administration concentration was recorded.


### 2.3. RNA Quality Control, Library Preparation, and Sequencing

In our study, the FastPure Cell/Tissue Total RNA Isolation Mini Kit (Vazyme, Nanjing, China) was applied to extract the total RNA of the cell lines. NanoDrop 2000 and Agilent 2100 Bioanalyzer (Agilent, USA) were used to obtain the OD260/280; the value is 1.8~2.0. The RIN (RNA Integrity Number) values and concentration of all total RNA samples were quantified; RIN values reached 8. mRNA libraries were constructed by the MGIEasy kit; cDNA libraries with an insert size of 200-300 bp were prepared according to standard BGISEQ protocol with total RNA samples. Paired-end sequencing with 100 bp read length was sequenced using the BGISEQ-500 instrument and BGISEQ-500RS high-throughput sequencing kit.


### 2.4. RNA-seq Clean Data Preparation and Quality Checking

After sequencing, raw data were obtained in the fastq format. FastQC has been used to perform the quality control detecting the quality of sequencing data [37] according to the standards: filtering the low-quality reads with low base-call scores, the adapter sequences, the reads including N, and a shift from the expected GC content. The generated high-quality data were moved forward to the alignment step. Afterwards, the read sequences were trimmed and filtered by Trimmomatic software (v.0.36), which was included in the Trinity package [38, 39]; more detailed information about quality checking is shown in Supplemental Table 2; the generated clean data has been uploaded to NCBI Sequence Read Archive (SRA). Take multiple correlation analyses and screen the potential natural component for lung cancer therapy (Supplemental Table 1). Among these groups, groups 1-10 refer to the samples treated with TCM fractions, group 11 was mixed by samples treated with unclassified fractions, and groups 12 and 13 refer to positive controls and negative controls, respectively. Group 14 represents the samples treated with TCM which may include anticancer ingredient.


### 2.5. Alignment and Transcript Assembly

When quality control was finished, the general RNA-seq analyses would be carried out. There are four main steps that need to be done: aligning the reads to the reference; assembling the alignments on the alignment into a full-length transcript; quantitative expression of genes and transcripts; calculating the expression difference of all genes under different experimental conditions. The “new Tuxedo” package including HISAT, StringTie, and Ballgown has been used to perform this process. During this process, HISAT [40] has been used to align RNA-seq reads to the genome, and StringTie [41] is responsible for assembling transcripts and constructing isoforms to estimate gene expression. Ballgown [42] uses the results of StringTie splicing to calculate gene expression, then obtained the FPKM (Fragments Per Kilobase Million) results. The input data was generated by the BGISEQ-500 instrument; after running our pipeline, useful outputs were produced, including transcripts, gene expression values (FPKM), differentially expressed gene (DEG) list, and the merged statistical results. The detailed steps are shown in Table 2.Table 2
Detailed information of RNA-seq analysis pipeline.
Step
Analyses
Software
Script
Input
Main output
1
Prepare raw data
Trimmomatic ver.0.36
trimmomatic-0.36.jar PE -threads 12 -phred33 -trimlog xx.log xx.fq1 xx.fq2 -baseout
FASTQ file
FASTQ file
2
Quality control
FastQC
fastqc -t 10 -f fastq -o out xx.fq
FASTQ file from step 1
FASTQ file, multiqc_report.html
3
Alignment
HISAT
hisat2 -p 20 -x hg38.fa -1 xx_clean_1P.fq.gz -2 xx_clean_2P.fq.gz -S XX.sam
FASTQ file from step 2
Sam file
4
Sort
Samtools
samtools sort -@ 20 -O bam -o xx.bam xx.sam
Sam file from step 3
Bam file
5
Transcript assembly
StringTie
stringtie -e -p 30 -G hg38.gff -o xx.gtf xx.bam
Bam file from step 4
gtf file
6
Merged transcripts
StringTie-- merge
stringtie --merge -o merged.gtf gtflist
All gtf file from step 5
gtf file
7
Reads count
StringTie--eB
stringtie -B -e -p 30 -G merged.gtf -o xx.gtf xx.bam
gtf file from step 6 and Bam file from step 4
Gallgown input file and gtf file
8
Generated gene expression (FPKM)
Ballgown
gene_expression = gexpr(bg)
Output files from StringTie--eB
Gene expression tables
9
Differential expression analysis
DESeq2 Wrapper of TBtools
Selection criteria of DEGs:p.adjust≤0.05 and log2FC≥1
Gene expression tables
Volcano plot and differential gene list
10
KEGG pathway analysis
R package clusterProfiler andhttp://org.Hs.eg.dbdata set
Parameter:p.valueCutoff=0.05, p.AdjustMethod=“BH”, qvalueCutoff=0.1
Differential gene list
KEGG pathway tables
11
Gene expression visualization
TBtools
Default parameters
Expression of top genes
Heatmaps
12
General analysis of WGCNA
“WGCNA” R package
Power ofβ=7 (scale free R2=0.9)
FPKM table
Plots and forms of module connectivity
13
Selecting module and hub genes
“WGCNA” R package
cor
.
geneModuleMembership
>
0.8 and cor.geneTraitSignificance>0.2
Output files from WGCNA
Module and gene list
14
Network
Cytoscape
Default parametersThe criteria for edge filter: weight value>0.13 and weight value>0.06
Output files from WGCNA
Network plots
The table lists analysis steps, software, and main scripts in our pipeline. Starting from the input FASTQ files produced by sequencing and finally generating the results of candidate medicine and genes for NSCLC cancer research.


### 2.6. Differential Expression and KEGG Pathway Analysis

TBtools is a Toolkit for Biologists integrating various biological data-handling tools [43]; we applied it to perform differential expression analysis. Differentially expressed genes (DEGs) were filtered according to the criteria of p.adjust≤0.05 and ∣log2FC∣≥1, and a volcano plot was used to visualize the distribution of each gene. Subsequently, DEGs were imported to complete KEGG pathway analysis based on R package clusterProfiler and http://org.Hs.eg.dbdata set with the parameter of p.valueCutoff=0.05, p.AdjustMethod=“BH”, and qvalueCutoff=0.1. Based on KEGG results, the top 10 significant pathways and their enriched genes were selected as the candidate pathways and targets. They may be associated with inflammatory and metastasis processes in cancer. Furthermore, the expression of these genes was visualized by heat map through TBtools and the common genes were screened according to the high significance of differential expression. The analysis steps, software, and main scripts in our pipeline are listed in Table 2.


### 2.7. Weighted Gene Coexpression Network Analysis (WGCNA)

The “WGCNA” R package with default parameters was applied to construct the weighted coexpressed networks and identify the coexpression modules. The modified WGCNA pipeline used in this study is shown in Table2. The expression data generated by previous analyses were log2 normalized and imported into the “WGCNA” R package; sample information of 14 groups was summarized in the table. The soft thresholding was set according to the power of β=7 (scale free R2=0.9), and MEDissThres was set as 0.25 to merge similar modules, so as to ensure a scale-free network. Finally, edges were screened by the criteria of weight value, then inputted them into Cytoscape to visualize the coexpression network and identify the nodes and hub genes [44]. These genes, which may play an important role in tumorigenesis process, can be selected as the potential targets in the future research on cancer therapy.


## 2.1. Preparation of Medicinal Plant Extracts

All the 78 medicinal plants were provided by Guangxi Botanical Garden of Medicinal Plants. Plant materials were collected to dry and ground into fine powder for preparation of extracts. Some of the plant materials were boiled by hot water, then freeze dried; the obtained extracts were named as “W.” The extraction method named by “GX” refers to the medicinal plants extracted by 60% ethanol for 2 h. The extracted samples were vacuum concentrated and loaded on macroreticular resin column, and the product of interest was eluted by water. The eluted fractions were concentrated by vacuum evaporation, and the completely dried samples were used for drug preparation. The dried fruit ofMyristica fragrans was extracted using the CO2 supercritical extraction method, and the obtained extract was named as “C.” The remaining plant materials were firstly extracted by petroleum ether for 2 h and filtered, the remained residue was then extracted by ethyl acetate for another 2 h, and the extracts were vacuum concentrated to dry powder for drug preparation and named as “S.” More description of samples is shown in Supplemental Table 1. Ten conventional anticancer drugs purchased from various companies were set as positive control (Table 1). All the plant samples and drugs were dissolved in DMSO for further experiments, marked as fractions and positive control, respectively.Table 1
The list of positive control.
Drug
Source (company)
Sorafenib
Beijing Solarbio Science & Technology Co., Ltd
cis-Platinum
Beijing Solarbio Science & Technology Co., Ltd
Lenvatinib
Shanghai Topscience Bio-Technology Co., Ltd
Cabozantinib
Shanghai Yuanye Bio-Technology Co., Ltd
Doxorubicin hydrochloride
Shanghai Yuanye Bio-Technology Co., Ltd
Gefitinib
Shanghai Macklin Biological Co., Ltd
Paclitaxel
Shanghai Macklin Biological Co., Ltd
Docetaxel
Shanghai Macklin Biological Co., Ltd
Vinorelbine
Shanghai Aladdin Reagent Co., Ltd
Gemcitabine
Beijing OKA Bio-Technology Co., Ltd


## 2.2. Cell Culture and Drug Preparation

Purchased from the Shanghai Cell Bank of the Chinese Academy of Sciences, the A549 cell line was cultured in DMEM containing 10% fetal bovine serum and placed in a cell incubator with 5% CO2 at 37°C. Then, the cells with complete culture media when their growth density reached 80%-90% trypsin (0.125%) were subcultured. We selected logarithmic growth phase to test and screen the anticancer activity of A549 cell line, which was treated with 78 TCM and 10 conventional anticancer drugs (Supplemental Table 1, Supplemental Table 6, Table 1). The samples, which were treated with 10 conventional anticancer drugs, were considered as positive control, while 4 samples without any treatment were regarded as negative control. The fractions were extracted from the TCM species mainly including Leguminosae, Compositae, Araliaceae, Euphorbiaceae, and Rutaceae, all of which have been reported to have antitumor effects through previous studies. The positive control was conventional anticancer drugs purchased from the companies (Table 1).After the treatment experiments of fractions and positive control, a concentration of 100μg/ml and a final volume of 200 μl per well were prepared for first screening; then, the cell survival rate after 24 h of drug treatment was calculated. A 6-well plate with the same concentration was applied to screen the cells whose survival rate was more than 80%. If the survival rate was less than 80%, the 96-well plate screening was carried out again with a reduced concentration. A gradient of 2 times was set to screen the fractions; the screen concentration was set base on the cell survival rate during first screening. Finally, we selected the concentration at which the cell survival rate reached about 80% to carry out the 6-well plate screening. Subsequently, the cell growth conditions under a microscope were observed and this concentration as the final administration concentration was recorded.


## 2.3. RNA Quality Control, Library Preparation, and Sequencing

In our study, the FastPure Cell/Tissue Total RNA Isolation Mini Kit (Vazyme, Nanjing, China) was applied to extract the total RNA of the cell lines. NanoDrop 2000 and Agilent 2100 Bioanalyzer (Agilent, USA) were used to obtain the OD260/280; the value is 1.8~2.0. The RIN (RNA Integrity Number) values and concentration of all total RNA samples were quantified; RIN values reached 8. mRNA libraries were constructed by the MGIEasy kit; cDNA libraries with an insert size of 200-300 bp were prepared according to standard BGISEQ protocol with total RNA samples. Paired-end sequencing with 100 bp read length was sequenced using the BGISEQ-500 instrument and BGISEQ-500RS high-throughput sequencing kit.


## 2.4. RNA-seq Clean Data Preparation and Quality Checking

After sequencing, raw data were obtained in the fastq format. FastQC has been used to perform the quality control detecting the quality of sequencing data [37] according to the standards: filtering the low-quality reads with low base-call scores, the adapter sequences, the reads including N, and a shift from the expected GC content. The generated high-quality data were moved forward to the alignment step. Afterwards, the read sequences were trimmed and filtered by Trimmomatic software (v.0.36), which was included in the Trinity package [38, 39]; more detailed information about quality checking is shown in Supplemental Table 2; the generated clean data has been uploaded to NCBI Sequence Read Archive (SRA). Take multiple correlation analyses and screen the potential natural component for lung cancer therapy (Supplemental Table 1). Among these groups, groups 1-10 refer to the samples treated with TCM fractions, group 11 was mixed by samples treated with unclassified fractions, and groups 12 and 13 refer to positive controls and negative controls, respectively. Group 14 represents the samples treated with TCM which may include anticancer ingredient.


## 2.5. Alignment and Transcript Assembly

When quality control was finished, the general RNA-seq analyses would be carried out. There are four main steps that need to be done: aligning the reads to the reference; assembling the alignments on the alignment into a full-length transcript; quantitative expression of genes and transcripts; calculating the expression difference of all genes under different experimental conditions. The “new Tuxedo” package including HISAT, StringTie, and Ballgown has been used to perform this process. During this process, HISAT [40] has been used to align RNA-seq reads to the genome, and StringTie [41] is responsible for assembling transcripts and constructing isoforms to estimate gene expression. Ballgown [42] uses the results of StringTie splicing to calculate gene expression, then obtained the FPKM (Fragments Per Kilobase Million) results. The input data was generated by the BGISEQ-500 instrument; after running our pipeline, useful outputs were produced, including transcripts, gene expression values (FPKM), differentially expressed gene (DEG) list, and the merged statistical results. The detailed steps are shown in Table 2.Table 2
Detailed information of RNA-seq analysis pipeline.
Step
Analyses
Software
Script
Input
Main output
1
Prepare raw data
Trimmomatic ver.0.36
trimmomatic-0.36.jar PE -threads 12 -phred33 -trimlog xx.log xx.fq1 xx.fq2 -baseout
FASTQ file
FASTQ file
2
Quality control
FastQC
fastqc -t 10 -f fastq -o out xx.fq
FASTQ file from step 1
FASTQ file, multiqc_report.html
3
Alignment
HISAT
hisat2 -p 20 -x hg38.fa -1 xx_clean_1P.fq.gz -2 xx_clean_2P.fq.gz -S XX.sam
FASTQ file from step 2
Sam file
4
Sort
Samtools
samtools sort -@ 20 -O bam -o xx.bam xx.sam
Sam file from step 3
Bam file
5
Transcript assembly
StringTie
stringtie -e -p 30 -G hg38.gff -o xx.gtf xx.bam
Bam file from step 4
gtf file
6
Merged transcripts
StringTie-- merge
stringtie --merge -o merged.gtf gtflist
All gtf file from step 5
gtf file
7
Reads count
StringTie--eB
stringtie -B -e -p 30 -G merged.gtf -o xx.gtf xx.bam
gtf file from step 6 and Bam file from step 4
Gallgown input file and gtf file
8
Generated gene expression (FPKM)
Ballgown
gene_expression = gexpr(bg)
Output files from StringTie--eB
Gene expression tables
9
Differential expression analysis
DESeq2 Wrapper of TBtools
Selection criteria of DEGs:p.adjust≤0.05 and log2FC≥1
Gene expression tables
Volcano plot and differential gene list
10
KEGG pathway analysis
R package clusterProfiler andhttp://org.Hs.eg.dbdata set
Parameter:p.valueCutoff=0.05, p.AdjustMethod=“BH”, qvalueCutoff=0.1
Differential gene list
KEGG pathway tables
11
Gene expression visualization
TBtools
Default parameters
Expression of top genes
Heatmaps
12
General analysis of WGCNA
“WGCNA” R package
Power ofβ=7 (scale free R2=0.9)
FPKM table
Plots and forms of module connectivity
13
Selecting module and hub genes
“WGCNA” R package
cor
.
geneModuleMembership
>
0.8 and cor.geneTraitSignificance>0.2
Output files from WGCNA
Module and gene list
14
Network
Cytoscape
Default parametersThe criteria for edge filter: weight value>0.13 and weight value>0.06
Output files from WGCNA
Network plots
The table lists analysis steps, software, and main scripts in our pipeline. Starting from the input FASTQ files produced by sequencing and finally generating the results of candidate medicine and genes for NSCLC cancer research.


## 2.6. Differential Expression and KEGG Pathway Analysis

TBtools is a Toolkit for Biologists integrating various biological data-handling tools [43]; we applied it to perform differential expression analysis. Differentially expressed genes (DEGs) were filtered according to the criteria of p.adjust≤0.05 and ∣log2FC∣≥1, and a volcano plot was used to visualize the distribution of each gene. Subsequently, DEGs were imported to complete KEGG pathway analysis based on R package clusterProfiler and http://org.Hs.eg.dbdata set with the parameter of p.valueCutoff=0.05, p.AdjustMethod=“BH”, and qvalueCutoff=0.1. Based on KEGG results, the top 10 significant pathways and their enriched genes were selected as the candidate pathways and targets. They may be associated with inflammatory and metastasis processes in cancer. Furthermore, the expression of these genes was visualized by heat map through TBtools and the common genes were screened according to the high significance of differential expression. The analysis steps, software, and main scripts in our pipeline are listed in Table 2.


## 2.7. Weighted Gene Coexpression Network Analysis (WGCNA)

The “WGCNA” R package with default parameters was applied to construct the weighted coexpressed networks and identify the coexpression modules. The modified WGCNA pipeline used in this study is shown in Table2. The expression data generated by previous analyses were log2 normalized and imported into the “WGCNA” R package; sample information of 14 groups was summarized in the table. The soft thresholding was set according to the power of β=7 (scale free R2=0.9), and MEDissThres was set as 0.25 to merge similar modules, so as to ensure a scale-free network. Finally, edges were screened by the criteria of weight value, then inputted them into Cytoscape to visualize the coexpression network and identify the nodes and hub genes [44]. These genes, which may play an important role in tumorigenesis process, can be selected as the potential targets in the future research on cancer therapy.


## 3. Results

### 3.1. Data Preprocessing and Differentially Expressed Gene (DEG) Screening

In this study, 92 cDNA libraries were sequenced and produced about 806.67 Gb raw data, consisting of 8,066,731,626 reads with an average read length of 100 bp (Supplemental Table3). Clean sequencing reads (785.50 Gb and 7,855,046,204 reads) are available at NCBI Sequence Read Archive. After quality control, obtained sequencing reads were mapped to reference genome using STAR and sorted using Samtools. Subsequently, the gene expression (FPKM, Fragments Per Kilobase Million) was generated using RSEM. Differential expression analysis was performed by DESeq2 Wrapper included in TBTools; the selection criteria for differentially expressed genes (DEGs) are p.adjust≤0.05 and ∣log2FC∣≥1. We obtained 92 DEG set; the set contained most DEGs is Treat53 vs. CK, including 6130 DEGs, 2479 upregulated, and 3651 downregulated. We list the top 10 set shown in Tables 3, and selected the DEG set of three TCM candidates, which were Breynia fruticosa (B. fruticosa), Argyreia pierreana (A. pierreana), and Senecio scandens (S. scandens), illustrated by a volcano plot (Figures 1(b), 1(e), and 1(h)). The detailed DEG information is shown in Supplemental Table 4.Table 3
Differential gene expression of top 10 sets.
Compare
Significant
Upregulated
Downregulated
Sample name
Note
Treat53_vs_CK
6130
2479
3651
Argyreia pierreana
Candidate TCM
Treat91_vs_CK
1865
966
899
Excoecaria cochinchinensis Lour.
Treat8_vs_CK
1651
867
784
Breynia fruticosa
Candidate TCM
Treat56_vs_CK
1541
566
975
Tabernaemontana divaricata
Treat63_vs_CK
1162
550
612
Senecio scandens
Candidate TCM
Treat6_vs_CK
956
350
606
Melia azedarach
Treat59_vs_CK
814
527
287
Phyllanthus reticulatus Poir. var. glaberMuell.-Age.
Treat44_vs_CK
745
266
479
Chromolaena odorata
Treat27_vs_CK
645
186
459
Sorafenib
Treat90_vs_CK
627
259
368
Excoecaria cochinchinensis Lour.Figure 1
Plant morphology, differentially expressed genes (DEGs), and KEGG pathways of samples treated with three candidate TCM. (a–c)Breynia fruticosa, (d–f) Argyreia pierreana, and (g–i) Senecio scandens. Plant morphology pictures cited from Flora of China (source: http://www.iplant.cn). In volcano plots, the horizontal line at ∣log2FC∣=1; vertical line at false discovery rateFDR=0.05. In bar plots of KEGG pathways, red block represented the key pathways.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)


### 3.2. KEGG Pathway Analysis and Candidate TCM Screening

In order to identify biologic pathways, functional categories, and networks of DEGs, we used the R package clusterProfiler andhttp://org.Hs.eg.dbdata set to complete KEGG pathway analysis. The results indicated four key pathways associated with the biological processes of proinflammatory effects, and the metastasis process is significantly enriched, including the MAPK signaling pathway, NF-kappa B signaling pathway, p53 signaling pathway, and TGF-beta signaling pathway (Figures 1(c), 1(f), and 1(i)). Among all, the MAPK signaling pathway was enriched by the three plant samples B. fruticosa, A. pierreana, and S. scandens. The NF-kappa B signaling pathway was enriched by B. fruticosa and S. scandens. It indicated these three TCM may significantly regulate the expression of key genes in MAPK and NF-kappa B signaling pathways; they may include the potential ingredient for NSCLC therapy and should be screened as the effective TCM candidates (Figure 1). On the other hand, the p53 signaling pathway was enriched by 4 TCM samples and 5 positive control samples, while the TGF-beta signaling pathway was enriched by 8 TCM samples. It is still unknown whether they are associated with NSCLC or not.The detailed pathway results of other TCM are shown in Supplemental Table5.Additionally,B. fruticosa possesses the pharmaceutical effect of rheumatism treatment, A. pierreana and S. scandens possess the pharmaceutical effect of heat-clearing (Supplemental Table 6), and it is worth to explore the potential mechanism of NSCLC with rheumatism treatment and heat-clearing effect. We investigated and surveyed previous research; excitingly, B. fruticosa was reported including anti-inflammatory effect in cancer, which is one of the main factors with carcinogenesis. B. fruticosa widely grows in South China and was used to cure chronic bronchitis, sore throat, wounds, and gastroenteritis and presented an anticancer effect [45, 46]. B. fruticosa contains almost 10 bioactive compounds; among them, zizyberanalic acid and isoceanothic acid possess strong cytotoxic activity against five human cancer cell lines [46]. Meanwhile, it reported that the main active constituents of B. fruticosa include tannins, triterpenes, sterol derivatives, and lignins. Some compounds such as breynins presented the pharmacologic action of reducing inflammation through inhibiting NF-kappa B DNA-binding activity and expression of iNOS and COX-2; it is similar with thiacremonone [47]. Furthermore, it reported that the chloroplast genome of B. fruticosa has been completed; it helps to take the further research on genomic resources of Breynia species [45]. S. scandens is known as “Qianliguang” in Chinese; the active constituents include flavonoids, terpenes, alkaloids, carotenoids, volatile oils, chlorogenic acids, phenolic acids, and vitamins [48]. Among them, the main constituent is PAs and the main typical ingredients are adonifoline [48, 49]. Additionally, previous studies stated that some components isolated from Senecio scandens presented various pharmacological activities, such as antitumoral, anti-inflammatory, mutagenic, antiviral, antioxidant, and abirritation activities [49, 50]. However, there is almost no research on A. pierreana; we will plan to study on its pharmacological activities especially anticancer in the next research.


### 3.3. Expression and Regulation of Top Genes Enriched on These Key Pathways

In this study, we list the genes which were regulated significantly in MAPK and NF-kappa B signaling pathways in three candidate TCM, as well as other TCM enriched on p53 and TGF-beta signaling pathway (Figure2). The common significantly regulated genes in three candidate TCM associated with the MAPK signaling pathway are IL1R1, DUSP4, EGFR, EREG, and MAP3K20. Among them, IL1R1 was downregulated and DUSP4 was upregulated in all three candidate TCM; other genes present different regulations. In the NF-kappa B signaling pathway, BIRC3, TRIM25, TAB3, PLAU, and IL1B are common significantly upregulated; only IL1R1 was downregulated, the same as in the MAPK signaling pathway. MDM2 was upregulated and CDKN1A was downregulated in the p53 signaling pathway; ID1, ID2, and ID4 present different regulations in the TGF-beta signaling pathway. The regulated results of the common genes are shown in Table 4; the expression of common genes in all the samples is shown in Figure 2.Figure 2
Expression and regulation of top genes enriched on the key pathways. (a) MAPK signaling pathway; (b) NF-kappa B signaling pathway; (c) p53 signaling pathway; (d) TGF-beta signaling pathway. CK: the negative control.
(a)
(b)
(c)
(d)Table 4
The regulated results of the common genes.
Pathway
Common genes in medicine
Regulate
MAPK signaling pathway
IL1R1
Downregulated
DUSP4
Upregulated
EGFR
Upregulated inB. fruticosa and A. pierreana, downregulated in S. scandens
EREG
Upregulated inB. fruticosa and A. pierreana, downregulated in S. scandens
MAP3K20
Upregulated inB. fruticosa and A. pierreana, downregulated in S. scandens
NF-kappa B signaling pathway
IL1R1
Downregulated
BIRC3
Upregulated
TRIM25
Upregulated
TAB3
Upregulated
PLAU
Upregulated
IL1B
Upregulated
p53 signaling pathway
MDM2
Upregulated
CDKN1A
Upregulated
TGF-beta signaling pathway
ID2
Downregulated
ID1
Upregulated inA. pierreana, downregulated in others
ID4
DownregulatedIt is worth noting that some of the common genes were demonstrated to play an important role in inflammatory or metastasis processes in cancers, especially in NSCLC. What is more, some of them interacted on the key pathways, such as EGFR, DUSP4, IL1R1, IL1B, MDM2, CDKNIA, and IDs; they all have been used as the target biomarkers for cancer diagnosis and therapy. The result indicated that the TCM candidates significantly inhibited the NSCLC by regulating various common genes of different signaling pathways. Overall, this study inferred that our TCM candidates probably possessed the significant curative effects for NSCLC. In the subsequent studies, identification of effective compounds and functional investigation of common genes would be considered.As shown in Figure2 and Table 4, IL1R1 was included in both MAPK and NF-kappa B signaling pathways; it is IL-1 receptor type I. IL-1 is a master regulator which contributed to inflammation, hematopoiesis, and innate immunity. IL-1α and IL-1β are the representative proinflammatory cytokines, and both of them bonded to IL1R1 which also possesses proinflammatory effects. Besides, nuclear factor-kB (NF-kappa B), p38, and MAPKs are the key transcription factors involved in inflammatory and immune response activities [51, 52]. The epidermal growth factor receptor (EGFR) is one of the mutated signaling proteins in NSCLC, and the EGFR gene is one of the first molecules to be selected for targeted gene therapy [53]. EGFR is the kinase inhibitor and has successfully been used in targeted therapies of various oncogenic driver mutations [54]. DUSP4 is a negative regulator of the MAPK pathway [55]. It reported that the loss of complex heterozygosity between DOK2 and DUSP4 leads to the occurrence of lung adenocarcinomas, with a short incubation period and high incidence rate. Their codeletion can activate MAPK signaling and promote cell proliferation [56]. On the other hand, DUSP4 is involved in negative feedback control of EGFR signaling, and it proved the loss of DUSP4 associated with p16/CDKN2A deletion in the lung cancer patients [57].MDM2 as an oncoprotein and a major negative regulator of the p53 pathway possesses the function of activating the p53 target gene and target p53 protein for the degradation of proteasome. A previous study proved that chromatin-bound MDM2 also plays the p53-independent role to control the transcriptional genes associated in cell fate and metabolism [58]. MDM2 also presented the function of promoting degradation for ubiquitination and proteasomal dependence of wild-type p53, which is the regulator associated with cell cycle, senescence, apoptosis, DNA repair, and angiogenesis-related pathways. Xing et al. found that TNFAIP8 can regulate p53, MDM2, and cyclin D1 to induce cell proliferation and tumor growth in NSCLC [59]. The cyclin-dependent kinase inhibitor 1A (CDKN1A) gene is reported associated with drug (e.g., gefitinib) resistance and regulated the cell cycle in cancer, as its function of involving cell differentiation, DNA repair, and apoptosis. Besides, CDKN1A’s activities are closely related to p53 status. For example, when p53 lost mutate function, CDKN1A overexpression will induce cells to present the aggressive phenotype to avoid cell block, senescence, and apoptosis. In Zamagni et al.’s study, it found that CDKN1A is an oncogene which can inhibit apoptosis and promote cancer cell proliferation, so it is used as a response indicator of NSCLC chemotherapy; the regulation of CDKN1A is a potential therapy to reverse the acquisition of resistance to drug such as gefitinib. What is more, for KRAS- and TP53-mutated NSCLC, CDKN1A also can be used as a predictive biomarker of response [60–62].ID belongs to the helix–loop–helix (HLH) transcription factor family, was known as the inhibitor of differentiation/DNA-binding, and reported that it is involved in tumorigenesis, angiogenesis, and invasiveness. The ID family includes four members: ID1, ID2, ID3, and ID4. ID family member proteins contributed to proliferation, invasion, differentiation, metastasis, apoptosis, and angiogenesis in various human cancers; they may provide new targets and biomarkers for the treatment and prognosis of lung cancer [63–65]. In this study, we found ID1, ID2, and ID4 present significant differences; they may be involved in the physiological activity of NSCLC cells. There are plenty of evidences to support it: Rollin et al. found that ID2 may be developed as the biomarker for the prognosis of poorly differentiated tumors as it is associated with dedifferentiation [5]; Li et al. demonstrated that ID1 may activate the NF-kappa B signaling pathway to promote the process of proliferation, migration, and invasion in NSCLC cells [66]; Qi et al. found that ID4 impacted on the p38 MAPK pathway to inhibit cisplatin-induced apoptosis [67]. These two pathways also were proved very important in NSCLC in our study.


### 3.4. Coexpression Network Constructing and Significant Modules Identifying

In order to investigate the connectivity between the pharmaceutical effect of TCM fractions and module eigengene (ME), the “WGCNA” R package has been applied to construct weighted coexpressed networks and identify coexpression modules. In this research, the power ofβ=7 (scale free R2=0.9) was set to ensure low mean connectivity and high scale independence (Figure 3(a)), setting the dissimilarity of modules as 0.2 and generating 27 modules in total (Figure 4(a)). The module trait relationship is shown in Figure 4(b), and group 3, group 5, and group 12 were exhibited higher connectivity in several modules. Group 3 and group 5 were treated with TCM; group 12 was treated with conventional anticancer drugs (positive control). In TCM groups, the colors of green and darkgrey modules were the deepest; it suggested that it is suitable to identify the hub genes from these two modules which may be associated with the staging of cancer. The interaction relationship of 27 modules is shown in Figure 3(b); it indicated a high-scale independence and differential gene expressions between the modules as all the modules were independent with each other. 27 modules enriched into five clusters based on the eigengene adjacency heat map (Figure 3(c)).Figure 3
Estimation of soft-thresholding values, network heat map plot, and clustered correlation. (a) Scale independence and mean connectivity, power ofβ=7 (scale free R2=0.9). (b) The topological overlap heat map of coexpression genes summarized the interaction relationship of selected genes. Each module is represented by different colors. The brightness of yellow and red represented the high connectivity of modules. (c) Module eigengene dendrogram and eigengene adjacency heat map presented the relationship of the modules generated by the clustering analysis.
(a)
(b)
(c)Figure 4
The gene set enrichment analysis and the identification of modules. (a) Dendrogram of dissimilarity between all filtered gene sets, which enriched according to a dissimilarity measure (1-TOM), and the cluster modules were marked as different colors. Dynamic tree cut algorithm generated the first set, then merged the correlated modules together. Each gene was represented by a branch, and each module was represented by one color. (b) Heat map of the correlation between the pharmaceutical effect of TCM fractions and MEs in NSCLC. The dark degree of the module color represented the significance of their relationship.
(a)
(b)


### 3.5. Hub Gene Identification in the Selected Modules

Generally, genes included in the coexpression modules and with high connectivity were selected as hub genes. In this study, 12 hub genes (8 in green module and 4 in darkgrey module) were obtained (Figure5, Supplemental Table 7). Among them, PLAU and DUSP4 were the common genes identified by differential analysis. What is more, PLAU was included in the MAPK signaling pathway and DUSP4 was in the NF-kappa B signaling pathway; both of them were identified as the key pathways in KEGG analysis. The edges signifying the correlations in the green module were filtered by the criteria that weight value>0.13; then, a total of 132 nodes were identified after importing to Cytoscape (Figure 5(a), Supplemental Table 6). Meanwhile, darkgrey module was filtered by a condition of the weight value>0.06; 81 nodes were identified (Figure 5(b), Supplemental Table 7).Figure 5
The coexpression network. (a) The coexpression network of the significant genes in the green module, including 132 nodes. (b) The coexpression network of the significant genes in the darkgrey module, including 81 nodes. The hub genes were represented in red color.
(a)
(b)


## 3.1. Data Preprocessing and Differentially Expressed Gene (DEG) Screening

In this study, 92 cDNA libraries were sequenced and produced about 806.67 Gb raw data, consisting of 8,066,731,626 reads with an average read length of 100 bp (Supplemental Table3). Clean sequencing reads (785.50 Gb and 7,855,046,204 reads) are available at NCBI Sequence Read Archive. After quality control, obtained sequencing reads were mapped to reference genome using STAR and sorted using Samtools. Subsequently, the gene expression (FPKM, Fragments Per Kilobase Million) was generated using RSEM. Differential expression analysis was performed by DESeq2 Wrapper included in TBTools; the selection criteria for differentially expressed genes (DEGs) are p.adjust≤0.05 and ∣log2FC∣≥1. We obtained 92 DEG set; the set contained most DEGs is Treat53 vs. CK, including 6130 DEGs, 2479 upregulated, and 3651 downregulated. We list the top 10 set shown in Tables 3, and selected the DEG set of three TCM candidates, which were Breynia fruticosa (B. fruticosa), Argyreia pierreana (A. pierreana), and Senecio scandens (S. scandens), illustrated by a volcano plot (Figures 1(b), 1(e), and 1(h)). The detailed DEG information is shown in Supplemental Table 4.Table 3
Differential gene expression of top 10 sets.
Compare
Significant
Upregulated
Downregulated
Sample name
Note
Treat53_vs_CK
6130
2479
3651
Argyreia pierreana
Candidate TCM
Treat91_vs_CK
1865
966
899
Excoecaria cochinchinensis Lour.
Treat8_vs_CK
1651
867
784
Breynia fruticosa
Candidate TCM
Treat56_vs_CK
1541
566
975
Tabernaemontana divaricata
Treat63_vs_CK
1162
550
612
Senecio scandens
Candidate TCM
Treat6_vs_CK
956
350
606
Melia azedarach
Treat59_vs_CK
814
527
287
Phyllanthus reticulatus Poir. var. glaberMuell.-Age.
Treat44_vs_CK
745
266
479
Chromolaena odorata
Treat27_vs_CK
645
186
459
Sorafenib
Treat90_vs_CK
627
259
368
Excoecaria cochinchinensis Lour.Figure 1
Plant morphology, differentially expressed genes (DEGs), and KEGG pathways of samples treated with three candidate TCM. (a–c)Breynia fruticosa, (d–f) Argyreia pierreana, and (g–i) Senecio scandens. Plant morphology pictures cited from Flora of China (source: http://www.iplant.cn). In volcano plots, the horizontal line at ∣log2FC∣=1; vertical line at false discovery rateFDR=0.05. In bar plots of KEGG pathways, red block represented the key pathways.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)


## 3.2. KEGG Pathway Analysis and Candidate TCM Screening

In order to identify biologic pathways, functional categories, and networks of DEGs, we used the R package clusterProfiler andhttp://org.Hs.eg.dbdata set to complete KEGG pathway analysis. The results indicated four key pathways associated with the biological processes of proinflammatory effects, and the metastasis process is significantly enriched, including the MAPK signaling pathway, NF-kappa B signaling pathway, p53 signaling pathway, and TGF-beta signaling pathway (Figures 1(c), 1(f), and 1(i)). Among all, the MAPK signaling pathway was enriched by the three plant samples B. fruticosa, A. pierreana, and S. scandens. The NF-kappa B signaling pathway was enriched by B. fruticosa and S. scandens. It indicated these three TCM may significantly regulate the expression of key genes in MAPK and NF-kappa B signaling pathways; they may include the potential ingredient for NSCLC therapy and should be screened as the effective TCM candidates (Figure 1). On the other hand, the p53 signaling pathway was enriched by 4 TCM samples and 5 positive control samples, while the TGF-beta signaling pathway was enriched by 8 TCM samples. It is still unknown whether they are associated with NSCLC or not.The detailed pathway results of other TCM are shown in Supplemental Table5.Additionally,B. fruticosa possesses the pharmaceutical effect of rheumatism treatment, A. pierreana and S. scandens possess the pharmaceutical effect of heat-clearing (Supplemental Table 6), and it is worth to explore the potential mechanism of NSCLC with rheumatism treatment and heat-clearing effect. We investigated and surveyed previous research; excitingly, B. fruticosa was reported including anti-inflammatory effect in cancer, which is one of the main factors with carcinogenesis. B. fruticosa widely grows in South China and was used to cure chronic bronchitis, sore throat, wounds, and gastroenteritis and presented an anticancer effect [45, 46]. B. fruticosa contains almost 10 bioactive compounds; among them, zizyberanalic acid and isoceanothic acid possess strong cytotoxic activity against five human cancer cell lines [46]. Meanwhile, it reported that the main active constituents of B. fruticosa include tannins, triterpenes, sterol derivatives, and lignins. Some compounds such as breynins presented the pharmacologic action of reducing inflammation through inhibiting NF-kappa B DNA-binding activity and expression of iNOS and COX-2; it is similar with thiacremonone [47]. Furthermore, it reported that the chloroplast genome of B. fruticosa has been completed; it helps to take the further research on genomic resources of Breynia species [45]. S. scandens is known as “Qianliguang” in Chinese; the active constituents include flavonoids, terpenes, alkaloids, carotenoids, volatile oils, chlorogenic acids, phenolic acids, and vitamins [48]. Among them, the main constituent is PAs and the main typical ingredients are adonifoline [48, 49]. Additionally, previous studies stated that some components isolated from Senecio scandens presented various pharmacological activities, such as antitumoral, anti-inflammatory, mutagenic, antiviral, antioxidant, and abirritation activities [49, 50]. However, there is almost no research on A. pierreana; we will plan to study on its pharmacological activities especially anticancer in the next research.


## 3.3. Expression and Regulation of Top Genes Enriched on These Key Pathways

In this study, we list the genes which were regulated significantly in MAPK and NF-kappa B signaling pathways in three candidate TCM, as well as other TCM enriched on p53 and TGF-beta signaling pathway (Figure2). The common significantly regulated genes in three candidate TCM associated with the MAPK signaling pathway are IL1R1, DUSP4, EGFR, EREG, and MAP3K20. Among them, IL1R1 was downregulated and DUSP4 was upregulated in all three candidate TCM; other genes present different regulations. In the NF-kappa B signaling pathway, BIRC3, TRIM25, TAB3, PLAU, and IL1B are common significantly upregulated; only IL1R1 was downregulated, the same as in the MAPK signaling pathway. MDM2 was upregulated and CDKN1A was downregulated in the p53 signaling pathway; ID1, ID2, and ID4 present different regulations in the TGF-beta signaling pathway. The regulated results of the common genes are shown in Table 4; the expression of common genes in all the samples is shown in Figure 2.Figure 2
Expression and regulation of top genes enriched on the key pathways. (a) MAPK signaling pathway; (b) NF-kappa B signaling pathway; (c) p53 signaling pathway; (d) TGF-beta signaling pathway. CK: the negative control.
(a)
(b)
(c)
(d)Table 4
The regulated results of the common genes.
Pathway
Common genes in medicine
Regulate
MAPK signaling pathway
IL1R1
Downregulated
DUSP4
Upregulated
EGFR
Upregulated inB. fruticosa and A. pierreana, downregulated in S. scandens
EREG
Upregulated inB. fruticosa and A. pierreana, downregulated in S. scandens
MAP3K20
Upregulated inB. fruticosa and A. pierreana, downregulated in S. scandens
NF-kappa B signaling pathway
IL1R1
Downregulated
BIRC3
Upregulated
TRIM25
Upregulated
TAB3
Upregulated
PLAU
Upregulated
IL1B
Upregulated
p53 signaling pathway
MDM2
Upregulated
CDKN1A
Upregulated
TGF-beta signaling pathway
ID2
Downregulated
ID1
Upregulated inA. pierreana, downregulated in others
ID4
DownregulatedIt is worth noting that some of the common genes were demonstrated to play an important role in inflammatory or metastasis processes in cancers, especially in NSCLC. What is more, some of them interacted on the key pathways, such as EGFR, DUSP4, IL1R1, IL1B, MDM2, CDKNIA, and IDs; they all have been used as the target biomarkers for cancer diagnosis and therapy. The result indicated that the TCM candidates significantly inhibited the NSCLC by regulating various common genes of different signaling pathways. Overall, this study inferred that our TCM candidates probably possessed the significant curative effects for NSCLC. In the subsequent studies, identification of effective compounds and functional investigation of common genes would be considered.As shown in Figure2 and Table 4, IL1R1 was included in both MAPK and NF-kappa B signaling pathways; it is IL-1 receptor type I. IL-1 is a master regulator which contributed to inflammation, hematopoiesis, and innate immunity. IL-1α and IL-1β are the representative proinflammatory cytokines, and both of them bonded to IL1R1 which also possesses proinflammatory effects. Besides, nuclear factor-kB (NF-kappa B), p38, and MAPKs are the key transcription factors involved in inflammatory and immune response activities [51, 52]. The epidermal growth factor receptor (EGFR) is one of the mutated signaling proteins in NSCLC, and the EGFR gene is one of the first molecules to be selected for targeted gene therapy [53]. EGFR is the kinase inhibitor and has successfully been used in targeted therapies of various oncogenic driver mutations [54]. DUSP4 is a negative regulator of the MAPK pathway [55]. It reported that the loss of complex heterozygosity between DOK2 and DUSP4 leads to the occurrence of lung adenocarcinomas, with a short incubation period and high incidence rate. Their codeletion can activate MAPK signaling and promote cell proliferation [56]. On the other hand, DUSP4 is involved in negative feedback control of EGFR signaling, and it proved the loss of DUSP4 associated with p16/CDKN2A deletion in the lung cancer patients [57].MDM2 as an oncoprotein and a major negative regulator of the p53 pathway possesses the function of activating the p53 target gene and target p53 protein for the degradation of proteasome. A previous study proved that chromatin-bound MDM2 also plays the p53-independent role to control the transcriptional genes associated in cell fate and metabolism [58]. MDM2 also presented the function of promoting degradation for ubiquitination and proteasomal dependence of wild-type p53, which is the regulator associated with cell cycle, senescence, apoptosis, DNA repair, and angiogenesis-related pathways. Xing et al. found that TNFAIP8 can regulate p53, MDM2, and cyclin D1 to induce cell proliferation and tumor growth in NSCLC [59]. The cyclin-dependent kinase inhibitor 1A (CDKN1A) gene is reported associated with drug (e.g., gefitinib) resistance and regulated the cell cycle in cancer, as its function of involving cell differentiation, DNA repair, and apoptosis. Besides, CDKN1A’s activities are closely related to p53 status. For example, when p53 lost mutate function, CDKN1A overexpression will induce cells to present the aggressive phenotype to avoid cell block, senescence, and apoptosis. In Zamagni et al.’s study, it found that CDKN1A is an oncogene which can inhibit apoptosis and promote cancer cell proliferation, so it is used as a response indicator of NSCLC chemotherapy; the regulation of CDKN1A is a potential therapy to reverse the acquisition of resistance to drug such as gefitinib. What is more, for KRAS- and TP53-mutated NSCLC, CDKN1A also can be used as a predictive biomarker of response [60–62].ID belongs to the helix–loop–helix (HLH) transcription factor family, was known as the inhibitor of differentiation/DNA-binding, and reported that it is involved in tumorigenesis, angiogenesis, and invasiveness. The ID family includes four members: ID1, ID2, ID3, and ID4. ID family member proteins contributed to proliferation, invasion, differentiation, metastasis, apoptosis, and angiogenesis in various human cancers; they may provide new targets and biomarkers for the treatment and prognosis of lung cancer [63–65]. In this study, we found ID1, ID2, and ID4 present significant differences; they may be involved in the physiological activity of NSCLC cells. There are plenty of evidences to support it: Rollin et al. found that ID2 may be developed as the biomarker for the prognosis of poorly differentiated tumors as it is associated with dedifferentiation [5]; Li et al. demonstrated that ID1 may activate the NF-kappa B signaling pathway to promote the process of proliferation, migration, and invasion in NSCLC cells [66]; Qi et al. found that ID4 impacted on the p38 MAPK pathway to inhibit cisplatin-induced apoptosis [67]. These two pathways also were proved very important in NSCLC in our study.


## 3.4. Coexpression Network Constructing and Significant Modules Identifying

In order to investigate the connectivity between the pharmaceutical effect of TCM fractions and module eigengene (ME), the “WGCNA” R package has been applied to construct weighted coexpressed networks and identify coexpression modules. In this research, the power ofβ=7 (scale free R2=0.9) was set to ensure low mean connectivity and high scale independence (Figure 3(a)), setting the dissimilarity of modules as 0.2 and generating 27 modules in total (Figure 4(a)). The module trait relationship is shown in Figure 4(b), and group 3, group 5, and group 12 were exhibited higher connectivity in several modules. Group 3 and group 5 were treated with TCM; group 12 was treated with conventional anticancer drugs (positive control). In TCM groups, the colors of green and darkgrey modules were the deepest; it suggested that it is suitable to identify the hub genes from these two modules which may be associated with the staging of cancer. The interaction relationship of 27 modules is shown in Figure 3(b); it indicated a high-scale independence and differential gene expressions between the modules as all the modules were independent with each other. 27 modules enriched into five clusters based on the eigengene adjacency heat map (Figure 3(c)).Figure 3
Estimation of soft-thresholding values, network heat map plot, and clustered correlation. (a) Scale independence and mean connectivity, power ofβ=7 (scale free R2=0.9). (b) The topological overlap heat map of coexpression genes summarized the interaction relationship of selected genes. Each module is represented by different colors. The brightness of yellow and red represented the high connectivity of modules. (c) Module eigengene dendrogram and eigengene adjacency heat map presented the relationship of the modules generated by the clustering analysis.
(a)
(b)
(c)Figure 4
The gene set enrichment analysis and the identification of modules. (a) Dendrogram of dissimilarity between all filtered gene sets, which enriched according to a dissimilarity measure (1-TOM), and the cluster modules were marked as different colors. Dynamic tree cut algorithm generated the first set, then merged the correlated modules together. Each gene was represented by a branch, and each module was represented by one color. (b) Heat map of the correlation between the pharmaceutical effect of TCM fractions and MEs in NSCLC. The dark degree of the module color represented the significance of their relationship.
(a)
(b)


## 3.5. Hub Gene Identification in the Selected Modules

Generally, genes included in the coexpression modules and with high connectivity were selected as hub genes. In this study, 12 hub genes (8 in green module and 4 in darkgrey module) were obtained (Figure5, Supplemental Table 7). Among them, PLAU and DUSP4 were the common genes identified by differential analysis. What is more, PLAU was included in the MAPK signaling pathway and DUSP4 was in the NF-kappa B signaling pathway; both of them were identified as the key pathways in KEGG analysis. The edges signifying the correlations in the green module were filtered by the criteria that weight value>0.13; then, a total of 132 nodes were identified after importing to Cytoscape (Figure 5(a), Supplemental Table 6). Meanwhile, darkgrey module was filtered by a condition of the weight value>0.06; 81 nodes were identified (Figure 5(b), Supplemental Table 7).Figure 5
The coexpression network. (a) The coexpression network of the significant genes in the green module, including 132 nodes. (b) The coexpression network of the significant genes in the darkgrey module, including 81 nodes. The hub genes were represented in red color.
(a)
(b)


## 4. Discussion and Conclusions

Guangxi Botanical Garden of Medicinal Plants (Nanning, China) has collected and conserved various plant species belonging to different families including Euphorbiaceae (Ricinus, Mallotus, Euphorbia, etc.), Araliaceae (Panax), Asteraceae (Acmella, Arctium, Gynura, etc.), Fabaceae, Malvaceae, and Solanaceae. There are almost 78 kinds of fractions (Supplemental Table 1) with remarkable medical effect, including promoting blood circulation and removing stasis; heat-clearing; rheumatism treatment; hemostatic; tonic; asthmatics, expectorants, and antitussives; water-disinhibiting and damp-percolating; Qi-regulating; heat-clearing astringent; and detoxifying, analgesia, antipruritic, etc. Some of the plant materials with potential anticancer effects were analyzed in previous studies, such as Senecio scandens, Breynia fruticosa, Salvia miltiorrhiza Bunge, Gynura procumbens, Euphorbia hirta, Forsythia suspensa, Senecio scandens, Polygonum perfoliatum, Malva verticillata var. crispa, Mallotus apelta, Pteris semipinnata, Acmella paniculate, Aristolochia tagala, Schisandra chinensis, Gelsemium elegans, and Ligustrum confusum [46, 47, 68–71]. Most of them were found to contain anticancer compounds like terpenoid, flavonoid, alkaloid, quercetin, epigallocatechin-3-gallate (EGCG), epicatechin, oleanolic acid, ursolic acid, and tanshinone, which are involved in multiple important pathways in cancers such as epithelial to mesenchymal transition (EMT), TGF-β, PTEN/PI3K/Akt, NF-kappa B, MAPK, p53 signaling pathway, and Wnt/β-catenin pathway [72–75].Some pathways which are mainly enriched by the TCM play an important role in cancer. Mitogen-activated protein kinase (MAPK) plays a crucial role in cellular signal transduction pathways as it involved not only cell proliferation, differentiation, and apoptosis but also the deregulation of cancer [76]. MAPK dysregulation may cause various cancer formations, including lungs, breast, oral, colorectal, ovarian, and thyroid [77]. It will improve the effectiveness of drugs through targeting MAPK pathways which contributed to conventional anticancer drugs [78]. The nuclear factor κB (NF-kappa B) complex was constructed from a family which induced transcription factors and can be found in almost all cells. Some inflammatory cytokine genes are IL-6, IL-8, and TNF-α; their expression would be improved by activating NF-kappa B [79]. It proved that NF-kappa B is involved in the initiation and progression of inflammation tumor tissues [80]; it even presented the function of harmonizing the key endogenous tumor promoter and inflammation. In addition to promoting tumor cell growth, NF-kappa B also activated the relative elements such as adhesion molecules, inflammatory cytokines, and angiogenic factors which expedited the cell proliferation in cancer [81, 82]. The NF-kappa B pathway has been regarded as the inflammation-mediated pathways, as it presented the capacity of escaping the apoptosis to cause emergency resistance to chemotherapy therapy [83].p53 is one of the most intensive tumor suppressor proteins, and TP53 is reported as the most commonly mutated gene in human cancer, and the mutations induced the high expression of mutant p53 proteins. Previous genomic studies indicated that the function of p53 was compromised by its altofrequency in different cancer cell lines [84]. The development of most cancers required perturbations in p53 signaling pathways; more and more evidences indicated that restorated or reactivated function of p53 will benefit to cancer therapy [85]. The p53 signaling pathway is the key pathway involved in cell division, cell proliferation, and signal transduction in NSCLC; Tu et al. proved it via comparative analysis of 221 DEGs and constructing the PPI network [86]. Transforming Growth Factor beta (TGF-beta) is a tumor-suppressive factor in early stage tumors, as aberrant activation of TGF-beta will induce angiogenesis, invasion, immunosuppression, and self-renewal of cancer-initiating [87]. TGF-beta is associated with poor prognosis and aggressive disease progression in NSCLC, as its ability of driving cell proliferation, metastasis, angiogenesis, emergence of drug resistance, and immune evasion. Some processes such as lung organogenesis, tissue remodeling after lung injury, and postnatal lung homeostasis all required the attendance of TGF-beta [54]; it is the reason why it has been reported as one of the most commonly indispensable and activated pathways in the metastasis process of various cancers. Particularly, activating TGF-β signaling induced progression and metastasis in NSCLC [54, 87, 88].In this research, we provide the better understanding and evidence to support the exploration of anticancer potential of the natural ingredient from TCM. Based on our results, we can probe the interaction between medicinal plant ingredients and NSCLC cells, to explore the differential gene activation mechanisms for pharmacodynamic substances to identify medicinal compositions and action mechanisms. Furthermore, we will be able to not only identify novel makers and screen out NSCLC therapeutic targets but also develop new natural drug groups against lung cancer, bringing new ideas for the development of new methods to screen out new drugs.



---

*Source: 5955343-2021-08-27.xml*