Osteosarcoma is the most common primary bone malignancy, and the lung is the most frequent site of metastasis. The limited understanding of the tumoral heterogeneity and evolutionary process of genomic alterations in pulmonary metastatic osteosarcoma impedes development of novel therapeutic strategies. Here we systematically illustrate the genomic disparities between primary tumors and corresponding pulmonary metastatic tumors by multiregional whole-exome and whole-genome sequencing in 86 tumor regions from 10 patients with osteosarcoma. Metastatic tumors exhibited a significantly higher mutational burden and genomic instability compared with primary tumors, possibly due to accumulation of mutations caused by a greater number of alterations in DNA damage response genes in metastatic tumors. Integrated analysis of the architecture and relationships of subclones revealed a dynamic mutational process and diverse dissemination patterns of osteosarcoma during pulmonary metastasis (6/10 with linear and 4/10 with parallel evolutionary patterns). All patients demonstrated more significant intertumoral rather than intratumoral heterogeneity between primary tumors and metastatic tumors. Mutated genes were enriched in the PI3K–Akt pathway at both the early and late stages of tumor evolution and in the MAPK pathway at the metastatic stage. Conversely, metastatic tumors showed improved immunogenicity, including higher neoantigen load, elevated PD-L1 expression, and tumor-infiltrating lymphocytes than the corresponding primary tumors. Our study is the first to report the dynamic evolutionary process and temporospatial tumor heterogeneity of pulmonary metastatic osteosarcoma, providing new insights for diagnosis and potential therapeutic strategies for pulmonary metastasis.

Significance:

High-throughput sequencing of primary and metastatic osteosarcoma provides new insights into the diagnosis of and potential clinical therapeutic strategies for pulmonary metastasis.

Osteosarcoma is the most common primary bone malignancy affecting the health of children, adolescents, and young adults (1). Despite the use of aggressive radical surgery and intensive (neo-)adjuvant chemotherapy, osteosarcoma relapse recurs in 30%–40% of patients (2). The most frequently relapsed organ is the lung, which accounts for approximately 85% of osteosarcoma metastases (3, 4). For patients with unresectable pulmonary metastases (metastatic tumors; MT), the prognosis is generally poor, with long-term survival rates of less than 20% (3, 4), even with delivery of active salvage chemotherapy and multiple targeted agents, such as FDA-approved pazopanib. The current status of osteosarcoma underlines the urgency and importance of understanding both the clinical biological behavior and genetic profiles of metastatic tumors, which is crucial for genetic diagnosis and disease management. Ultimately, this information can provide new insights into therapies to improve the prognosis of pulmonary metastatic osteosarcoma.

In recent years, several high-throughput genomic studies have characterized the genomes of patients with osteosarcoma with complex and a few consistent genomic alterations (5–9), unlike some sarcomas typified by specific chromosome translocations (10). These studies have revealed significant mutation events involved in the osteosarcoma etiology, including mutations in TP53, RB1, RECQL4, ATRX, DLG23, RUNX2, WRN, PTEN, and BRCA, the PI3K/mTOR pathway and IGF signaling–related genes (5–9). However, these potential osteosarcoma driver genes have failed to provide critical insight that can help improve osteosarcoma patient care, implying that targeting specific molecular aberrations may not be a rational strategy for the treatment of osteosarcoma.

On the other hand, studies examining the genomic traits of pulmonary metastatic osteosarcoma and those comparing primary tumors with metastatic tumors remain limited. Cancer metastasis is believed to be the result of survival adaptation of select tumor subclones to the metastatic organism microenvironment, constituting a major challenge for anticancer therapy and serving as the cause of malignant tumor-related deaths (11–12). Genetically, intratumoral heterogeneity (IaTH) due to the stochastic process of gene mutations and evolutionary progression determines the complexity of metastatic disease (13). Massive parallel sequencing of primary tumors and matched metastases has been performed for several cancers (14–17). However, a clear understanding of the ITH of osteosarcoma and the evolutionary history of pulmonary metastasis remains extremely limited. One important reason for this lack of understanding is that the genomic characteristics of osteosarcoma metastatic tumors have not yet been revealed, especially the molecular dynamics during metastasis, although the mutation landscape of primary osteosarcoma has been portrayed (6–9), which has impeded our understanding of the molecular and genetic mechanisms underlying the biological and clinical features to some extent, as well as the evolutionary trajectories of osteosarcoma metastatic tumors.

In this study, we performed comprehensive multiregional whole-exome sequencing (WES) on paired primary and pulmonary metastatic osteosarcomas and aimed to reveal the dynamic mutational evolutionary process during osteosarcoma pulmonary metastasis and the disparity of genomic alterations between primary tumors and metastatic tumors, which may shed light on novel therapeutic strategies and contribute to decisions in personalized treatment for patients with pulmonary metastatic osteosarcoma.

Study design

The initial cohort for screening included 321 patients with osteosarcoma diagnosed between September 2010 and September 2013 and 128 patients who developed pulmonary metastasis during their follow-up, 27 of whom received both primary and metastatic surgery at participating hospitals (Beijing Ji Shui Tan Hospital; National Cancer Center/Cancer Hospital, Chinese Academy of Medical Sciences & Peking Union Medical College; and Peking University Cancer Hospital, Beijing, China). Finally, the patients provided matched primary osteosarcoma and pulmonary metastasis samples with tumor cells occupying ≥50% of the content of the tumor tissue, and morphologically normal paracancerous tissues were also included. According to tumor size, 1–7 spatially separated tumor regions were obtained by microdissection. To detect the genomic aberrations and tumor heterogeneity, all specimens (n = 86) were utilized for WES (Supplementary Tables S1). A series of somatic mutations and copy number variants (CNV) were identified by bioinformatics analysis and were used to construct the phylogenetic tree and trace the evolutionary process during pulmonary metastatic osteosarcoma. Next, the genomic disparity between primary tumors and metastatic tumors, including single-nucleotide variants (SNV)/Indels, CNVs, and structural variations (SV), was described. Finally, the immunogenicity of each tumor region was evaluated, including the neoantigen burden, PD-L1 expression, and CNVs of immune-regulatory genes. One patient with PD-L1 expression and a high neoantigen burden was evaluated for response to pembrolizumab treatment. This study was approved by the ethics committees of all participating centers (201411-07, 2015KT13, and NCC2016G-074) and conducted according to CIOMS, and all enrolled patients provided written informed consent.

Tumor specimen processing

For the tumor specimens of the enrolled patients, we first cut off a piece of tumor tissue with a thickness of approximately 5 mm from the frozen tissue used to create formalin-fixed, paraffin-embedded (FFPE) samples. Then, an expert pathologist reviewed the hematoxylin and eosin slides of each tumor to ensure that the tumor cell contents were greater than 50%. Next, we obtained 10 cases with high-quality tumor tissues, including one case (OS07) with paired primary and metastatic FFPE samples and 9 cases with paired freshly frozen primary and metastatic samples. For the 9 paired frozen samples, 1–7 spatially separated tumor regions were obtained by microdissection according to the size of the tumor tissue. Each dissected region was approximately 5 × 2 × 2 mm, with at least 2 mm separating the edge of each region. Finally, DNA was extracted using the Qiagen Tissue kit according to the manufacturer's instructions, and the quality and integrity of the DNA were examined by agarose gel electrophoresis.

Genomic sequencing

For WES, a total of 0.6 μg of genomic DNA was used for sample preparation. Sequencing libraries were generated using an Agilent SureSelect Human All Exon Kit V5 (Agilent Technologies). Finally, 150-bp paired-end sequencing was performed on the Illumina HiSeq platform using the constructed DNA libraries for WES. All sequencing data were deposited in the Sequence Read Archive (SRA) under accession code SRP110150.

Bioinformatics data analysis

First, data quality control was performed, and all downstream bioinformatics analyses were based on high-quality clean data, in which reads containing an adapter, reads containing poly-N, and low-quality reads were removed.

