Experimental Study of Somatic Variants of Osteosarcoma by Whole-Exome Sequencing.

BACKGROUND This study aimed to investigate the role of gene mutation site distribution, biological function, pathway enrichment, and gene association analysis in the occurrence, development, and migration of osteosarcoma. MATERIAL AND METHODS Somatic mutation screening was performed using the whole-exome sequencing of osteosarcoma samples, and the distribution of mutations was demonstrated by Circos diagrams. Metascape was used to analyze the GO and KEGG signal pathway enrichment of the genes harboring protein coding alterations, and GeneMANIA was used to analyze the interaction of mutated genes. RESULTS The results showed that the protein coding alterations were found throughout the whole genome in 3 osteosarcoma samples. A large number of identical or related biological processes and pathways were found in osteosarcoma samples. The GeneMANIA analysis of the 10 mutations shared by 3 samples showed that the target gene minichromosome maintenance complex component 4 (MCM4) and 3 lateral genes were most functional, and were all related to DNA replication. The analysis of GO and KEGG signal pathway enrichment showed that the mutated genes were involved mainly in tumor-related metabolic pathways. Three mutated genes were involved in the cell process, and 2 mutated genes were involved in the metabolic process. Known driver gene mutations were also observed in the samples. CONCLUSIONS The gene analysis confirmed that patients with osteosarcoma had a wide range of common gene mutations related to each other, which are involved in tumor-related metabolic pathways. These findings provide a basis for further gene-targeted therapy and pathway research.

Background

Osteosarcoma is a primary malignant tumor of the bone that occurs in the metaphysis of long bones in children and adolescents. The treatment of osteosarcoma has remained almost unchanged for nearly 30 years. With the standardized implementation of neoadjuvant chemotherapy, the 5-year survival rate of patients without metastasis has increased to 60–70%. However, about 15–20% of patients have different degrees of metastatic organ diseases at the time of definite diagnosis, and the 5-year survival rate of patients with metastasis is still less than 30% [1-3]. It is clinically urgent to find a safe and effective treatment method to detect the disease at an early stage and improve the survival rate of patients with metastasis. Genomic studies and targeted treatments have become important areas of research in recent years.

At present, the mechanisms underlying the involvement of genes, pathways, and biological processes in the occurrence, development, and migration of osteosarcoma are unclear. Moreover, no effective biomarkers are available for the early detection and treatment of osteosarcoma. This study started with the genomic study of osteosarcoma. In this stage, whole-exome sequencing and data analysis were performed for the pathological tissue and peripheral blood samples of 3 patients with osteosarcoma. This study aimed to find genes and biological processes related to osteosarcoma, providing new directions for the further study of the occurrence, development, metastasis, and prognosis of osteosarcoma.

Material and Methods

Sample selection and genomic DNA extraction

Pathological tissue and peripheral blood samples were obtained from 3 patients with typical osteosarcoma from the Department of Bone and Soft Tissue Oncology, Affiliated Cancer Hospital of Zhengzhou University. Genomic DNA was extracted for whole-exome capture and sequencing. After genomic DNA samples were randomly broken into fragments of 180–280 bp, end repair, addition of dA, and connection of the connector were performed to construct a library. After exon capture was completed using an Agilent SureSelect Human All exon V6 kit, the library was quantified and then sequenced on the Illumina HiSeq platform. The sequencing depth of the genomic library of the tumor tissues was 500×, and the sequencing depth of the genomic library of the blood samples was 100×. The original sequencing results were compared with the human reference genome Build 37 using a Burrows–Wheeler Aligner. The recognition of variants was implemented using SAMtools, annotations were completed using ANNOVAR, v2013Aug23, and recognition of somatic variant was performed using MuTect.

Data processing

The variants in the 3 patients with osteosarcoma were counted based on single-nucleotide variation (SNV), insertion or deletion of small-fragment genes (Indel), and changes in CNV. The distribution of mutant loci was displayed using a Circos plot. GO analysis and KEGG signal pathway enrichment analysis were performed using the Metascape software. Interaction of common mutant genes was analyzed using GeneMANIA.

Selection of genetic mutations

Somatic oncogenes and tumor suppressor genes were searched by checking the PubMed database. The variant loci of somatic cells specific to osteosarcoma samples were searched by the sequence alignment of pathological tissue and peripheral blood samples. Meaningful mutant loci and common mutant genes were screened according to the information on mutant loci.

Results

Basic clinical information and pathological sections of the patients

Table 1 and Figure 1.

Table 1

Basic clinical information of 3 patients with osteosarcoma.

