N6-methyladenosine Modification-Related Long Non-Coding RNAs are Potential Biomarkers for Predicting the Prognosis of Patients With Osteosarcoma.

Background: The role of N6-methyladenosine (m6A)-related long non-coding RNAs (lncRNAs) in osteosarcoma (OS) has not been fully studied yet. We aimed to identify m6A-related lncRNAs that could act as prognostic biomarkers for OS.
Methods: Pearson correlation was performed to identify m6A-related lncRNAs. Univariate and multivariate Cox regression analyses were performed to construct the risk model and assess whether the risk score was an independent prognostic factor for patients with OS. Gene Set Enrichment Analysis (GSEA) was performed to analyze the functions of genes in high-risk and low-risk groups. StarBase and Cytoscape were used to construct a competing endogenous RNA (ceRNA) network based on m6A-related prognostic lncRNA signature. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed to analyze the function of genes involved in the ceRNA network.
Results: We extracted 122 common lncRNAs from TCGA and Gene Expression Omnibus (GEO) databases. Pearson correlation results revealed 59 significant m6A-related lncRNAs in The Cancer Genome Atlas (TCGA) database, from which 2 were screened to construct a risk signature in TCGA dataset, which was then validated in the GEO dataset. A corresponding risk score was calculated and shown to be an independent prognostic factor for patients with OS. Enrichment analysis indicated that cell proliferation-related biological processes were more common in the high-risk group, while immune-related biological processes were more common in the low-risk group. Moreover, we established a nomogram that had a good ability to predict the overall survival of patients with OS. Additionally, a ceRNA network based on small nucleolar RNA host gene 7 (SNHG7) and small nucleolar RNA host gene 12 (SNHG12) was constructed, with genes that were enriched in hepatocellular carcinoma, gastric cancer, and non-small-cell lung cancer pathways.
Conclusion: Our study revealed the prognostic role of m6A-related lncRNAs in OS and identified SNHG7 and SNHG12 as potential biomarkers for predicting the prognosis of patients with OS. These findings have enriched our understanding of the role of m6A modification in the dysregulation of lncRNAs in OS.

Introduction

Osteosarcoma (OS) mainly occurs in the metaphysis of long bones and is one of the most malignant bone tumors in children and adolescents. OS shows a worldwide incidence of 3.4 per million people per year. Comprehensive therapies including surgical resection, neoadjuvant chemotherapy, and radiotherapy are currently used in the treatment of patients with OS. The 5-year survival rate of patients with OS has increased from 20% to 68%.[4,5] However, the overall survival rate has remained markedly unchanged owing to its potent invasiveness and high rate of local recurrence and distant metastasis. Furthermore, patients with OS and distant metastasis at the first clinical visit or chemoresistance fail to achieve a favorable 5-year survival rate. Therefore, new biomarkers and novel therapeutic targets based on the molecular mechanisms of OS need to be developed for early diagnosis and treatment.

Long non-coding RNAs (lncRNAs) are defined as transcripts larger than 200 nucleotides with little or no protein-coding ability. Accumulating evidence shows that lncRNAs function as pivotal regulators of gene expression at the level of chromatin modification, transcription, and posttranscriptional processing.[9,10] LncRNA alterations are associated with changes in biological processes and in the development of diseases.[11,12] Aberrant expression of lncRNAs contributes to carcinogenesis and progression of OS tumors. Multiple studies have shown that lncRNAs are differentially expressed in OS tissues and are implicated in its initiation, invasion, and migration.[13-15] Moreover, abnormally expressed lncRNAs can serve as valuable biomarkers in predicting the prognosis of patients with OS.[16,17] Deeper insights into the functions of lncRNAs hold great promise for the identification of novel therapeutic targets and development of effective drugs.