Valid sequencing data were mapped to the reference human genome (UCSC hg19) using the Burrows-Wheeler Aligner (BWA; ref. 18). Somatic SNVs were detected with muTect (19). For the 9 pairs of cases with frozen tissues, a minimum of 14 reads covering a site in the tumor and 8 reads in the normal sample were required for mutation calling. Somatic Indels were detected using Strelka (20). High false-positive calls are usually present in FFPE specimens (21); therefore, only mutation sites with a minimum of 3 variant reads and a variant allele frequency (VAF) >0.08 were used for further analysis. To validate the variants, we performed ultra-deep WES (averaging 744 ×) on all samples from three patients (OS06, OS09, and OS10), and we viewed each variant using the integrative genomics viewer and filtered out the low credibility variants for other samples (22). To precisely determine the presence or absence of mutations in distinct samples from one patient, we used the second somatic mutation detection method “force calling” to detect SNVs in each sample based on the aggregate set of somatic events in one patient (22). ANNOVAR (23) was performed to annotate the variant call format (VCF). Polymorphisms, somatic SNVs and indels, referenced in the 1000 Genomes Project or Exome Aggregation Consortium (ExAC) with a minor allele frequency over 1% were removed. All of the detail-filtered SNVs and Indels are presented in Supplementary Fig. S1. FACETS (24) software was used to detect somatic CNVs. B-allele frequencies were used to infer loss-of-heterozygosity events. The GISTIC (25) algorithm was used to infer recurrently amplified or deleted genomic regions. To accurately detect the somatic breakpoint, we utilized CREST (26) software, which is based on soft-clipped approaches rather than paired-ends. Furthermore, all gene enrichment analyses were performed by clusterProfiler (27) using default parameters.

Spectrum and signatures of somatic mutations

Somatic SNVs were divided into 96 trinucleotides (mutated base plus its sequence context; refs. 28, 29) by 16 possible flanking nucleotide contexts. For each patient, we classified each mutation as a “primary tumor–specific mutation,” an “metastatic tumor–specific mutation” or a “shared mutation” based on whether it was located on the trunk or branch of the phylogenetic tree. The variants were also split according to their occurrence in primary tumors and metastatic tumors. Student t test was used to compare the mutation spectra of the six mutation types (C>A, C>G, C>T, T>A, T>C, and T>G) and 96 trinucleotides between these variants. The 30 known signatures were analyzed with deconstructSigs (30).

Phylogenetic tree construction

All “force calling” somatic events (SNVs/Indels) were used to build the osteosarcoma patient phylogenetic trees. Binary presence/absence matrices were built using the R Bioconductor package phangorn (31) from the regional distribution of variants within the patient. To find potential cancer driver genes, we first used MuSiC and MutSigCV to identify significantly and highly mutated genes in osteosarcoma genomes [false discovery rate (FDR) < 0.1]. Second, we constructed a list of potential driver genes (n = 833), which was comprised of all genes identified in the COSMIC cancer gene census (June 2016; ref. 32), DDR genes, plus those identified in previous osteosarcoma studies (5–7). Using the unsupervised hierarchical clustering method, we built the phylogenetic tree based on the genome-wide CNVs.

Cancer cell fraction estimation and cluster analysis

For each somatic mutation, the VAF was calculated using the number of reads supporting the variant allele (Rmut) and the number of reads supporting the reference allele (Rnorm; namely, VAF = Rmut/(Rmut + Rnorm)). Then, we calculated the value of cancer cell fractions (CCF) as follows: VAF = p*CCF/(CPNnorm (1-p) + p*CPNmut), where CPNmut indicates the local copy number in the tumor, CPNnorm indicates the local copy number in the normal controls (usually assumed to be 2 except for sex chromosomes), and p indicates the tumor purity in each sequenced sample. We applied ABSOLUTE to infer tumor purity. The VAF was defined as the VAF of each somatic mutation. All of the detailed CCF-filtered mutations in each patient are presented in Supplementary Fig. S2. This fraction is represented as a distribution between 0 and 1; a CCF value of 1 corresponds to a mutation present in 100% of the cancer cells in a sample, while a CCF value < 1 indicates that the mutation is present in a subset of the cancer cells in a sample and thus is subclonal. To portray the subclonal architecture, validated SNVs in 3 patients (OS06, OS09X, and OS10) were applied to infer the clone cluster according to PyClone (a Bayesian clustering method; ref. 33).

Identification of putative neoantigens

For each patient, we first applied the Polysolver algorithm to infer the individual HLA type, and then we used the identified nonsilent mutations to generate a list of mutant peptides of approximately 9–11 amino acids in length (34–35). Putative neoantigens were identified when the binding affinity of the mutant peptides was higher than its corresponding wild-type peptide and with a predicted binding strength of < 500 nmol/L (34, 35).

Statistical analysis

All statistical analyses were implemented in the R environment, and Fisher exact tests were two-tailed.

IHC

FFPE tissue slices were baked in a dryer at 60°C for 2 hours. Sections were then deparaffinized three times with xylene for 15 minutes, rehydrated with 100%, 95%, and 75% ethanol successively for 5 minutes, washed with PBS for 5 minutes, and then heated in citrate buffer at 95°C for 20 minutes and cooled to room temperature. Sections were then treated with 3% hydrogen peroxide solution for 10 minutes. Slices were incubated overnight at 4°C with antibodies specific to PD-L1 (Spring Bioscience SP142, diluted 1:30)/CD8 (Beijing Zhongshan Golden Bridge Biotechnology SP16, diluted 1:20) and visualized on a PV-9000 Enhanced Polymer System according to the manufacturer's instructions (Beijing Zhongshan Golden Bridge Biotechnology). Finally, slices were counterstained with hematoxylin.

Systemic treatments of patient OS06

Patient O206 received ifosfamide (3 g/m2, d1–4) plus doxorubicin (40 mg/m2, d1) for 12 cycles as first-line treatment, gemcitabine (1250 mg/m2, d1, 8) plus docetaxel (60 mg/m2, d1) for 4 cycles as second-line treatment, and pembrolizumab (240 mg, q3W × 6 cycles, self-purchased) as third-line treatment.

Genomic architecture of pulmonary metastatic tumors and matched primary tumors

We performed multiregional WES (median sequencing depth of 275 ×) on 86 tumor regions from resected primary tumors and metastatic tumors from 10 patients with osteosarcoma, with 9 patients providing matched freshly frozen samples and one providing paired FFPE specimens (Supplementary Fig. S3; Supplementary Table S2). The median number of fractionated regions of each tumor was 4 (range from 1 to 7). The WES results were further validated using ultra-deep WES, with a mean depth of 744 × in all tumor regions from three representative patients and a validation rate of 96.8% (patients OS06, OS09, and OS10). The nonsynonymous mutations in each sample were identified, with a median number of 28 (ranging from 3 to 67) in primary tumor regions and 99 (ranging from 42 to 177) in metastatic tumor regions (Supplementary Fig. S1).

We sought to illustrate the distributions of potential driver nonsynonymous SNVs/Indels in multiple regional samples of each patient. All patients were identified as harboring potential driver oncogenic or suppressor gene mutations (2–7 mutations, Supplementary Fig. S1). No overwhelming recurrent gene mutations were observed across different patients with osteosarcoma, except for TP53 mutations in three patients (OS04, OS05, and OS10), and the mutations showed significant distribution diversity among the primary tumor and metastatic tumor samples (Supplementary Fig. S1). These results suggest it is impossible to identify common targets and develop specific targeted drugs for osteosarcoma. To further explore osteosarcoma-related molecular events, genomic CNVs were analyzed: 6p12.3, 8q24.3, 12q13.11, 17p11.2, and 19q12 were the most common chromosome regions with gain copy numbers in both primary tumors and metastatic tumors, and 2q21.2 and 6p22.31 were the most commonly deleted regions in both primary tumors and metastatic tumors (Fig. 1A and B). The most recurrent genes with previously reported CNVs were also detected in our research, including loss of TP53 (3/10), RB1 (6/10), and CDKN2A (3/10) and gain of MYC (6/10) and MDM2 (6/10). Most of the CNVs were detected in all regions of one tumor as the early events, such as TP53 loss in OS09_MT, RB1 loss in OS06_MT3, MDM2 gain in OS01_PT, OS04_PT and OS03_MT. Some of the CNV events were detected in some but not all regions of one tumor (Supplementary Fig. S4). Generally, most recurrent CNV regions were shared in all primary tumor and/or metastatic tumor samples, indicating that CNVs are early molecular events in tumorigenesis and disease progression.