Sample number	Sex	Age	Part	Pathologic figure
T18011106001	Male	20	Distal left femur	Figure 1A
T18011706001	Male	34	Proximal Right humerus	Figure 1B
T18022706001	Female	15	Distal right radius	Figure 1C

Figure 1

(A–C) Pathology of HE staining in 3 patients with osteosarcoma.

Statistics of mutant types and distribution of mutant loci with changed coding

The comprehensive analysis of the results in Table 2 and Figure 2A and 2B showed that meaningful mutant loci were not concentrated in a single or a minority region, but they were located in the whole-exon region. Additionally, Figure 2A and 2B show that a small number of mutant genes were present in the region of sex chromosomes, while a large number of mutant genes were distributed on chromosomes 1, 8, 9, 18, 19, and 22, which might be associated with the occurrence, development, and migration of osteosarcoma. However, this requires further investigation.

Table 2

Somatic variation measured in 3 osteosarcoma samples.

Simple name	Mutation type
	SNV	CNV	Indel
T18011106001	510	68	25
T18011706001	712	60	19
T18022706001	2621	203	27
Total	3843	331	71

Figure 2

Figure A and B represent the distribution diagram of the mutant loci of somatic cells and the distribution diagram of meaningful mutant loci excluding introns plotted by Circos, which were hg19 whole genome T18011106001, T18011706001, and T18022706001 from the outer circle to the inner circle. Figure C is a gene correlation diagram among the 3 samples with functional annotation and non-intron regions screened using Metascape software.

Analysis of biological functions and pathway enrichment

Figure 2C shows that the blue lines were distributed widely and densely, suggesting a large number of functionally related biological processes or relevant pathways in the 3 patients. The mutant genes were the most abundant in sample T18022706001. The distribution of purple lines indicated numerous genetic overlaps between sample T18011106001 and sample T18022706001, followed by genetic overlaps between sample T18011706001 and sample T18022706001; the least genetic overlaps were between sample T18011106001 and sample T18011706001. A large number of genetic overlaps suggested numerous identical or related biological processes and pathways between samples. Therefore, the occurrence, development, migration, and pathogenesis of osteosarcoma could be explored by analyzing the common biological processes or related pathways in the 3 patients.

Figure 3 shows that sample T18022706001 had more enrichment, mainly involving the following biological processes or pathways: cell cycle, cellular responses to external stimuli, apoptotic signaling pathway, adaptive immune system, and signaling by interleukins. Sample T18011106001 was enriched mainly in the positive regulation of defense response and the response of cells to external stimuli. Sample T18011706001 was enriched mainly in the positive regulation of defense response.

Figure 3

Cluster analysis was performed for the genetic data of the 3 samples using Metascape software. The convergence and difference among the 3 samples were displayed by biological functions and pathway enrichment.

In Table 3, annotations simultaneously enriched to the gene sequence of the 3 samples are counted separately. The count column shows the number and percentage of genes enriched to the annotation in this gene sequence.

Table 3

The statistical results of GO annotation were simultaneously enriched to 3 samples.

ID	Description	Count	%
R-HSA-1280218	Adaptive immune system	142	7.23
R-HSA-1640170	Cell cycle	121	6.16
R-HSA-449147	Signaling by interleukins	119	6.06
GO: 0031349	positive regulation of defense response	104	5.29
R-HSA-556833	Metabolism of lipids	129	6.56
R-HSA-8953897	Cellular response to external stimuli	108	5.50
R-HSA-3700989	Transcriptional Regulation by TP53	73	3.72
GO: 0045859	Regulation of protein kinase activity	129	6.56
GO: 001903827	Regulation of cellular protein localization	94	4.78

The comparison in Figure 4A–4C revealed which biological processes or related pathways and to what extent genes in the 3 samples were enriched. Figure 4A shows that 4 biological processes – apoptotic signaling pathway, positive regulation of proteolysis, cell cycle, and adaptive immune system – were the pivotal links in the network diagram. Figure 4B shows that the degree of enrichment of the positive regulation of the defense response was the maximum. Figure 4C shows that genes of sample T1802276001 were involved in almost all the biological processes. The genetic distribution of sample T18011106001 was mainly in the positive regulation of defense response, positive regulation of proteolysis, cell cycle, and regulation of protein kinase activity. Moreover, the degree of enrichment was high. Genes of sample T18011706001 were distributed mainly in the positive regulation of defense response and regulation of cellular protein localization.

Figure 4