N6-methyladenosine (m6A) modification occurs at the N6 position of adenosine and has been identified as the most prevalent and conserved posttranscriptional modification of messenger RNAs (mRNAs). m6A modification is involved in regulating gene expression via modulating the processing, translation, splicing, and stability of mRNAs.[19-22] It is a dynamic and reversible process that is facilitated by writers (eg METTL3/14 and WTAP), reverted by erasers (eg ALKBH5 and FTO), and recognized by reader proteins (eg YTHDF1/2/3 and IGF2BP1). Alteration of m6A regulators is implicated in the initiation, migration, and progression of various cancers. Notably, recent studies have revealed that m6A modification is strongly associated with tumor metastasis and prognosis of patients with OS.[25-27] For instance, it was found that WTAP is highly expressed in OS tissues and is involved in the proliferation and metastasis of OS by inhibiting the expression of homeobox-containing 1 (HMBOX1) in an m6A-dependent manner. Interestingly, a considerable amount of research has been carried out on m6A modification of lncRNAs. The m6A-related lncRNAs may play critical roles in the pathogenesis of OS. However, few studies have focused on the prognostic value and potential mechanisms of m6A-related lncRNAs in OS.

In this study, we collected transcriptome data and clinicopathological features of patients with OS from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases and performed bioinformatics and statistical analyses to identify m6A-related lncRNAs and their prognostic value in OS. Furthermore, to enrich the regulatory network of prognostic m6A-related lncRNAs, a network of competing endogenous RNAs (ceRNAs) associated with the occurrence and progression of OS was created to show the microRNA (miRNAs) targets of m6A-related lncRNAs and the corresponding mRNA targets of these miRNAs. The biological functions of mRNAs in the ceRNA network were analyzed to better understand the underlying mechanisms.

Materials and Methods

Data Acquisition