Figure 1.

The genomic disparity between primary tumor (PT) and metastatic tumor (MT). A, The overall picture of the CNVs in each sample in 9 patients based on the WES data; blue, the deletion copy number; red, the amplification copy number. The bar on the left indicates the primary tumor samples (blue) and metastatic tumor samples (yellow). The color bar on the right represents the different patients. B, The high-frequency analysis of CNVs from WES data and regions with q values < 0.05 are selectively labeled. Top, the amplification of primary and metastatic tumors; bottom, deletions. C, The tumor nonsynonymous mutation burden, including SNVs and small insertions and deletions (Indels) in each tumor of 10 patients. D, The average tumor nonsynonymous mutation burden of primary tumors and metastatic tumors. ***, significant difference (P < 0.001). E and F, The copy number of amplifications or deletions in each sample. G, Top, the tumor nonsynonymous mutation burden in each sample. The central heatmap shows alterations of DDR pathway–related genes in each sample, including nonsilent SNVs/Indels (blue) and frameshifts (yellow). The histogram on the right presents the number of alterations on every gene. The bar on the bottom presents the primary tumor (blue) samples and metastatic tumor (yellow) samples.

Figure 1.

The genomic disparity between primary tumor (PT) and metastatic tumor (MT). A, The overall picture of the CNVs in each sample in 9 patients based on the WES data; blue, the deletion copy number; red, the amplification copy number. The bar on the left indicates the primary tumor samples (blue) and metastatic tumor samples (yellow). The color bar on the right represents the different patients. B, The high-frequency analysis of CNVs from WES data and regions with q values < 0.05 are selectively labeled. Top, the amplification of primary and metastatic tumors; bottom, deletions. C, The tumor nonsynonymous mutation burden, including SNVs and small insertions and deletions (Indels) in each tumor of 10 patients. D, The average tumor nonsynonymous mutation burden of primary tumors and metastatic tumors. ***, significant difference (P < 0.001). E and F, The copy number of amplifications or deletions in each sample. G, Top, the tumor nonsynonymous mutation burden in each sample. The central heatmap shows alterations of DDR pathway–related genes in each sample, including nonsilent SNVs/Indels (blue) and frameshifts (yellow). The histogram on the right presents the number of alterations on every gene. The bar on the bottom presents the primary tumor (blue) samples and metastatic tumor (yellow) samples.

Close modal

Genomic disparities between primary osteosarcoma and the corresponding metastatic tumors

To examine the genetic variation after pulmonary metastasis, we systemically analyzed the disparity of mutations and genomic instability between primary tumors and metastatic tumors. The nonsynonymous mutational burden showed consistent elevation in metastatic tumors compared with that in primary tumors, except patient OS07, and the mean mutational burden in the metastatic tumor cohort was significantly higher than that in the primary tumor cohort (fold change = 3.0, P = 0.0003; Fig. 1C and D), indicating mutation accumulation during pulmonary metastasis of osteosarcoma. Furthermore, the variation of the total mutational burden in the exome regions between primary tumors and metastatic tumors was smaller than that of the nonsynonymous mutations (Supplementary Fig. S5), implying a more important role for nonsynonymous mutations during metastasis. We further analyzed the CNVs of the WES data and found steadily increased CNV frequencies in the metastatic tumors, especially deletion events (Fig. 1A, B, E, and F). In addition, a set of high-confidence SVs in the exon regions was found, including TP53 fusion as a driver event in patient OS08 and MLH3 fusion in patients OS02_PT_2, OS03_PT_3, and OS08_PT_3. MLH3 translocation in two samples (OS02_PT_2 and OS08_PT_3) had a common translocation partner of RABGAP1L at the same site. SV numbers were higher in metastatic tumors than in primary tumors in several patients (OS04, OS09, OS12, and OS02; Supplementary Fig. S1), although there was no significant difference in the mean SV number between primary tumors and metastatic tumors (P = 0.23; Supplementary Fig. S6). These results demonstrate that metastatic tumors have significantly higher genomic instability than primary tumors.

A deficiency of DNA damage response (DDR) genes is believed to be associated with a high mutational burden and genomic instability (36). On the basis of the results of our above analyses, we speculated that an association exists between DDR gene alterations and tumorigenesis and the development of osteosarcoma. Here, we identified more DDR gene deficiencies in metastatic tumors than in primary tumors (Fig. 1G). These results are in accordance with the trends of mutational burden, suggesting cascade reactions of the mutational burden and DDR gene inactivation during the metastasis process of osteosarcoma. In addition, all patients were found to harbor germline mutations in DDR genes (Supplementary Tables S3), including mismatch repair (MMR) gene mutations in three patients: MSH2 (OS04 and OS07), MSH6 (OS08), MLH1 (OS08), and PMS2 (OS04 and OS07). These results suggest that mutation accumulation due to DDR gene deficiency may contribute to the process of pulmonary metastasis and osteosarcoma tumorigenesis.

Evolutionary history of pulmonary metastatic osteosarcoma

To identify the key molecular events and their evolutionary process during tumorigenesis and metastasis, phylogenetic relation trees were constructed on the basis of the filtered nonsynonymous SNVs/Indels (Fig. 2A; Supplementary Fig. S1; refs. 31, 37). Potential driver events were labeled on the phylogenetic trees, displaying broad spatial and temporal tumor heterogeneity, which were identified according to the validated genetics (5) and high-throughput sequencing (6–9) studies of osteosarcoma and the Catalog of Somatic Mutations in Cancer (COSMIC) cancer gene census (38). To clearly illustrate the dynamic genetic evolution during tumor progression and pulmonary metastasis, five mutation categories were defined according to the distributions of gene alterations in different regions of phylogenetic relation trees, including “Trunk”, “PTun trunk”, “MTun trunk”, “primary tumor branch” and “metastatic tumor branch” (Fig. 2A). “PTun Trunk” and “PT branch” were classified into “PTun”, and “MTun Trunk” and “metastatic tumor branch” were classified into “MTun” (Fig. 2A). By multiregional deep sequencing in fractional tumor samples, early and late SNV/Indel events contributing to tumorigenesis (Trunk), disease progression (PTun trunk and MTun trunk), and pulmonary metastasis (MTun trunk and metastatic tumor branch) during osteosarcoma evolution were clearly displayed. Linear evolution (6/10 patients, OS01, OS04, OS05, OS06, OS08, and OS09) and parallel evolution (4/10 patients, OS02, OS03, OS07, and OS10) were the two major evolutionary patterns during pulmonary metastasis of osteosarcoma. Patients with the linear evolutionary pattern had more shared nonsynonymous mutations between primary tumors and metastatic tumors, in which the ancestral cells of the metastatic tumor might deviate at the later stage during the development of primary tumor. For patients with the parallel evolutionary pattern, a metastatic tumor may form at a very early stage during formation of the primary tumor, with extremely limited common mutations and no shared clonal mutations between the primary tumors and metastatic tumors (Supplementary Fig. S2). The survival interval from primary tumors resection to initial pulmonary metastasis was longer in patients with the linear evolution pattern than that in patients with the parallel evolution pattern (median: 26 months vs. 8.5 months, P = 0.012; Fig. 2A).

Figure 2.