Biological function processes or networks of related pathways enriched by each sample were constructed using the Metascape software, and the results are shown in Figure 4. Figure A shows that biological processes or related pathways enriched in the 3 samples were closely correlated. In Figure B, colors were applied based on the enrichment degree of biological processes or related pathways. Based on the shade of colors, it was seen that function-related classes had a large number of enrichment pathways. In Figure C, colors were applied based on different genetic sequences, which displayed the distribution of genes of the 3 samples in various biological processes or related pathways.

Association analysis among genes

The common mutant genes with changed coding in the 3 patients with osteosarcoma were screened, and 11 genes were obtained: CD8b molecule (CD8B), inner membrane mitochondrial protein (IMMT), minichromosome maintenance complex component 4 (MCM4), protein kinase, DNA-activated, catalytic subunit (PRKDC), heterogeneous nuclear ribonucleoprotein K (HNRNPK), rho-associated coiled-coil-containing protein kinase 1 pseudogene 1 (ROCK1P1), eukaryotic translation elongation factor 1 alpha 1 (EEF1A1), actin beta (ACTB), protein phosphatase 1 catalytic subunit beta (PPP1CB), tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein zeta (YWHAZ), and dedicator of cytokinesis 8 (DOCK8). Among them, ROCK1P1 was a pseudogene, and its interaction with other genes was not considered. Figure 5A depicts the output after inputting the 10 genes into GeneMANIA, clearly showing that the target genes DOCK8 and CD8B were not correlated with other genes. Functions involving the target gene MCM4 and the lateral genes replication factor C subunit 2 (RFC2), minichromosome maintenance complex component 2 (MCM2), and minichromosome maintenance complex component 6 (MCM6) were the most abundant, and these genes were all correlated with DNA replication.

Figure 5

The interaction network diagram of genes was obtained after inputting the screened 10 genes into GeneMANIA. In Figure A, functions involved by target genes and lateral genes are marked by different colors. The functions related to the gene were obtained by the color distribution of this gene. In Figure B–D, the middle circles represent target genes, and the outer circles represent the outer genes. Figure B displays the network of physical interaction between genes, which was constructed based on whether genetic products had interactions. Figure C shows the co-expression network between genes, which was constructed based on the expression level of genes. Figure D shows the predicted functional relationship between genes, usually referring to the interaction between gene products (proteins).

Figure 5B shows physical interactions between the YWHAZ gene and the 6 target genes, as well as between the YWHAZ gene and the 11 lateral genes. The EEF1A1 gene physically interacted with the 5 target genes and the 5 lateral genes. Thus, it could be concluded that the EEF1A1 and YWHAZ genes are important in the target gene network diagram. According to the thickness of the lines in the diagram, it was concluded that the physical interaction between the PPP1CB gene and the mevalonate kinase (MVK) gene was the largest. Physical interactions between the PPP1CB gene and the density-regulated re-initiation and release factor (DENR) gene and between the EEF1A1 gene and the phosphofructokinase platelet (PFKP) gene were the second largest. Figure 5C shows that the YWHAZ, PRKDC, and MCM4 genes were co-expressed with more genes. More co-expression between target genes, as well as co-expression of YWHAZ and EEF1A1 genes with other target genes, was observed, which was consistent with the results in Figure 5A. Figure 5D shows that strong functional interactions may exist between the PPP1CB gene and the neurofibromin 2 (NF2) gene and between the EEF1A1 gene and the eukaryotic translation elongation factor 1 beta 2 (EEF1B2) gene. In summary, it was reasonable to believe that these genes are functionally correlated, which lays the foundation for further exploration of common biological processes or related pathways.

WEGO: Statistical analysis of GO annotations

A total of 11 common mutant genes had changed the coding in the 3 samples. Among them, ROCK1P1 was a pseudogene without biological functions. According to the remaining 10 genes, including the known results of GO annotations of the 7 genes (CD8B, IMMT, PRKDC, HNRNPK, EEF1A1, ACTB, and YWHAZ) and no GO annotation results of the 3 genes (MCM4, PPP1CB, and DOCK8), Figure 6 was obtained using WEGO. It was concluded that 3 common mutant genes were involved in the cell process, and 2 common mutant genes were involved in the metabolic process. Changes in these biological processes might lead to the production and proliferation of cancer cells.

Figure 6

Three aspects of biological functions, including cell component, molecular functions of the genes, and biological processes are displayed by a histogram using WEGO.

Driver genes