We downloaded the level 3 RNA-seq expression data comprising lncRNAs and 21 m6 A regulators, namely METTL3, METTL14, LRPPRC, WTAP, KIAA1429, ZC3H13, CBLL1, RBM15, RBM15B, FTO, ALKBH5, YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, IGF2BP1, FMR1, ELAVL1, HNRNPC, and HNRNPA2B1, of 88 patients with OS from TCGA database (https://portal.gdc.cancer.gov). Subsequently, the latest clinical data of patients with OS were downloaded from the TARGET database (https://ocg.cancer.gov/programs/target). Finally, both RNA-seq expression data and valid clinical information of 85 patients were used as a training set for analysis in this study. The clinical information of patients with OS in TCGA is shown in Supplementary Table S1. For the validation set, the RNA-sequencing data and corresponding clinicopathological features of 37 patients with OS were downloaded from the GSE39055 database.

Construction and Evaluation of Risk Score

Common lncRNAs were first extracted from TCGA and GSE39055 datasets. Pearson correlation analysis was then performed to find m6A-related lncRNAs with cor > 0.3 and P < .05, and univariate Cox regression analysis was performed to filter genes significantly related to survival. Consequently, five m6A-related lncRNAs were identified and an m6A-related lncRNA prognostic signature was established using the multivariate Cox regression analysis. Finally, the risk score was calculated using the following formula: where Coef is the coefficient and Exp is the normalized expression value of each signature gene. The risk score system of 2 m6A-related lncRNAs was built in the TCGA-TARGET-OS training set and evaluated in the GSE39055 validation set. The patients were divided into high-risk and low-risk groups based on the median risk score. Univariate and multivariate Cox regression analyses were then performed to assess whether the risk score was an independent prognostic factor for the survival of patients with OS.

The Relationship Between Risk Score and Immune Cell Infiltration in OS

Single-sample gene set enrichment analysis (ssGSEA) was performed to analyze the enrichment score of 24 immune cell types in each sample: dendritic cells (DCs), activated DCs (aDCs), B cells, CD8 T cells, cytotoxic T cells, eosinophils, immature DCs (iDCs), macrophages, mast cells, neutrophils, natural killer (NK) cells, NK CD56 bright cells, NK CD56 dim cells, plasmacytoid DCs (pDCs), T cells, T helper cells, central memory T (Tcms) cells, effector memory T (Tems) cells, follicular helper T (TFH) cells, gamma delta T (Tgds) cells, T helper 1 (Th1) cells, T helper 2 (Th2) cells, T helper 17 (Th17) cells, and regulatory T (TRegs) cells. The relationship between risk score and immune cell infiltration was calculated by Pearson correlation.

Establishment of the Nomogram

A nomogram was built using multivariate Cox regression analysis, and calibration curves were created to show the concordance between the nomogram-predicted probability and observed overall survival of patients with OS. Decision curves were plotted to evaluate the clinical utility of the nomogram by estimating the net benefits at a range of threshold probabilities.

Construction of the ceRNA Network

First, the StarBase database was used to predict the target miRNAs of small nucleolar RNA host gene 7 (SNHG7) and small nucleolar RNA host gene 12 (SNHG12). Then, the corresponding target mRNAs of these miRNAs were selected from StarBase, and the lncRNA-miRNA and miRNA-mRNA regulatory relationships were integrated to construct the ceRNA network using Cytoscape.

Functional Analysis

Gene set enrichment analysis (GSEA) was performed to investigate the most common pathways in the high-risk or low-risk groups, and gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of the target mRNAs in the ceRNA network were performed using the clusterProfiler R package.

Statistical Analysis

The difference in overall survival between the high-risk and low-risk groups was assessed using the Kaplan–Meier method and log-rank test, and the predictive accuracy of the risk model was determined using the receiver operating characteristic (ROC) curves and area under the curve (AUC) values. Student's t-test was performed to compare the risk score values between subgroups according to the following clinicopathological features: age (<14 or >14 years), gender (male or female), and metastatic status (metastatic or non-metastatic). All data were analyzed using R software and results with a P value <.05 were considered as statistically significant.

Results

Identification of m6A-Related lncRNAs in Patients with OS

First, lncRNA data of patients with OS were downloaded from TCGA-TARGET-OS and GSE39055, and a total of 122 common lncRNAs were identified (Supplementary Table S2). Subsequently, data on the expression of 21 m6 A regulators were obtained from TCGA-TARGET-OS. Pearson correlation analysis was performed to identify m6A-related lncRNAs, the expression levels of which were correlated with one or more of the 21 m6 A regulators. Consequently, 59 m6 A-related lncRNAs were identified as significant in the TCGA dataset (Supplementary Table S3). LncRNAs with the highest correlation with each m6A regulator are shown in Figure 1. These significant m6A-related lncRNAs were further examined in the study.

Figure 1.

Heat map representing the correlation between lncRNAs and m6A-related genes.

Construction of Risk Score Signature Using m6A-Related lncRNAs with Prognostic Value

We further explored the role of these m6A-related lncRNAs in predicting the prognosis of patients with OS by performing univariate Cox regression analysis in TCGA training set. We found that five m6A-related lncRNAs were related to patient survival (P < .05) and all were risk factors (HR > 1) for OS (Figure 2A). Furthermore, multivariate Cox regression analysis was performed to increase the predictive robustness of the five m6A-related lncRNAs in TCGA training set. Two m6A-related lncRNAs were screened to construct the risk model, in which the risk score was calculated as 0.062 × expression value of SNHG7 + 0.115 × expression value of SNHG12 (Figure 2B). Based on the median value of the risk score, we divided the patients from TCGA training set into high-risk and low-risk groups. The expression values of SNHG7 and SNHG12 are represented in the heat map (Figure 2C). The distributions of risk score and survival status of patients with OS in TCGA training set are represented in Figure 2D. The results of the Kaplan–Meier analysis showed that the high-risk group had a worse clinical outcome (Figure 2E), and the ROC curve indicated that the two m6A-related lncRNA signatures had high accuracy in predicting the overall survival in TCGA training set with an AUC value of > 0.65 (Figure 2F). Additionally, we investigated the relationship between the risk score value and clinicopathological features of patients with OS in TCGA training set. Although female patients aged >14 years and with metastatic disease had high-risk score values, the differences were not statistically significant (Figure 2G).

Figure 2.

Screening for prognostic biomarkers and risk model construction. (A) Forest plot of univariate Cox regression analysis to screen for prognosis-related lncRNAs. (B) Forest plot of multivariate Cox regression analysis to screen for prognosis-related lncRNAs. (C) Heat map showing SNHG7 and SNHG12 expression levels in the training set. (D) Risk curve in the training set. (E) Kaplan–Meier curve in the training set. (F) ROC curve in TCGA training set. (G) Correlation between risk scores and clinical outcomes.

Validation of m6A-Related lncRNA Signature in the GEO Dataset

We further validated the risk model in the GSE39055 dataset. The patients were assigned into low- and high-risk groups according to the median value of the risk score, which was calculated using the same formula as mentioned previously. A heat map was generated to show the expression values of SNHG7 and SNHG12 in the low- and high-risk groups (Figure 3A). Similar to the results in TCGA training set, we found that the mortality status was mainly distributed in patients with higher risk score values (Figure 3B), and patients in the high-risk group had shorter survival time (Figure 3C). The ROC curve also demonstrated that the risk model had a strong efficiency for predicting the prognosis of patients with OS in the GSE39055 dataset (AUC = 0.671, Figure 3D).

Figure 3.

Validation of the risk model in the GEO dataset. (A) Heat map showing SNHG7 and SNHG12 expression levels in the validation set. (B) Risk curve in the validation set. (C) Kaplan–Meier curve in the validation set. (D) ROC curve in the GEO validation set.

Identification of Pathways Related to the Risk Model

To identify the pathways related to the risk model, we performed GSEA to investigate the potential biological processes and pathways in high- and low-risk groups. We found that the B-cell receptor signaling and NK cell pathways were enriched in patients in the low-risk group (Figure 4A), while RNA polymerase and cell cycle pathways were enriched in patients in the high-risk group (Figure 4B), suggesting that in OS, the m6A-related lncRNA risk signature may be related to immune system-related and cell proliferation. Thus, we calculated the infiltration levels of 24 immune cell types in each sample (Supplementary Table S4). The results of Pearson correlation showed that the risk score was significantly negatively correlated with iDCs (P value < .05, r = −0.43), macrophages (P value < .05, r = −0.4), neutrophils (P value < .05, r = −0.47), and NK cells (P value < .05, r = −0.41) (Figure 5).

Figure 4.

Identification of pathways related to the risk model. (A) GSEA analysis in low-risk group. (B) GSEA analysis in high-risk group.

Figure 5.

Relationship between risk score and immune cell infiltration in OS. ssGSEA was performed to analyze the enrichment score of 24 immune cell types in each sample. Pearson correlation between risk score and each immune cell type was calculated and visualized in a scatter plot.

Generation of Nomogram for Patients with OS

We then performed univariate and multivariate Cox regression analyses to explore whether the risk score and metastasis were independently associated with prognosis (Figure 6A). The univariate Cox regression analysis revealed that the risk score and metastasis were strongly associated with prognosis, and the multivariate Cox regression analysis further confirmed that the risk score was an independent prognostic factor (P < .001, Figure 6B). A nomogram for predicting the 1/3/5-year survival of patients with OS was created based on independent prognostic factors (risk score and metastasis) (Figure 6C). We found a good concordance between the observed and predicted 1- and 3-year overall survival rates, as shown in the calibration plot (Figure 6D). In addition, we plotted decision curves to assess the clinical utility of the nomogram, and found that it yielded more net benefit for predicting the 3- and 5-year survival rates (with the threshold probability > 0.2 and > 0.25, respectively) than the treat-all and treat-none strategies (Figure 6E, F).

Figure 6.

Construction of nomogram for patients with OS. (A) Forest plot of univariate Cox regression analysis to screen for factors associated with prognosis. (B) Forest plot of multivariate Cox regression analysis to screen for factors associated with prognosis (C) Nomogram constructed with risk scores and metastasis. (D) Calibration plot of nomogram. (E) Decision curves for predicting 3-year survival. (F) Decision curves for predicting 5-year survival.

Construction and Analysis of the ceRNA Network

To further elucidate how the m6A-related lncRNAs (SNHG7 and SNHG12) regulate mRNA expression via sponging miRNAs in OS, we searched for miRNAs interacting with SNHG7 and SNHG12, and their target mRNAs in StarBase and constructed the ceRNA network (Figure 7A). Moreover, we investigated the biological functions of these target mRNAs and found that these genes were enriched in ficolin-1-rich granule, ficolin-1-rich granule lumen (GO cellular component), DNA polymerase binding, nuclear receptor binding (GO molecular function), estrogen signaling pathway, and hepatocellular carcinoma (KEGG pathway), as shown in Figures 7B and 7C.

Figure 7.

Establishment and analysis of the ceRNA network. (A) Construction of ceRNA network based on SNHG7 and SNHG12. (B) GO analysis of mRNAs in the ceRNA network. (C) KEGG analysis of mRNAs in the ceRNA network.

Discussion

In recent years, m6A modification has been reported as a common internal modification of eukaryotic mRNAs and can regulate cell proliferation, invasion and apoptosis, thereby modulating the progression of OS. m6A-related methyltransferases, including METTL3, METTL14, and WTAP, have shown critical roles in OS cell growth and metastasis by modifying the expression of corresponding downstream genes.[25,32-35] YTHDF2, a key reader protein, was found to serve as the main factor involved in abnormal m6A modification of tripartite motif-containing 7, thus, promoting tumorigenesis and chemoresistance in OS. Moreover, current evidence indicates that m6A modification of lncRNAs is involved in the tumorigenesis of OS. However, it is still largely unknown how m6A modification acts in a lncRNA-dependent manner during OS development. Hence, a comprehensive bioinformatics analysis was performed to identify m6A-related lncRNAs and determine their prognostic value and potential underlying mechanisms in OS.

Several studies have suggested that dysregulated m6A-related lncRNAs are strongly associated with tumorigenesis and clinical outcomes of patients with tumors. In the current study, data of 122 patients with OS from TCGA and GSE39055 datasets were analyzed to investigate the significance of m6A-related lncRNAs, and the pooled results showed that the expression of 59 lncRNAs was significantly correlated with the expression of one or more m6A-related regulators, suggesting that these lncRNAs may be modified in an m6A-dependent manner. Furthermore, we explored the role of these m6A-related lncRNAs in the prognosis of OS by performing univariate Cox regression analysis; five m6A-related lncRNAs, namely SNHG1, SNHG6, SNHG7, SNHG8, and SNHG12, were related to the overall survival of patients with OS.

Similarly, previous studies have confirmed that these lncRNAs are involved in the oncogenesis and metastasis of OS. SNHG1 was upregulated in OS tissues and promoted cell proliferation, invasion, and migration by serving as a ceRNA.[38-40] Highly expressed SNHG6 was positively associated with poor overall survival of patients with OS and could enhance OS cell proliferation by regulating the miR-26a-5p/ULK1 axis and expression of p21 and Kruppel-like factor 2 (KLF2).[41,42] It has been reported that SNHG7 was aberrantly expressed in OS tissues and elevated levels of SNHG7 could facilitate tumor growth and epithelial-mesenchymal transition (EMT) by regulating p53 expression and miR-34a levels.[43,44] Recent studies have confirmed that SNHG8 could stimulate the growth and migration of OS cells by acting as a ceRNA to sponge miR-876-5p and miR-542-3p, thereby suggesting that SNHG8 may be a novel therapeutic target for the treatment of OS.[45,46] SNHG12 was also found to be overexpressed in OS tissues and silencing SNHG12 may attenuate Notch2- and insulin-like growth factor 1 receptor (IGF1R)-induced tumorigenesis and metastasis via regulating miR-195-5p.[47,48] Multiple studies have indicated that upregulation of these five lncRNAs are believed to promote the proliferation, migration, and inhibit apoptosis of tumr cells, thereby enhancing chemoresistance.[49-53] Their dysregulation is positively correlated with poor clinical outcomes of patients with OS,[41,54] which is consistent with our findings. Furthermore, m6A-related lncRNAs may influence the expression of m6A regulators via cis-or trans-acting regulation, or by functioning as ceRNAs, thereby forming an integrated regulatory network. Hence, additional studies are required to investigate the function of m6A-related lncRNAs in the initiation and development of OS, to identify novel prognostic biomarkers and therapeutic targets.

Furthermore, multivariate Cox regression analysis was performed to identify m6A-related lncRNAs that could predict the prognosis of OS more accurately and reliably. SNHG7 and SNHG12 were found to have such prognostic value and were used to create a risk model which was confirmed to be effective in predicting the clinical outcomes of patients with OS. Notably, patients enrolled from TCGA were divided into high-risk and low-risk groups based on the median value of the risk score. The expression levels of both SNHG7 and SNHG12 were elevated in the high-risk group, suggesting that these two m6A-related lncRNAs could be risk factors for OS. Moreover, the relationship between the risk score and the clinicopathological features of patients with OS was explored. The results showed that the risk score values were higher in female patients aged >14 years with metastatic disease when compared to each respective control group. However, no statistical significance was observed. These nonsignificant findings may be attributed to small sample sizes. Nevertheless, the univariate and multivariate regression analyses confirmed that the risk score based on SNHG7 and SNHG12 was an independent prognostic factor. Furthermore, a nomogram for predicting the 1/3/5-year survival of patients with OS was created based on independent prognostic factors. Both the calibration plot and decision curves showed that this nomogram indicated a good prognostic value.

To date, few studies have focused on the molecular mechanisms of these two lncRNAs in the progression of OS,[43,44,47,48,54,56] and further research regarding their role in the pathophysiological process of OS is needed. Nevertheless, there are no previous studies reporting the important roles that these lncRNAs might play in OS development via cross-talk with their downstream targets in an m6A-dependent manner. Hence, using bioinformatics analysis, a ceRNA network was created to search for potential miRNA and mRNA targets. Multiple OS-related m6A regulator-lncRNA-miRNA-mRNA networks were identified. Moreover, these genes were enriched in several cancer-related pathways including estrogen signaling pathway, hepatocellular carcinoma, gastric cancer, and necroptosis. Furthermore, the risk model based on SNHG7 and SNHG12 was enriched in immune system-related and cell proliferation pathways. The ssGSEA results further demonstrated the close relationship between risk score and immune cell infiltration. Considering these findings, these two m6A-related lncRNAs and their downstream target molecules may be positively associated with the tumorigenesis and prognosis of OS.

Nevertheless, there are some limitations in this study. First, transcriptome data and clinicopathological features were obtained from TCGA and GEO databases. However, there were only 85 patients in the training set and 37 in the validation set. The prognostic value of these m6A-related lncRNAs should be validated in more independent cohorts. Second, the underlying mechanisms of these m6A-related lncRNAs were evaluated only by bioinformatics-based prediction. These findings should be further corroborated through in vitro and in vivo experiments. Third, although higher risk score values were observed in female patients aged >14 years with metastatic disease, the differences were not significant. Studies with larger sample sizes are needed to further validate these results.

Conclusion

In summary, to our knowledge, this study is the first to comprehensively screen for m6A-related lncRNAs that influence the prognosis of patients with OS. Our results showed that SNHG7 and SNHG12 may be valuable prognostic biomarkers for OS, and a risk model that could accurately predict the survival of patients with OS was established based on these two m6A-related lncRNAs. A ceRNA network was created to gain further understanding of these two lncRNAs targeting downstream molecules in the pathogenesis of OS. These findings may provide novel insights into the molecular mechanisms of OS by giving a new perspective on m6A modification.

Click here for additional data file.

Supplemental material, sj-xlsx-1-tct-10.1177_15330338221085354 for N6-methyladenosine Modification-Related Long Non-Coding RNAs are Potential Biomarkers for Predicting the Prognosis of Patients With Osteosarcoma by Kun Yang, Fengyan Wang, Ke Li, Guoxuan Peng, Hua Yang, Hong Xu, Yang Xiang and Hong Sun in Technology in Cancer Research & Treatment

Click here for additional data file.

Supplemental material, sj-xlsx-2-tct-10.1177_15330338221085354 for N6-methyladenosine Modification-Related Long Non-Coding RNAs are Potential Biomarkers for Predicting the Prognosis of Patients With Osteosarcoma by Kun Yang, Fengyan Wang, Ke Li, Guoxuan Peng, Hua Yang, Hong Xu, Yang Xiang and Hong Sun in Technology in Cancer Research & Treatment

Click here for additional data file.

Supplemental material, sj-xlsx-3-tct-10.1177_15330338221085354 for N6-methyladenosine Modification-Related Long Non-Coding RNAs are Potential Biomarkers for Predicting the Prognosis of Patients With Osteosarcoma by Kun Yang, Fengyan Wang, Ke Li, Guoxuan Peng, Hua Yang, Hong Xu, Yang Xiang and Hong Sun in Technology in Cancer Research & Treatment

Click here for additional data file.

Supplemental material, sj-xlsx-4-tct-10.1177_15330338221085354 for N6-methyladenosine Modification-Related Long Non-Coding RNAs are Potential Biomarkers for Predicting the Prognosis of Patients With Osteosarcoma by Kun Yang, Fengyan Wang, Ke Li, Guoxuan Peng, Hua Yang, Hong Xu, Yang Xiang and Hong Sun in Technology in Cancer Research & Treatment