The phylogenetic tree and intratumor heterogeneity in 10 pulmonary metastatic osteosarcomas. A, The integral phylogenetic trees of primary tumors (PT) and metastatic tumors (MT), the branch length of which is proportional to the number of nonsynonymous mutations (the detailed data can be found in Supplementary Fig. S1). The integral phylogenetic tree depicts the genetic flow and full lineage relationship between different tumor regions in one patient and the intratumor heterogeneity in primary tumor and metastatic tumor. All potential driver gene mutations are labeled on the phylogenetic tree, and red font indicates fusion genes. The interval on the bottom indicates the time interval from primary tumor resection to initial detection of pulmonary metastasis tumors. B, Gene set enrichment in the categories “Trunk,” “PTun,” and “MTun.” The size of the circle indicates the gene numbers mapped to the pathway, and the different colors indicate the q value. C, The percentage of trunk mutations and heterogeneity mutations in primary tumor and metastatic tumor. A higher percentage of heterogeneous mutations indicates greater tumor heterogeneity. The tumors of patients OS07_PT, OS07_MT, OS01_MT, and OS05_MT were sampled in only one region, and their data are not available.

Figure 2.

The phylogenetic tree and intratumor heterogeneity in 10 pulmonary metastatic osteosarcomas. A, The integral phylogenetic trees of primary tumors (PT) and metastatic tumors (MT), the branch length of which is proportional to the number of nonsynonymous mutations (the detailed data can be found in Supplementary Fig. S1). The integral phylogenetic tree depicts the genetic flow and full lineage relationship between different tumor regions in one patient and the intratumor heterogeneity in primary tumor and metastatic tumor. All potential driver gene mutations are labeled on the phylogenetic tree, and red font indicates fusion genes. The interval on the bottom indicates the time interval from primary tumor resection to initial detection of pulmonary metastasis tumors. B, Gene set enrichment in the categories “Trunk,” “PTun,” and “MTun.” The size of the circle indicates the gene numbers mapped to the pathway, and the different colors indicate the q value. C, The percentage of trunk mutations and heterogeneity mutations in primary tumor and metastatic tumor. A higher percentage of heterogeneous mutations indicates greater tumor heterogeneity. The tumors of patients OS07_PT, OS07_MT, OS01_MT, and OS05_MT were sampled in only one region, and their data are not available.

Close modal

On the basis of the mutation categories, we performed gene function enrichment analyses to better understand the potential impact of the signaling pathways on tumorigenesis and pulmonary metastasis (Fig. 2B; Supplementary Fig. S7). The PI3K–Akt signaling pathway was significantly enriched in all mutation categories except for the “metastatic tumors branch,” indicating that activation of the PI3K–Akt pathway may contribute to initial osteosarcoma tumorigenesis and early pulmonary metastasis. In contrast, the MAPK signaling pathway was enriched in the mutation categories “Trunk”, “MTun”, and “MT branch” but not “PTun” or “primary tumor branch”, suggesting that this pathway may be critical in later stages of pulmonary metastasis. In addition, distinct from the “PTun” category, the “MTun” category demonstrated more significant enrichment in chemo-drug resistance-related pathways, such as “platinum drug resistance” and “chemical carcinogenesis,” as well as tumor development pathways, such as the “Hippo signaling pathway.” Most of the mutations enriched in chemo-resistance were in the “MTun trunk” category, including PIK3R5 (OS06_MT3), TP53 (OS10_MT), TOP2B (OS06_MT2, OS06_MT3), and MGST1 (OS10_MT), and others were in the “MTun branch” category, including BBC3 (OS04_MT_1), TOP2B (OS04_MT_2), and ABCC2 (OS10_MT_3).

Tumor heterogeneity and clonal evolution of pulmonary metastatic osteosarcoma

To evaluate the tumor heterogeneity of pulmonary metastatic osteosarcoma, we analyzed the intertumoral heterogeneity (IrTH) between primary tumors and metastatic tumors and the IaTH of metastatic tumors or primary tumors. The percentage of trunk mutations showed an increasing tendency from primary tumors (mean, 53.82%) to metastatic tumors (mean, 74.14%), indicating that the proportion of heterogeneous alterations decreases (P = 0.05) after metastasis (Fig. 2C). To further evaluate the spatial and temporal heterogeneity, Jaccard similarity coefficients (39) were estimated among all tumor regions in each patient, revealing more significant IrTH rather than IaTH (Supplementary Fig. S8). The nonsynonymous mutations maintained a higher degree of similarity among different fractioned regions in each metastatic tumor lesion than those in matched primary tumor lesion.

To better understand the IaTH and genetic variation dynamics during pulmonary metastasis, we applied an n-dimensional Bayesian Dirichlet process (33) to portray the clonal architecture and to analyze the genetic evolutionary trajectory of subclones in each patient. The mutations in the trunk, PTun trunk, and MTun trunk categories tended to be clones and those in the primary tumor–private and metastatic tumor–private categories were subclones. Next, we analyzed the clone status of potential driver mutations (Supplementary Fig. S9), more than half of which exhibited subclones. Potential driver genes on the trunk, PTun trunk, and MTun trunk, such as FNACC, FANCI, and TP53, exhibited clones in all regions of the tumor. The results of clonality analysis further confirmed the above investigation of the phylogenetic tree, which showed that mutations on the trunk were probably earlier genetic alteration events, followed by those on the PTun trunk and MTun trunk, and the primary tumor–private and metastatic tumor–private mutations occurred later during tumorigenesis and osteosarcoma development. In patients with the linear evolution pattern, all subclones in the primary tumors arose from a common ancestor clone, which achieved metastatic potential disseminating from the primary tumor to the metastatic tumor (such as in patients OS06 and OS09) and subsequently from the metastatic tumor to another metastatic tumor (such as in patient OS06; Fig. 3A–D). The aberrant driver gene-defining subclones, such as NRF1, MDM2, PRAMEF4, HLF, and FNDC8 in patient OS06 and ERC1, CACNG8, and IKZF1 in patient OS09, may contribute to metastatic progression. For patients with the parallel evolutionary pattern, parallel polyclone origination was the main formation pattern of metastatic tumors, and the major subclone of metastatic tumors was already present in the primary tumor independent from the major clone, even though it was a minor clone (Fig. 3E and F). Interestingly, the parallel evolutionary pattern also existed in the late stage of pulmonary metastasis from metastatic tumors to metastatic tumors in linear evolutionary patients, such as subclone HLA-G and PNPO in patient OS06, suggesting that linear and parallel evolutionary patterns may coexist in one patient.

Figure 3.

Subclonal evolution model during pulmonary metastasis of patient OS06, OS09, OS10. A, C, and E, The CCF containing one mutation in the different regions. One representative region of each tumor is shown. The different colors indicate different clusters, and the potential driver genes are labeled. B, D, and F, The models of subclone evolution during pulmonary metastasis; different colors represent different subclones.

Figure 3.

Subclonal evolution model during pulmonary metastasis of patient OS06, OS09, OS10. A, C, and E, The CCF containing one mutation in the different regions. One representative region of each tumor is shown. The different colors indicate different clusters, and the potential driver genes are labeled. B, D, and F, The models of subclone evolution during pulmonary metastasis; different colors represent different subclones.

Close modal

Dynamic mutational spectra and signatures during pulmonary metastasis

To determine the mutational spectra and signatures associated with pulmonary metastasis, we evaluated the base substitution spectrum based on the “Trunk”, “PTun” and “MTun” mutation categories. “PTun” displayed a similar mutational spectrum with “MTun,” both of which showed a significantly higher percentage of C>A (P = 0.0003) but fewer percentages of A>G (P = 0.0012) and C>T (P = 0.0005). These results suggest that C>T and A>G transitions may be important events associated with tumorigenesis and that C>A transversions may contribute to tumor progression (Fig. 4A).

Figure 4.

The temporal dynamic somatic mutation spectrum and signature. A, The six categories of possible base substitution distributions of unique primary tumor (PT) mutations, primary tumor and metastatic tumor (MT) shared mutations and unique metastatic tumor mutations in the 9 patients; significant differences according to P values are labeled. B, The known mutation signature enrichment in the “Trunk,” “PTun,” and “MTun” categories; the medium bar indicates the significant enrichment signature with a medium value. C, Comparison of the six categories of possible base substitution distributions of primary tumor mutations and metastatic tumor mutations in the 9 patients, and the distributions of lung adenocarcinoma (LUAD) and lung squamous carcinoma (LUSC) from samples in the The Cancer Genome Atlas database.