The known driver genes (Table 4) were screened in the 3 patients, and 3 known driver genes were obtained in sample T18011106001 [MYC proto-oncogene, BHLH transcription factor (MYC), patched 1 (PTCH1), and semaphorin 4D (SEMA4D)]. Moreover, 2 candidate driver genes were obtained von Hippel-Lindau tumor suppressor (VHL) and catenin delta 1 (CTNND1). Two known driver genes were present in sample T18011706001 [tumor protein P53 (TP53) and GRB2-associated binding protein 1 (GAB1)]. There were driver gene mutations in sample T18022706001, including 4 known driver genes: MYC, RB transcriptional corepressor 1 (RB1), TP53, and NF2. In addition, 9 candidate driver genes were present: platelet-derived growth factor receptor alpha (PDGFRA), disks large MAGUK scaffold protein 2 (DLG2), cyclin E1 (CCNE1), E1A-binding protein P300 (EP300), epidermal growth factor receptor (EGFR), MutY DNA glycosylase(MUTYH), NKD inhibitor of Wnt signaling pathway 2(NKD2), GNAS complex locus (GNAS), and BLM RecQ-like helicase (BLM). The candidate driver genes were involved in cell fate, cell survival, or genomic stability, which might be associated with the occurrence, development, and migration of osteosarcoma.

Table 4

Known driver genes of tumors.

Gene name	Description
TP53	Tumor suppressor gene that can cause apoptosis of cancer cells
WIF1	Tumor suppressor gene that is involved in the Wnt signaling pathway
NF2	Tumor suppressor gene that can maintain cytoskeletal stability
BRCA1/2	Tumor suppressor gene that can regulate cell proliferation and DNA damage repair
PRKAR1A	Tumor suppressor gene that encodes protein kinase A regulatory subunit, and is involved in the cAMP signaling pathway
APC	Tumor suppressor gene that regulates cell proliferation and migration and maintains chromosome stability
PTCH1	Tumor suppressor gene involved in the hedgehog signaling pathway
RB1	Tumor suppressor gene that can promote cells to the classification stage, leading to cell proliferation
PTEN	Tumor suppressor gene that can inhibit the growth and migration of tumor cells
BAP1	Tumor suppressor gene that regulates cell cycle
RUNX3	Tumor suppressor gene that regulates the growth and proliferation of gastric epithelial cells
CDKN2A/B	Tumor suppressor gene involved in regulating the cell cycle
P16	Tumor suppressor gene involved in regulating the cell cycle
GRB10	Tumor suppressor gene involved in muscle growth, regeneration, and repair
ZNF217	Oncogene, transcriptional regulator
ZNF592	Oncogene involved in the regulation of cerebellar development
GAB1	Oncogene involved in cell growth, transformation, and apoptosis
FOS	Oncogene that regulates cell proliferation, differentiation, and transformation
Notch1-4	Oncogene that mediates T cell development and differentiation
MYC	Oncogene involved in cell proliferation and apoptosis
JUN	Oncogene related to cell proliferation
TWIST	Oncogene involved in embryonic development that can promote the invasion and metastasis of tumor cells
RET	Oncogene involved in cell growth, differentiation, and migration
SEMA4D	Oncogene involved in PI3K-Akt-mTOR and MAPK signaling pathways
SEMA6D	Oncogene involved in PI3K-Akt-mTOR and MAPK signaling pathways
CAPRIN1	Oncogene that can promote metastasis and anti-apoptosis and participate in Akt and ERK 1/2 signaling pathways
AKT2	Oncogene that regulates metabolism, cell proliferation, and growth

Discussion

High-throughput sequencing in the 3 patients showed a large number of genetic mutations and overlaps. Using various statistical methods, we proved that genetic mutations in the 3 patients were involved in the biological process and pathways, and these genes were correlated. These biological changes and correlation analyses of genes were closely correlated with the characteristics of osteosarcoma, such as the occurrence, development, migration, and pathogenesis.

Pediatric patients with osteosarcoma reported by Chen et al. had mainly SNV, regardless of metastasis [4]. Mutations in the 3 patients were also SNV, which might be related to tumor heterogeneity. Moreover, in the study by Chen et al., a few somatic chromosomal diseases were also involved, and the distribution of mutant genes in diagnostic osteosarcoma and metastatic osteosarcoma was slightly different. However, most of these genes were concentrated in autosomes such as chromosomes 1, 8, 18, and 19, while a few mutant genes were distributed on sex chromosomes, which was basically consistent with the findings of the present study.