Figure 4.

The temporal dynamic somatic mutation spectrum and signature. A, The six categories of possible base substitution distributions of unique primary tumor (PT) mutations, primary tumor and metastatic tumor (MT) shared mutations and unique metastatic tumor mutations in the 9 patients; significant differences according to P values are labeled. B, The known mutation signature enrichment in the “Trunk,” “PTun,” and “MTun” categories; the medium bar indicates the significant enrichment signature with a medium value. C, Comparison of the six categories of possible base substitution distributions of primary tumor mutations and metastatic tumor mutations in the 9 patients, and the distributions of lung adenocarcinoma (LUAD) and lung squamous carcinoma (LUSC) from samples in the The Cancer Genome Atlas database.

Close modal

We further performed deconstructSigs (30) analyses to sort the base substitution spectra in the “Trunk”, “PTun,” and “MTun” mutation categories based on known mutational signatures (Fig. 4B; Supplementary Fig. S10; refs. 28, 29). Signature 3 (associated with failure of DNA double-strand break repair by homologous recombination) was significantly enriched in all stages during tumor evolution, suggesting that DNA double-strand break repair by homologous recombination may be related to osteosarcoma tumorigenesis. Increased enrichment of signature 6 (defective DNA mismatch repair) was present in the early stage (Trunk), whereas signatures 4 (smoking), 15 (defective DNA mismatch repair), 24 (exposure to aflatoxin), and 29 (tobacco chewing) were present in the late stage (PTun and MTun). These results suggest that failure of DNA damage repair might be crucial for mutation accumulation during tumorigenesis and osteosarcoma development.

The spectra and signatures of all mutations in the primary tumors and metastatic tumors were also analyzed (Fig. 4C; Supplementary Fig. S11). Signatures 1, 3, 4, 6, 12, 24, and 29 were enriched in both the primary tumors and metastatic tumors, and no significant mutation signature was observed in the metastatic tumors. To determine the relationships between the mutational spectra and tumor organ sites, we analyzed the spectra of primary tumors, metastatic tumors and primary lung cancer. C>T and C>A were the most common base substitutions, which were similar between primary tumors and metastatic tumors, but significantly different from primary lung cancer (from The Cancer Genome Atlas database; Fig. 4C; Supplementary Fig. S12). These results suggest that C>T and C>A base substitutions may be related to osteosarcoma tumorigenesis and that osteosarcoma cells maintain their primary mutational traits after pulmonary metastasis.

Improved immunogenicity of metastatic tumors could raise the effectiveness of the immune checkpoint blockade

Currently, pulmonary metastatic osteosarcoma is refractory to chemotherapy and lacks an efficient therapeutic strategy (1, 40). The elevated mutational load and genomic instability observed in pulmonary metastatic tumors compared with primary tumors indicated improved immunogenicity (41–42), indicating the potential application value of immune checkpoint blockades (ICBs) such as PD-1 antibodies, etc. We further analyzed several indices previously reported to be associated with the efficacy of PD-1/PD-ligand-1 (PD-L1) antibodies (43–44). First, we estimated the neoantigen burden in each tumor (35). The metastatic tumor cohort (median = 147) had a significantly higher neoantigen number than the primary tumor cohort (median = 45, fold change = 3.3, P < 0.001, Fig. 5A and B). Second, metastatic tumors had higher positive rates of PD-L1 expression in tumor cells and tumor-infiltrating lymphocytes (TIL) than primary tumors (7 vs. 3 patients for PD-L1–positive, and 10 vs. 6 patients for TIL-positive, Fig. 5C; Supplementary Fig. S13). Third, the immune-regulatory genes covered by the gain regions based on the genomic CNV data (45) were more frequently found in the metastatic tumors cohort than the primary tumor cohort (Supplementary Fig. S14). These results provide evidence that the use of immune checkpoint inhibitors, such as PD-1 antibodies, could be a promising strategy for pulmonary metastatic osteosarcoma. Correspondingly, one patient (OS06) who relapsed after initial primary tumor surgery and subsequent resection of pulmonary metastasis with positive PD-L1 in metastatic tumors, but not in primary tumors received pembrolizumab treatment (240 mg/3 weeks) for 6 cycles after failure to respond to multiple lines of chemotherapy. The pulmonary metastasis showed a partial response, with some pulmonary legions fading away (Supplementary Fig. S15).

Figure 5.

The immunogenicity across primary tumors (PT) and metastatic tumors (MT). A, The number of all predicted neoantigens in each tumor of 10 patients. B, The difference in the mean number of potential neoantigens in primary tumors and metastatic tumors. ***, significant difference (P < 0.001). C, IHC results represent the expression of PD-L1 in tumor tissue and CD8+ tumor-infiltrating lymphocytes in the primary tumor and metastatic tumor in patient OS06.

Figure 5.

The immunogenicity across primary tumors (PT) and metastatic tumors (MT). A, The number of all predicted neoantigens in each tumor of 10 patients. B, The difference in the mean number of potential neoantigens in primary tumors and metastatic tumors. ***, significant difference (P < 0.001). C, IHC results represent the expression of PD-L1 in tumor tissue and CD8+ tumor-infiltrating lymphocytes in the primary tumor and metastatic tumor in patient OS06.

Close modal

Our comprehensive investigation of the mutational landscapes of 86 tumor regions yielded a detailed genetic evolutionary history of the tumorigenesis, progression, and pulmonary metastasis of osteosarcoma.

In this study, which was limited by small patient numbers, both linear and parallel evolutionary patterns were observed. The linear evolutionary pattern was verified in 6 of 10 patients with osteosarcoma during pulmonary metastasis in our study, in which the primary tumors and metastatic tumors shared substantial genetic alterations and had a common ancestor. A critical assumption of the linear progression model is that only cancer cells with a genetic evolutionary advantage can effectively disseminate and colonize, which underlies the importance of targeting trunk events in this subset of patients (13). In contrast to the linear evolution pattern, in 4 of 10 patients, we did not identify any shared clone mutations between the primary tumor and metastatic tumor. While considering the results of mutational signatures and clinical pathology, which revealed that all four metastatic tumors had osteosarcoma pulmonary metastasis, we surmised that the primary tumors and metastatic tumors in these four cases may have followed a parallel progression model, originating from independently synchronous tumor cells at a much earlier time during tumorigenesis, which has rarely been reported in other cancers (14–17, 46), indicating a distinct mechanism in osteosarcoma pulmonary metastasis. Theoretically, patients following the linear evolutionary pattern should exhibit longer disease-free survival (the time from surgery to pulmonary metastasis) than those with the parallel evolutionary pattern, which was confirmed in our patient cohort, further proving that parallel evolution may be associated with earlier pulmonary metastasis than linear evolution. Moreover, considering the early onset of pulmonary metastasis in a substantial proportion of patients with osteosarcoma, screening for such metastases should be performed at initial diagnosis, highlighting the potential clinical applicability of differentiating tumor evolution patterns in patients with osteosarcoma.