The statistical results in Table 3 show that the biological processes involved in mutant genes in the 3 patients were the adaptive immune system, transcriptional regulation by the TP53 gene, metabolism of lipids, and cell cycle. Changes in these biological processes can result in a cell cycle disorder, leading to the indefinite proliferation of cells, which is correlated with the infinite reproductive property of osteosarcoma cells. Among them, the transcriptional regulation by the TP53 gene and cell cycle were proven to be correlated with the occurrence and development of osteosarcoma. In addition, interleukin-8 can promote tumor invasion and formation of new blood vessels [5,6]. Hence, signaling by interleukins may be related to the occurrence and development of osteosarcoma. Regulation of protein kinase activity plays important roles in cell growth, proliferation, and differentiation. For example, receptor tyrosine kinase (RTK) is involved in many biological processes such as cell growth, tumor metastasis, and malignant transformation of tumors [7]. The biological process of signaling by regulation of protein kinase activity co-enriched in the 3 samples may be related to the characteristics of osteosarcoma, such as metastasis. Figure 4A–4C suggests that the biological processes or pathways in which the 3 samples were involved with a high degree of enrichment included mainly positive regulation of defense response and cellular responses to external stimuli, which might be related to the infinite replication property of osteosarcoma cells.

The YWHAZ gene is overexpressed in tumors such as lung cancer, liver cancer, and breast cancer. A large number of studies have shown that the YWHAZ gene can promote the proliferation, migration, and invasion of cancer cells [8-10]. However, Zeng et al. showed that overexpression of the YWHAZ gene did not affect the positive regulatory effect of the tripartite motif–containing 21 (TRIM21) gene on the proliferation of osteosarcoma cells [11]. The role of the YWHAZ gene in osteosarcoma remains controversial. Blanch et al. found that another important gene, EEF1A1, is anti-apoptotic. Overexpression of the EEF1A1 gene can specifically inhibit the TP53 gene [12], while TP53 is a well-known tumor suppressor gene. A study in 2018 proved that the EEF1A1 gene is involved in the occurrence of osteosarcoma through biological processes by constructing a PPI network and by GO analysis [13]. Figure 5B–5D shows that the EEF1A1 and YWHAZ genes were important genes in the network diagram of target genes. Figure 5D shows a strong functional correlation between the PPP1CB gene and the NF2 gene. The NF2 gene is one of the known driver genes of tumors. Hence, the PPP1CB gene may to some extent be associated with the occurrence and development of osteosarcoma. A study demonstrated that PPP1CB is a regulator of endothelial cell migration, participating in the process of angiogenesis [14], thus creating conditions for tumorigenesis. Hassan et al. analyzed the differences in the expression of EEF elongation factors in cancer cells, suggesting that the EEF1A1 and EEF1B2 genes play a role in tumorigenesis and affect survival in a cancer-specific manner [15]. The eEF1 complex plays a crucial role in the de novo synthesis of proteins. EEF1A constitutes a central functional component, and the EEF1B2 gene is one of the remaining subunits of the complex [16]. It is reasonable to believe that a strong functional correlation exists between the EEF1A1 gene and the EEF1B2 gene. The PRKDC gene plays a role in the expression of osteosarcoma cell line MG63 by participating in DNA damage repair [17]. MCM4 is a gene involved in cell division, affecting the regulation of the osteosarcoma cell cycle [18]. These are all important genes, resulting in osteosarcoma formation. Their functional correlation should be further investigated.

Previous studies on whole-exome sequencing demonstrated that the TP53 and RB1 genes played a major role in recurrent changes in osteosarcoma (80–90% and 10–39%, respectively) [4-19]. About 50% of patients with osteosarcoma have inactivating mutations in both the TP53 and RB1 genes. The synergistic effect of TP53/Rb1 mutations shortens the average latency of malignant osteosarcoma and also causes the migration of osteosarcoma cells [20]. In the present study, 2 samples had the MYC gene, while the other sample had the TP53 gene.

Conclusions

The occurrence and development of osteosarcoma is the result of joint regulation by multiple genes, multiple steps, and multiple factors. In the present study, the mutation types of osteosarcoma and distribution of mutant loci, analysis of biological functions and pathway enrichment, genetic association analysis, GO annotation, and driver genes were explored in 3 patients with osteosarcoma using whole-exome sequencing technology. Widespread gene mutations were found in patients, which were involved in tumor-related metabolic pathways and provide a basis for further gene-targeted therapy and pathway research. However, our study was limited to the gene level, and further studies are needed to investigate the downstream events. In addition, this study used whole-exome sequencing, which has a certain limitation in investigating the structural mutations of the genome. Further, osteosarcoma is a rare disease with a very low incidence rate and our study had a small sample size. More samples need to be used for improved osteosarcoma genome research.