Tumor heterogeneity is one common factor that may influence the clinical value of small biopsy specimens. Zhang and colleagues reported that the majority of driver mutations could be detected in biopsy specimens of lung cancer (47), suggesting that genotyping based on biopsy specimens was generally sufficient for guiding lung cancer treatment. Moreover, our microdissection research provided a more integrated mutational profile than single-sample analysis and evaluated tumor heterogeneity in pulmonary metastatic osteosarcoma. The results suggested that the molecular characteristics were similar among different intrametastatic tumor spatial regions. We also discovered that genetic aberrations in the trunk regions accounted for a higher percentage of total gene mutations in metastatic tumors than primary tumors. These results, combined with our verified significant disparity between primary tumors and metastatic tumors, indicate that biopsy samples for genotyping are sufficient for reflecting the whole-body genetic characteristics of pulmonary metastatic tumors, which are necessary for differentiation from the primary tumor. However, the significant genomic disparity emphasizes the necessity of real-time genetic analysis after pulmonary metastasis. Furthermore, in our study, most patients received neoadjuvant chemotherapy before primary lesion resection as part of routine clinical practice, which may affect the genomic profile and lead to an elevated tumor nonsynonymous mutation burden (TMB). However, the median period of four cycles (approximately 2–3 months) of neoadjuvant chemotherapy did not seem to significantly change the genome profiles of primary tumors. Importantly, only 3 of 10 patients underwent chemotherapy before surgery for their metastatic tumors, and all of these patients underwent only 2 cycles of chemotherapy, indicating that chemotherapy may have a more significant effect on primary tumors versus metastatic tumors in terms of genomic instability. However, metastatic tumors were found to have a significantly greater TMB and neoantigen load than primary tumors. Thus, we speculated that the genomic disparities might be more significant if no neoadjuvant chemotherapy was delivered. We think the genomic sequencing results of the primary osteosarcoma tumors and matched metastatic tumors are representative and reflect the real status in current clinical practice. On the other hand, the results of the mutational spectrum and the observation of signatures mostly enriched in mutations related to the DDR rather than chemotherapy also implied the slight effect of chemotherapy compared with the DNA repair system on genomic instability.

In addition, no predominantly driver oncogenic mutations were found except for TP53, consistent with previous reports (5–9), suggesting that current target therapies may not be available in most metastatic osteosarcomas. In our study, the variants of TP53-related genes were found in 7 of 10 patients (OS02, OS04, OS05, OS07, OS08, OS09, and OS10), which is accordance with some prior reports (8, 9) but lower than those in other studies (6–7). This finding may be due to the limited sample size. However, critical signaling pathways were identified as enriched in different evolutionary stages. On the basis of our findings, the PI3K–Akt signaling pathway may play an important role at both the early and late stages during osteosarcoma evolution, and the MAPK signaling pathway may be important at the metastatic stage, which implies that blocking the PI3K–Akt pathway may be one therapeutic strategy for osteosarcoma, in accordance with a previous report (6), and MAPK inhibitors may be important for pulmonary metastatic osteosarcoma. Moreover, germline mutations of DDR genes were found in all patients in this cohort, such as DNA damage repair genes, ATM and BRCA2, and MMR genes, MSH2, MSH6, MLH1, and PMS2. The incidence of germline DNA damage repair gene mutations had been reported significantly higher in patients with metastatic prostate than those with localized prostate cancer (48). However, the association of germline DDR gene mutation with osteosarcoma metastasis has not been well established. To the best our knowledge, our study represents the first report of the potential association in osteosarcoma. We did not detect a significantly high Indel mutational burden and microsatellite instability in these three cases with germline MMR mutations (Supplementary Fig. S16), although the deleterious score was very high, with an SIFT score of 0 and a polyphen2 score of 1 for these mutations (Supplementary Table S3). MMR germline mutations may not lead to microsatellite instability in osteosarcoma, and these mutations were correlated with tumor metastasis as reported in prostate cancer (48). Given the limited sample size of our study, the results need to be confirmed in a larger cohort. However, the consistent observations in our study suggest that mutation accumulation due to DDR deficiency might contribute to the process of osteosarcoma pulmonary metastasis.

Previous studies involving genomic base substitutions with a limited number of osteosarcoma cases have demonstrated that C>T is the most common mutation (28). Here, the mutational spectrum analyses illustrated that C>T and C>A were the most common base substitutions of osteosarcoma, which were similar between primary tumors and metastatic tumors. These data suggest that C>T may be the pivotal point mutation contributing to osteosarcoma tumorigenesis. However, further studies are needed to verify these findings. The mutational spectrum of pulmonary metastasis was different from lung cancers, suggesting that osteosarcoma cells maintain their primary mutational traits after pulmonary metastasis with limited influence from the lung microenvironment. On the other hand, the DNA damage response-related signatures, such as signatures 3, 6, and 15, were significantly enriched in both the early and late stages of osteosarcoma development. In particular, signature 3 was associated with failure of DNA double-strand break repair by homologous recombination, which may explain the recurrence of chromosomal instability in osteosarcoma (7). These results indicate that DDR dysfunction may play an important role during osteosarcoma tumor evolution and imply new therapeutic strategies for pulmonary metastatic osteosarcoma.

One pivotal finding of this study is the significant biological disparity identified between pulmonary metastasis and primary osteosarcoma, which mainly includes tumor immunogenicity and genomic instability. These results strongly indicate that the disparity should be considered sufficiently before formulating treatment regimens for pulmonary metastasis. Sarcomas have generally been identified to harbor a relatively low mutational load on the basis of previous reports (28) and are not considered candidates for immunotherapy. Although the mutation accumulation effect is not rare in metastatic tumors of other cancer types (14–17, 46), such prevalent and high augmentation of mutations as observed in pulmonary metastatic osteosarcoma have not been reported. At present, sarcomas have generally been identified to harbor a relatively low mutational load on the basis of previous reports (28) and are not considered candidates for immunotherapy. However, we provide solid evidence demonstrating dramatically improved immunogenicity after pulmonary metastasis compared with primary osteosarcoma, as suggested by a series of potential immune checkpoint inhibitor–related biomarkers (36, 43, 49), including elevated levels of the mutational load, a neoantigen burden, PD-L1 expression and T-cell infiltration. The tumor microenvironment was determined by both the intrinsic traits of cancer and the surrounding environment. Non–small cell lung cancer treatment has been reported to significantly benefit from immune checkpoint inhibitors (49–50), which may be due to the lung environment to some degree. However, no reports are available regarding related research on pulmonary metastatic tumors. PD-L1 expression has been considered one of the most important biomarkers for checkpoint inhibitor therapy and has been utilized in many studies (50–51), even in the clinical trials for an anti-PD-1 inhibitor (50). Therefore, we evaluated the expression of PD-L1 and not PD-1 in this study. Although the common conclusion that pulmonary metastasis is associated with higher levels of PD-L1 expression and TILs than those in primary tumors could not be confirmed on the basis of our limited samples, we believe that our convincing observations are representative to some degree. Further studies with larger sample sizes are warranted and should be performed to confirm our results. On the other hand, improved immunogenicity was one result of genomic instability (e.g., frequent SNVs/indels and CNVs), which was identified due to DDR deficiency (52). In addition to the high frequency of somatic mutations in DDR-related genes, germline mutations in MMR genes were also observed. We therefore speculate that mutation accumulation due to germline MMR mutations is initiated ahead of tumor emergence, and subsequent DDR deficiency could cause the cancer cells to be more aggressive and contribute to the subsequent process of pulmonary metastasis (48). These results strongly indicate that the disparity should be considered sufficiently before formulating treatment regimens for pulmonary metastasis.

Considering the limitations of this study, which were its relatively small sample size, our study should be interpreted carefully and regarded as descriptive observations. Therefore, studies with much larger cohorts (ideally with comprehensive real-time analyses of potentially predictive biomarkers and repeat biopsies at pulmonary relapse) are needed to fully unveil the general correlation between primary and pulmonary metastatic osteosarcoma and the immunogenicity of osteosarcoma and to explore the clinical relevance of immunotherapy.

In conclusion, this study describes the detailed evolutionary genomic landscape of and disparity between osteosarcoma and the corresponding metastatic tumors. The evidence for genetic and immunogenic evolution during pulmonary metastasis identified in our study strongly indicates that metastatic tumors should be treated separately from their original tumors, contributing to personalized therapy for pulmonary metastatic osteosarcoma and emphasizing the necessity of real-time genetic analysis after pulmonary metastasis. Despite the limited sample size, this study provides evidence suggesting that immune checkpoint inhibitors, such as PD-1 antibodies, represent promising clinical strategies for the treatment of pulmonary metastatic osteosarcoma.

No potential conflicts of interest were disclosed.

Conception and design: J. Wang, D. Wang, X. Niu, Z. Wang, Z. Huang, X.S. Xie, T. An, S. Gao, L. Wang

Development of methodology: D. Wang, C.-L. Song, X.S. Xie, Y. Tian, L. Wang

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): D. Wang, X. Niu, Z. Wang, C.-L. Song, Z. Huang, K.-N. Chen, J. Duan, J. Zhao, M. Zhuo, X. Kang, Y. Tian, J.-F. Han, T. An, Y. Sun, J. Zhao

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): D. Wang, Z. Wang, C.-L. Song, J. Duan, Y. Wang, L. Cai, J.-F. Han

Writing, review, and/or revision of the manuscript: D. Wang, Z. Wang, C.-L. Song, J. Duan, J. Xu, J.-F. Han, J. Ying, L. Wang

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): D. Wang, X. Niu, Z. Wang, K.-N. Chen, H. Bai, M. Zhuo, L. Cai, S. Gao, J. Ying, L. Wang, J. He

Study supervision: J. Wang, K.-N. Chen, T. An, L. Wang

We thank the patients and their families who agreed to provide samples to support this study, and the physicians and hospital staff who collected these samples. J. Wang received the grant from the National Natural Sciences Foundation Key Program of China (81630071 and 81330062), the National Key Research and Development Project of Precision Medicine Special Research of China (2016YFC0902300), National High Technology Research and Development Program 863 (SS2015AA020403), the CAMS Innovation Fund for Medical Sciences (CIFMS 2016-I2M-3-008). Z. J. Wang received the grant from the National Natural Sciences Foundation, China (81472206/81871889), and Beijing Natural Science Foundation, China (7172045), Beijing Nova Program Grants (Z141107001814051, Z181100006218130). D. Wang, Y. H. Tian, and J. Wang were supported in part by the Peking-Tsinghua Center Life Sciences. K. Chen received the grant from Beijing Municipal Administration of Hospitals Incubating Program (PX2018044). X. Kang received the grant from the Excellent Talent Training Project of Department of Beijing Municipal Organization of Communist Party of China (2016000021469G185). X. Kang and K. Chen received the grant from National Key R&D Program of China (2017YFC0907500).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

1.
Marina
N
,
Gebhardt
M
,
Teot
L
,
Gorlick
R
. 
Biology and therapeutic advances for pediatric osteosarcoma
.
Oncologist
2004
;
9
:
422
41
.
2.
Kager
L
,
Zoubek
A
,
Potschger
U
,
Kastner
U
,
Flege
S
,
Kempf-Bielack
B
. 
Primary metastatic osteosarcoma: presentation and outcome of patients treated on neoadjuvant Cooperative Osteosarcoma Study Group protocols
.
J Clin Oncol
2003
;
21
:
2011
8
.
3.
Harting
MT
,
Blakely
ML
. 
Management of osteosarcoma pulmonary metastases
.
Semin Pediatr Surg
2006
;
15
:
25
9
.
4.
Bacci
G
,
Rocca
M
,
Salone
M
,
Balladelli
A
,
Ferrari
S
,
Palmerini
E
. 
High grade osteosarcoma of the extremities with lung metastases at presentation: treatment with neoadjuvant chemotherapy and simultaneous resection of primary and metastatic lesions
.
J Surg Oncol
2008
;
98
:
415
20
.
5.
Martin
JW
,
Squire
JA
,
Zielenska
M
. 
The genetics of osteosarcoma
.
Sarcoma
2012
;
2012
:
1
11
.
6.
Perry
JA
,
Kiezun
A
,
Tonzi
P
,
Van Allen
EM
,
Carter
SL
,
Baca
SC
. 
Complementary genomic approaches highlight the PI3K/mTOR pathway as a common vulnerability in osteosarcoma
.
Proc Natl Acad Sci U S A
2014
;
111
:
E5564
73
.
7.
Chen
X
,
Bahrami
A
,
Pappo
A
,
Easton
J
,
Dalton
J
,
Hedlund
E
. 
Recurrent somatic structural variations contribute to tumorigenesis in pediatric osteosarcoma
.
Cell Rep
2014
;
7
:
104
12
.
8.
Kovac
M
,
Blattmann
C
,
Ribi
S
,
Smida
J
,
Mueller
NS
,
Engert
F
. 
Exome sequencing of osteosarcoma reveals mutation signatures reminiscent of BRCA deficiency
.
Nat Commun
2015
;
6
:
8940
.
9.
Behjati
S
,
Tarpey
PS
,
Haase
K
,
Ye
H
,
Young
MD
,
Alexandrov
LB
. 
Recurrent mutation of IGF signalling genes and distinct patterns of genomic rearrangement in osteosarcoma
.
Nat Commun
2017
;
8
:
15936
.
10.
Oda
Y
,
Yamamoto
H
,
Kohashi
K
,
Yamada
Y
,
Iura
K
,
Ishii
T
. 
Soft tissue sarcomas: from a morphological to a molecular biological approach
.
Pathol Int
2017
;
67
:
435
46
.
11.
Gupta
GP
,
Massague
J
. 
Cancer metastasis: building a framework
.
Cell
2006
;
127
:
679
95
.
12.
Greaves
M
,
Maley
CC
. 
Clonal evolution in cancer
.
Nature
2012
;
481
:
306
13
.
13.
Naxerova
K
,
Jain
RK
. 
Using tumour phylogenetics to identify the roots of metastasis in humans
.
Nat Rev Clin Oncol
2015
;
12
:
258
72
.
14.
Kim
TM
,
Jung
SH
,
An
CH
,
Lee
SH
,
Baek
IP
,
Kim
MS
. 
Subclonal genomic architectures of primary and metastatic colorectal cancer based on intratumoral genetic heterogeneity
.
Clin Cancer Res
2015
;
21
:
4461
72
.
15.
Brastianos
PK
,
Carter
SL
,
Santagata
S
,
Cahill
DP
,
Taylor-Weiner
A
,
Jones
RT
. 
Genomic characterization of brain metastases reveals branched evolution and potential therapeutic targets
.
Cancer Discov
2015
;
5
:
1164
77
.
16.
Hong
MKH
,
Macintyre
G
,
Wedge
DC
,
Van Loo
P
,
Patel
K
,
Lunke
S
. 
Tracking the origins and drivers of subclonal metastatic expansion in prostate cancer
.
Nat Commun
2015
;
6
:
6605
.
17.
Vignot
S
,
Frampton
GM
,
Soria
JC
,
Yelensky
R
,
Commo
F
,
Brambilla
C
. 
Next-generation sequencing reveals high concordance of recurrent somatic alterations between primary tumor and metastases from patients with non-small-cell lung cancer
.
J Clin Oncol
2013
;
31
:
2167
72
.
18.
Li
H
,
Durbin
R
. 
Fast and accurate long-read alignment with Burrows-Wheeler transform
.
Bioinformatics
2010
;
26
:
589
95
.
19.
Cibulskis
K
,
Lawrence
MS
,
Carter
SL
,
Sivachenko
A
,
Jaffe
D
,
Sougnez
C
. 
Sensitive detection of somatic point mutations in impure and heterogeneous cancer samples
.
Nat Biotechnol
2013
;
31
:
213
9
.
20.
Saunders
CT
,
Wong
WS
,
Swamy
S
,
Becq
J
,
Murray
LJ
,
Cheetham
RK
. 
Strelka: accurate somatic small-variant calling from sequenced tumor-normal sample pairs
.
Bioinformatics
2012
;
28
:
1811
7
.
21.
Yost
SE
,
Smith
EN
,
Schwab
RB
,
Bao
L
,
Jung
H
,
Wang
X
. 
Identification of high-confidence somatic mutations in whole genome sequence of formalin-fixed breast cancer specimens
.
Nucleic Acids Res
2012
;
40
:
e107
.
22.
Stachler
MD
,
Taylor-Weiner
A
,
Peng
S
,
McKenna
A
,
Agoston
AT
,
Odze
RD
. 
Paired exome analysis of Barrett's esophagus and adenocarcinoma
.
Nat Genet
2015
;
47
:
1047
55
.
23.
Wang
K
,
Li
M
,
Hakonarson
H
. 
ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data
.
Nucleic Acids Res
2010
;
38
:
e164
.
24.
Shen
R
,
Seshan
VE
. 
FACETS: allele-specific copy number and clonal heterogeneity analysis tool for high-throughput DNA sequencing
.
Nucleic Acids Res
2016
;
44
:
e131
.
25.
Boeva
V
,
Popova
T
,
Bleakley
K
,
Chiche
P
,
Cappo
J
,
Schleiermacher
G
. 
Control-FREEC: a tool for assessing copy number and allelic content using next-generation sequencing data
.
Bioinformatics
2012
;
28
:
423
5
.
26.
Wang
J
,
Mullighan
CG
,
Easton
J
,
Roberts
S
,
Heatley
SL
,
Ma
J
. 
CREST maps somatic structural variation in cancer genomes with base-pair resolution
.
Nat Methods
2011
;
8
:
652
4
.
27.
Yu
G
,
Wang
LG
,
Han
Y
,
He
QY
. 
clusterProfiler: an R package for comparing biological themes among gene clusters
.
OMICS
2012
;
16
:
284
7
.
28.
Alexandrov
LB
,
Nik-Zainal
S
,
Wedge
DC
,
Aparicio
SAJR
,
Behjati
S
,
Biankin
AV
. 
Signatures of mutational processes in human cancer
.
Nature
2013
;
500
:
415
21
.
29.
Helleday
T
,
Eshtad
S
,
Nik-Zainal
S
. 
Mechanisms underlying mutational signatures in human cancers
.
Nat Rev Genet
2014
;
15
:
585
98
.
30.
Rosenthal
R
,
McGranahan
N
,
Herrero
J
,
Taylor
BS
,
Swanton
C
. 
deconstructSigs: delineating mutational processes in single tumors distinguishes DNA repair deficiencies and patterns of carcinoma evolution
.
Genome Biol
2016
;
17
:
31
.
31.
Schliep
KP
. 
phangorn: phylogenetic analysis in R
.
Bioinformatics
2011
;
27
:
592
3
.
32.
Futreal
PA
,
Coin
L
,
Marshall
M
,
Down
T
,
Hubbard
T
,
Wooster
R
. 
A census of human cancer genes
.
Nat Rev Cancer
2004
;
4
:
177
83
.
33.
Nik-Zainal
S
,
Van Loo
P
,
Wedge
DC
,
Alexandrov
LB
,
Greenman
CD
,
Lau
KW
. 
The life history of 21 breast cancers
.
Cell
2012
;
149
:
994
1007
.
34.
Shukla
SA
,
Rooney
MS
,
Rajasagi
M
,
Tiao
G
,
Dixon
PM
,
Lawrence
MS
. 
Comprehensive analysis of cancer-associated somatic mutations in class I HLA genes
.
Nat Biotechnol
2015
;
33
:
1152
8
.
35.
Lundegaard
C
,
Lamberth
K
,
Harndahl
M
,
Buus
S
,
Lund
O
,
Nielsen
M
. 
NetMHC-3.0: accurate web accessible predictions of human, mouse and monkey MHC class I affinities for peptides of length 8–11
.
Nucleic Acids Res
2008
;
36
:
W509
12
.
36.
Diaz
LJ
Le DT PD-1 blockade in tumors with mismatch-repair deficiency
.
N Engl J Med
2015
;
373
:
1979
.
37.
Nixon
KC
. 
The parsimony ratchet, a new method for rapid parsimony analysis
.
Cladistics
1999
;
15
:
407
14
.
38.
Forbes
SA
,
Beare
D
,
Gunasekaran
P
,
Leung
K
,
Bindal
N
,
Boutselakis
H
. 
COSMIC: exploring the world's knowledge of somatic mutations in human cancer
.
Nucleic Acids Res
2015
;
43
:
D805
11
.
39.
Makohon-Moore
AP
,
Zhang
M
,
Reiter
JG
,
Bozic
I
,
Allen
B
,
Kundu
D
. 
Limited heterogeneity of known driver gene mutations among the metastases of individual patients with pancreatic cancer
.
Nat Genet
2017
;
49
:
358
66
.
40.
Reed
DR
,
Hayashi
M
,
Wagner
L
,
Binitie
O
,
Steppan
DA
,
Brohl
AS
. 
Treatment pathway of bone sarcoma in children, adolescents, and young adults
.
Cancer
2017
;
123
:
2206
18
.
41.
Riaz
N
,
Morris
L
,
Havel
JJ
,
Makarov
V
,
Desrichard
A
,
Chan
TA
. 
The role of neoantigens in response to immune checkpoint blockade
.
Int Immunol
2016
;
28
:
411
9
.
42.
Pritchard
AL
,
Burel
JG
,
Neller
MA
,
Hayward
NK
,
Lopez
JA
,
Fatho
M
. 
Exome sequencing to predict neoantigens in melanoma
.
Cancer Immunol Res
2015
;
3
:
992
8
.
43.
Rizvi
NA
,
Hellmann
MD
,
Snyder
A
,
Kvistborg
P
,
Makarov
V
,
Havel
JJ
. 
Cancer immunology. Mutational landscape determines sensitivity to PD-1 blockade in non-small cell lung cancer
.
Science
2015
;
348
:
124
8
.
44.
Zou
W
,
Wolchok
JD
,
Chen
L
. 
PD-L1 (B7-H1) and PD-1 pathway blockade for cancer therapy: mechanisms, response biomarkers, and combinations
.
Sci Transl Med
2016
;
8
:
324r
8r
.
45.
Hou
Y
,
Guo
H
,
Cao
C
,
Li
X
,
Hu
B
,
Zhu
P
. 
Single-cell triple omics sequencing reveals genetic, epigenetic, and transcriptomic heterogeneity in hepatocellular carcinomas
.
Cell Res
2016
;
26
:
304
19
.
46.
Robinson
DR
,
Wu
Y
,
Lonigro
RJ
,
Vats
P
,
Cobain
E
,
Everett
J
. 
Integrative clinical genomics of metastatic cancer
.
Nature
2017
;
548
:
297
303
.
47.
Zhang
J
,
Fujimoto
J
,
Zhang
J
,
Wedge
DC
,
Song
X
,
Zhang
J
. 
Intratumor heterogeneity in localized lung adenocarcinomas delineated by multiregion sequencing
.
Science
2014
;
346
:
256
9
.
48.
Pritchard
CC
,
Mateo
J
,
Walsh
MF
,
De Sarkar
N
,
Abida
W
,
Beltran
H
. 
Inherited DNA-repair gene mutations in men with metastatic prostate cancer
.
N Engl J Med
2016
;
375
:
443
53
.
49.
Reck
M
,
Rodríguez-Abreu
D
,
Robinson
AG
,
Hui
R
,
Csőszi
T
,
Fülöp
A
. 
Pembrolizumab versus chemotherapy for PD-L1–positive non–small-cell lung cancer
.
N Engl J Med
2016
;
375
:
1823
33
.
50.
Forde
PM
,
Chaft
JE
,
Smith
KN
,
Anagnostou
V
,
Cottrell
TR
,
Hellmann
MD
. 
Neoadjuvant PD-1 blockade in resectable lung cancer
.
N Engl J Med
2018
;
378
:
1976
86
.
51.
Gibney
GT
,
Weiner
LM
,
Atkins
MB
. 
Predictive biomarkers for checkpoint inhibitor-based immunotherapy
.
Lancet Oncol
2016
;
17
:
e542
51
.
52.
Supek
F
,
Lehner
B
. 
Differential DNA mismatch repair underlies mutation rate variation across the human genome
.
Nature
2015
;
521
:
81
4
.