Expression and Gene Regulation Network of SOX9 in Skin Metastasis of Gastric Carcinoma Derived From One Special Patient


 BackgroundAn exceptional case of a patient with advanced gastric cancer is presented in this study, treated with multiple chemotherapy, radiotherapy and immune regimens, who exhibited regression of one metastatic lesion with concomitant progression of the other lesion during a treatment of PD-1 antibody period. MethodsUsing whole-exome sequencing, TBM, MSI/MSS, PD-1 positive or negative respond measurement, gene function, tumor expression, clinical stage, survival curve, pathopoiesia gene prediction, neoantigen score and immunogenomic approaches. ResultsWe found that SOX9 might mainly participate into the response process of PD-1 antibody in the right skin metastasis. Then, we used sequencing data from the Cancer Genome Atlas database and Gene Expression Omnibus, analyzing SOX9 expression and gene regulation networks in gastric carcinoma (GC). Expression and CNVs were analyzed using Oncomine and Gene Expression Profiling Interactive Analysis tools, at the same time, SOX9 alterations and related pathway were identified using cBioPortal. LinkedOmics and GeneMANIA were used to identify differential mRNA expression with SOX9 and to analyze Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathways. Gene enrichment analysis examined target networks of kinases, miRNAs, transcription factors and PPIs. The results show that SOX9 is overexpressed and the alteration type of SXO9 is mutation in GC. Expression of this gene coding protein is associated with biological interaction networks involving the cell dedifferentiation and WNT signaling process. Functional network analysis indicates that SOX9 mRNA level regulates the RNA process, DNA replication and cell cycle via pathways involving several cancer-related kinases, miRNAs and transcription factors, like casein kinase 2 alpha 1, cyclin dependent kinase 2, Mir296 and Mir214. ConclusionsOur results finally demonstrate that bioinformation analysis efficiently reveals online information of SOX9 expression and potential regulatory networks in GC, laying a foresight for further study of the role of SOX9 and a new PD-1 treatment predictor in gastric carcinogenesis.


Introduction
PD-1 antibody is a new kind of immunotherapy, which competitively inhibits the PD-1/PD-L1 pathway and potentially augment endogenous antitumor responses. On the base of clinical trial KEYNOTE-059, whose was a multicenter, nonrandomized, open-label trial, pembrolizumab was granted its approval by FDA for the treatment of recurrent locally advanced or metastatic gastric cancer patients in September 2017 [1]. On the other hand, in 2019 CSCO, third-treatment of PD-1 blockade for metastatic gastric cancer has been recommended from III to II in guideline [2]. The evidence for the latest change of guideline is supported by ATTRACTION trial, which indicates that no matter PD-L1 positive or negative, terminal gastric cancer (GC) patients all can have bene ts from treatment of Nivolumab and KEYNOTE-059 trials, which suggests that treatment of Pembrolizumab by PD-L1 positive patients could receive better results for those whose gastric cancer cell has tolerated the treatment of rst-line and second-line or third-line chemotherapy.
Advanced gastric cancer commonly indicates poor prognosis and causes metastasis to the peritoneal cavity, lymph nodes, liver, and lung. However, metastasizing to patients' skin is a very rare site with an incidence rate of 1%, which is generally a late clinical presentation indicating systemic dissemination with short lifetime [3]. Generally, cutaneous metastasis from gastric cancer indicates inoperable state, and few effectively therapeutic methods in these cases have so far been reported.
SOX9 is an essential transcription factor, which can regulate gender and bone development. In 1994, Wagner et al. rst found the presence of SOX9 in patients with skeletal development defects [4]. The SOX family consists of 20 members, which are involved in the embryonic development, gender determination, human genetic syndrome and malignancy process [5]. The SOX family can be divided into 9 subgroups, based on homology part of the HMG domain and other structural characteristics [6]. Among these factors, SOX9 plays a key role in embryogenesis, male sexual development, organ development, chondrocyte differentiation and stem cell characteristics [7][8]. At the same time, SOX9 has in uence on neurodevelopment and pancreatic development [9][10].
In our studies with this special patient samples, in which a patient diagnosed with gastric cancer who developed resistance to chemotherapy was found sensitive to treatment with the addition of PD-1 monoclonal antibody [11][12]. Importantly we have obtained cutaneous metastatic tissue exhibiting different sensitivity to PD-1 antibody treatment for immunological and histological analysis. Thus, we explore different genes' expression and mutations from this patient and analysis the possibility reasons through Genome Atlas (TCGA) and various public databases. Using multi-dimensional analysis methods, we found that genomic alterations and functional networks was very likely related to this process. Therefore, our results might potentially reveal new targets and strategies for advanced GC diagnosis and treatment in the future.

Materials And Methods
Data resource: A 59-year-old man was referred to our hospital for frequent micturition in September 2017. During routine checkups, an abdominal CT scan found multiple enlarged lymph nodes in the porta hepatis and perigastric regions with slightly enhanced linings of the antrum. Subsequent gastroscope and pathology con rmed diagnosis of poor differentiated adenocarcinoma in the anterior wall of the gastric angular notch with suspected vascular invasion. Immunohistochemical examination showed that the tumor biopsy sample was positive for Ki-67(70%), carcinoembryonic antigen (CEA), CA221, CK(AE1/AE3), caudal type homeobox transcription factor 2 (CDX-2), but negative for P53, CK20, CD45(LCA), S-100. His serum CEA level was 10.4ng/ml (Figure 2A). The gastric cancer was classed T4N2-3M0 according to the AJCC classi cation. After receiving three rounds of FOLFOX, (Oxaliplatin + Calcium Folinatc + Fluorouracil) adjuvant chemotherapy, and 1 round of fractionated radiotherapy(PTV D95 45Gy/25F and PGTV D95 50Gy/25F) abdominal enhanced CT and liver MRI checkups suggested partial remission, with signi cant reduction in lymph node size in the porta hepatis and peri-gastric regions, but no signi cant changes in the primary tumor site. Serum CEA level was reduced to 1.6ng/ml. The patient received distal gastrectomy in Dec 2017. Postoperative histopathological examination indicated residual cancerous mucous glands with positive staining of CK, P53 and Ki67. Tumor cells positive for P53 was found scattered in the muscularis propria and in the 8th group of lymph nodes. Tumor regression at the primary site was over 90%. Immunohistochemical staining were positive for Ki-67, CK18, P53, CK(AE1/AE3), caudal type homeobox transcription factor 2 (CDX-2), but negative for EMA stains. After surgery, the patient continued to receive 5 rounds of FOLFOX chemotherapy.
By July 2018, before the 4th round of FOLFOX, abdominal enhanced CT found an enlarged lymph node in the retroperitoneum, at the level of the left renal vein. Serum CEA reached 6ng/ml. Multi-disciplinary team (MDT) discussion suggested disease progression and recommended the addition of radiotherapy (PTV D95 45Gy/25F and PGTV D95 50Gy/25F). Approximately 4 months later, in Oct 2018, physical examination of the patient found multiple, reddish-colored lesions in the left and right alar area, which were identical to cutaneous metastasis ( Figure 2C). Two of the largest lesions measured up to 2.7*1.5cm (left) and 0.8*0.5cm (right). Moreover, palpable enlarged lymph nodes were found in the left and right underarm and supraclavicular areas. While abdominal CT showed shrinkage of the retroperitoneum lymph node, CEA level continued to rise. The overall examination results indicated further progression and resistance to FOLFOX treatment. We next treated the patient with 3 rounds of paclitaxel and S-1(D1/D8) plus Apatinib. The regimen was well tolerated but CEA levels rose to 121ng/ml, with no observable changes in the cutaneous lesions.
In March 2019, MDT discussion recommended chemotherapy with paclitaxel and S-1 in combination with PD-1 antibody. Fortunately, after the rst round of treatment, the size of the cutaneous metastatic lesion in the left alar area reduced to 1.5*0.8cm while the other lesion remained the same. Serum CEA level decreased to 69.6ng/ml. The two lesions were subsequently resected by surgery, and the regimen was continued. The two cutaneous lesions were subjected to further pathology analysis. (Figure 2J and 2K) Other data were collected from TCGA and other public databases.
Whole-exome sequencing, TBM, MSI and PD-1 respond measurement The above four parts measurement were conducted by the company of biomed union, in Hangzhou, Zhejiang province. (Figure 1 and 2) GeneMANIA analysis GeneMANIA (http://www.genemania.org) is a exible, user-friendly web interface for constructing proteinprotein interaction (PPI) network, generating hypotheses about gene function, analyzing gene lists and prioritizing genes for functional assays [13]. The website can set the source of the edge of the network, and it features several bioinformatics methods: physical interaction, gene co-expression, gene colocation, gene enrichment analysis and website prediction. We used GeneMANIA to visualize the gene networks and predict function of genes from the results of Whole-exome sequencing.
We then used GeneMANIA to predict the most possible relationship with SOX9 on the way of coexpression. The performed GO and KEGG pathway enrichment analyses of the top 50 genes (the highest score) were committed by the Cytoscape app. The GO annotation had three parts: cellular component (CC), biological process (BP), and molecular function (MF).
KOBAS analysis KOBAS 2.0 is an update of KOBAS (KEGG Orthology Based Annotation System). It can identify statistically signi cantly enriched pathways, human diseases, and functional terms for an input set of genes using biological knowledge from well-known pathway databases, disease databases, and Gene Ontology [14]. KOBAS was used to analysis KEGG pathway enrichment of the both side metastasis genes.
Polyphen2, Provean/SIFT, Mutation Assessor, VEST3 and FATHMM analysis PolyPhen-2 is a new development of the PolyPhen tool for annotating coding nonsynonymous SNPs. Some of the highlights of the new version are: High quality multiple sequence alignment pipeline; Probabilistic classi er based on machine-learning method; Optimized for high-throughput analysis of the next-generation sequencing data [15].
PROVEAN was developed to predict whether a protein sequence variation affects protein function.
PROVEAN is able to provide predictions for any type of protein sequence variations including the following. Single or multiple amino acid substitutions. Single or multiple amino acid insertions. Single or multiple amino acid deletions [16].
Mutation Assessor server predicts the functional impact of amino-acid substitutions in proteins, such as mutations discovered in cancer or missense polymorphisms. The functional impact is assessed based on evolutionary conservation of the affected amino acid in protein homologs. The method has been validated on a large set (60k) of disease associated (OMIM) and polymorphic variants [17].
VEST3 is a method that predicts the functional signi cance of somatic missense mutations observed in the genomes of cancer cells, allowing mutations to be prioritized in subsequent functional studies, based on the probability that they give the cells a selective survival advantage [18].
FATHMM is capable of predicting the functional effects of protein missense mutations by combining sequence conservation within hidden Markov models (HMMs), representing the alignment of homologous sequences and conserved protein domains, with "pathogenicity weights", representing the overall tolerance of the protein/domain to mutations [19].
The above four prediction model were utilized to analysis the morbigenous score of the two metastasis sides genes and conduct mutual veri cation.
NetChop-3.1, NetCTL-1.2, NetMHC-4.0 and NetMHCpan-4.0 analysis The NetChop server produces neural network predictions for cleavage sites of the human proteasome, which has been trained on human data only and will therefore presumably have better performance for prediction of the cleavage sites of the human proteasome. However, since the proteasome structure is quite conserved, the server is able to produce reliable predictions for at least the other mammalian proteasomes. The Netchop 3.0 version has two different network methods that can be used for prediction. C-term 3.0 and 20S 3.0. C-term 3.0 network is trained with a database consisting of 1260 publicly available MHC class I ligands (using only C-terminal cleavage site of the ligands). 20S network is trained with in vitro degradation data published in Toes, et al. and Emmerich et al. C-term 3.0 network performs best in predicting the boundaries of CTL epitopes [20]. NetCTL 1.2 conducts prediction of CTL epitopes in protein sequences. The version 1.2 expands the MHC class I binding prediction to 12 MHC supertypes including the supertypes A26 and B39 and the prediction of proteome cleavage has been updated and is now identical to the predictions obtained by the NetChop-3.0 server. The new version has been trained on a set of 886 known MHC class I ligands. The method integrates prediction of peptide MHC class I binding, proteasomal C terminal cleavage and TAP transport e ciency. MHC class I binding and proteasomal cleavage is performed using arti cial neural networks. TAP transport e ciency is predicted using weight matrix and the MHC peptide binding is predicted using neural networks trained as described for the NetMHC server. The proteasome cleavage event is predicted using the version of the NetChop neural networks trained on C terminals of known CTL epitopes as describe for the NetChop-3.0 server. The TAP transport e ciency is predicted using the weight matrixbased method describe by Peters et al. [21].
NetMHC-4.0 has been mainly trained for 81 different Human MHC alleles including HLA-A, -B, -C and -E. Furthermore, predictions for 41 animals (Monkey, Cattle, Pig, and Mouse) alleles are available. Predictions can be conducted for peptides of any length now [22].
The above four prediction model were utilized to analysis the neoantigen producing score of the two metastasis sides genes, which indicated immune reaction and conduct mutual veri cation.

GEPIA analysis
GEPIA is a newly developed interactive web server for analyzing the RNA sequencing expression data of 9,736 tumors and 8,587 normal samples from the TCGA and the GTEx projects, using a standard processing pipeline. GEPIA provides customizable functions such as tumor/normal differential expression analysis, pro ling according to cancer types or pathological stages, patient survival analysis, similar gene detection, correlation analysis and dimensionality reduction analysis. This tool is developed by Zefang Tang, Chenwei Li and Boxi Kang of Zhang Lab, Peking University [24]. The RNA-Seq datasets GEPIA used is based on the UCSC Xena project (http://xena.ucsc.edu), which are computed by a standard pipeline. We used this tool to present the clinical expression, stage and survival curve of the two skin metastasis genes.
Timer analysis TIMER web is a comprehensive resource for systematical analysis of immune in ltrates across diverse cancer types. The abundances of six immune in ltrates (B cells, CD4 + T cells, CD8 + T cells, Neutrophils, Macrophages, and Dendritic cells) are estimated by TIMER algorithm [25]. We also utilized this web to explore the immune analysis of skin metastasis in this patient, including relation, survival and SCNA levels.

Oncomine analysis
The mRNA expression and DNA copy number of SOX9 in GC were analyzed through the Oncomine 4.5 database and it is the world's largest oncogene chip database, which integrated data mining platform, containing 715 gene expression data sets and data from 86,733 cancer tissues and normal tissues [26].
The analysis of SOX9 drew on a series of GC studies, including Derrico GC, Cho GC, Chen GC, Cui GC, Wang GC, Deng GC and TCGA GC studies. UALCAN analysis UALCAN uses TCGA level 3 RNA-seq and clinical data from 31 cancer types, which is an interactive webportal to perform in-depth analyses of TCGA gene expression data [27]. We used UALCAN to analysis the subgroups of patients with GC, strati ed based on gender, age and other criteria of SOX9 expression.

c-BioPortal analysis
The cBio Cancer Genomics Portal (http://cbioportal.org) is an open-access resource for interactive exploration of multidimensional cancer genomics data sets, currently containing 225 cancer studies [28]. c-BioPortal was utilized to analyze SOX9 alterations in the TCGA GC sample. The search content included mutation, CNVs, and mRNA expression. The tab OncoPrint displayed an overview of genetic alterations per sample in SOX9. The tab Pathway visualized the match score of different pathways related with SOX9 derived from public pathway databases [29].

LinkedOmics analysis
The LinkedOmics database (http://www.linkedomics.org/login.php) is a Web-based platform for analyzing 32 TCGA cancer-associated multi-dimensional datasets [30]. The LinkFinder module of LinkedOmics was used to study genes differentially expressed in correlation with SOX9 in the TCGA GC cohort. Results were analyzed statistically using Pearson's correlation coe cient. The LinkFinder also created statistical plots for individual genes. All results were graphically presented in volcano plots, heat maps or scatter plots. Meanwhile, the GSEA method identi ed as being enriched in gastric cancer: kinase, mi-RNA, transcription factor and PPI.

Statistical analysis
The obtained results were expressed as means ± S.D. Statistical analyses were performed using one-way ANOVA with Bonferroni post hoc test or the Student's t-test. The results were considered to be statistically significant at a two-tailed P value of less than 0.05.

Results
Patient's primary data analysis Phylogenetic analysis of somatic mutations in primary and metastasis tumors We performed whole-exome sequencing of the resected tumors to identify somatic mutations in the primary tumor and both skin metastases. Of all samples, we detected that the number of gene mutation in both skin metastases is surpassing than primary tumor. Among these dates, mutation of gene FAT1 was found in both primary tumor and right metastasis tumor. Meanwhile, gene EGR mutated at two DNA points, 923 and 929 points respectively, in both skin metastases. To estimate immune respond in all samples, we used TMB to measure and the data show that the highest number of TMB is right skin metastasis tumor and the primary tumor was the lowest one ( Figure 2D). However, all tumor samples were MSS and PD-1 negative respond, which indicates poor respond to PD-1 antibody drug ( Figure 2J and 2M). In all, both metastases were genetically more heterogeneous and harbored more mutations than primary tumor. (Figure 2) Tumor cell data analysis Phylogenetic analysis of gastric cancer progression between skin metastasis tumors The patient produced different response to PD-1 antibody between left and right skin metastasis tumor during the clinical treatment period. In order to point out the reason, we analysis the whole-exome sequencing results. Gene function, tumor expression, clinical stage, survival curve, pathopoiesia gene prediction and neoantigen score were conducted to study these genes. First, we contrast the gene data of two skin metastases through the TCGA data analysis of GEPIA. Just as showed in Figure 3(A, B, C and D), ve genes expression (CDC27, EP300, IGF2, KAT6A, SOX9) contain clinical signi cance, compared with normal tissues, in gastric cancer tissues, drawing from TCGA database. In addition, ERG and KMT2C protein expression have clinical signi cance during clinical stages. ERG and KMT2C are highly expressing in stage IV and III respectively ( Figure 3B). While in survival aspect, only KEAP1 and LRP1B could be received as indicator to predict survival condition ( Figure 3D).
Then, in gene function part, we use gene interaction module at the site of GeneMANIA and visualize their network as in Figure 3(E and H). We nd that EP300 and KAT6A, KAT6A and IGF2 are co-expressing in gastric cancer. In addition, the analysis of signi cantly enriched gene ontology (GO) and KEGG terms indicated that genes of left skin metastasis tumor are harbored with cell proliferation, migration and angiogenesis. While the right-side tumor cell genes reveal immune relative pathway, excepted with tumor growth and migration. (Figure 3F, G, I and J) At last, to check out the potential genes, which might mainly lead to gastric cancer, we use ve predictive modules, (Polyphen2, Provean/SIFT, Mutation Assessor, VEST3 and FATHMM) to compare the two sides tumors. Just as showed in Figure 3K-P, morbigenous scores of right skin metastasis tumor are lower than the left side tumor, which indicates that the left skin metastases might display more malignance character than the other side tumor. Among these data, PTCH1 and FGFR3 gens are more possibility leading to cancer. While predicting neoantigen producing, which means T cell response, the score of the right-side tumor is higher than the left skin metastasis, through the module of NetChop 3.1, NetCTL 1.2, NetMHC 4.0 and NetMHCpan 4.0, presenting that the possibility of generated neoantigen may be checked out easily in right skin metastasis tumor ( Figure 3Q, R, S and T). Gene mutation of EP300 is the highest potential mutation to produce neoantigen.

Immune cell data analysis Phylogenetic analysis of immune cell between skin metastasis tumors
In order to experience the reason why the two side tumors produce different response to PD-1 antibody. We utilize six kinds of immune cells to compare both side tumors at relation, survival and SCNA levels by the analysis of TIMER. In Figure 4A and 4B, we nd that the left side metastasis tumor is more related to B cells and neutrophil cells, while CD4 and CD8 cells are more possibly in ltrating in the right skin metastasis tissue, and SOX9 is more related with CD4 and CD8 in ltrated. However, survivorship curve is no difference in both tumors ( Figure 4E and 4F). Moreover, we analysis the SCNA data of both sides and the result show that there are more arm-level deletion and arm-level gain occurrence rate in right skin metastasis than the left one ( Figure 4C and 4D). We still nd that arm-lever deletion or gain of SOX9 and SNCAIP are more frequent with immune cells SCNA. At last, we compare the relationship between both tumors with immune reaction and more pertinence could be detected in the right skin metastasis tumor. We detect that SOX9 and SNCAIP of right-side tumor and FGFR3 of left side tumor have signi cance of immune respond.
Expression and gene regulation analysis of SOX9 SOX9 expression in GC Through the above two parts analysis of cancer cells and immune cell, SOX9 was found to mainly participate into the response of PD-1 antibody in the right skin metastasis and could be utilized as a further predictor for PD-1 antibody treatment. In order to estimate the possibility mechanism of SOX9. We initially evaluated SOX9 transcription expressions in multiple GC studies from TCGA and the Gene Expression Omnibus (GEO). Data of the Oncomine 4.5 database revealed that mRNA expression and DNA copy number variation (CNV) of SOX9 were signi cantly higher in GC tissues than in normal tissues (p 0.01) ( Figure 5A-5P). Further sub-group analysis of multiple clinic pathological features of gastric carcinoma (GC) samples in the TCGA consistently showed high transcription of SOX9. The transcription level of SOX9 was signi cantly higher in GC patients than healthy people in subgroup analyses based on gender, age, ethnicity, disease stages and tumor grade ( Figure 5Q-5X). Therefore, SOX9 expression may serve as a potential diagnostic marker in GC.

Frequency and type of SOX9 alterations in GC
We further utilized the cBioPortal wed to assess the types and frequency of SOX9 alterations in GC based on sequencing data from GC patients. SOX9 was altered in 5.68% (440 patients) GC patients ( Figure 5Y). These alterations contained mRNA deep deletion in 2 cases (0.45%), ampli cation in 8 cases (1.82%), mutation in 14 cases (3.18%), and multiple alterations in 1 case (0.23%). Hence, mutation is the most common type of SOX9 CNV in GC.

Biological interaction network of SOX9 protein level in GC
We sequentially wanted to determine the biological interaction network of SOX9 protein level in GC. So, we used the co-location, co-expression and pathway functions in GeneMANIA to show SOX9 relative genes coding proteins ( Figure 6E). The relative genes of SOX9 with the highest score were SLC12A2(0.00059), KRT23(0.00054) and CEACAM1(0.00049). The 50 highest scored genes were analyzed by gene ontology (GO) terms and the outcomes indicated that these genes encode proteins localized mainly to the circulation and cell cytoplasm ( Figure 6A). These proteins are primarily involved in negative regulation of cell differentiation and cell development, while they also serve as cell surface signal pathway and WNT cell-cell signaling ( Figure 6B and 6C). Meanwhile, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis showed enrichment in bio-rhythmic process and receptor tyrosine kinase signaling ( Figure 6D). Therefore, the biological interaction network of SOX9 alterations is involved in the cell dedifferentiation and WNT signaling process.
GO and KEGG pathway analyses of co-expression genes at the SOX9 mRNA level in GC The Function module of LinkedOmics was used to analyze mRNA sequencing data from GC patients in the TCGA. As shown in the volcano plot ( Figure 7A), 8,451 genes (dark red dots) showed signi cant positive correlations with SOX9, whereas 8,305 genes (dark green dots) showed signi cant negative correlations (false discovery rate [FDR] < 0.01). The 50 signi cant gene sets positively and negatively correlated with SOX9 as shown in the heat map ( Figure 5B and 5C). This outcome suggested a widespread impact of SOX9 on the transcriptome. SOX9 expression showed a strong positive association with expression of LLGL2 (positive rank #1, Pearson correlation=0.601, p= 4.10E-42), and ERBB3 (positive rank #2, Pearson correlation=0.59, p= 2.28E-40), which re ected changes in asymmetric cell division, epithelial cell polarity, cell migration and epidermal growth process. Signi cant GO term analysis by gene set enrichment analysis (GSEA) displayed that genes differentially expressed in correlation with SOX9 were located mainly in the chromosomal region, nucleolar part and spliceosomal complex, where they participated primarily in catalytic activity on RNA and ATPase activity. They also act as ncRNA processing, epidermis development and DNA replication ( Figure 7D, F, H). KEGG pathway analysis showed enrichment in ribosome biogenesis in eukaryotes, cell cycle and splicesome ( Figure 7J). Meanwhile, GSEA analysis of the above four parts were mainly emphasized on corni ed envelope, chromosome segregation, DNA secondary structure binding and maturity onset diabetes of the young ( Figure 7E, G, I and K). In addition, pathway of the gene set enriched for SOX9 was involved with the WNT pathway, which was the highest match score in the c-BioPortal web ( Figure 9B). The alterations of relative genes in this pathway were mRNA deep deletion, ampli cation, mutation, fusion and multiple alterations ( Figure 9A and 9C).

Discussion
Cutaneous manifestation of gastric cancer is very rare and associated with a poor prognosis and widespread metastasis in patient's body, which indicates short lifetime for the patient [31]. Furthermore, skin metastasis may be the rst clinically sign for underlying systemic malignancy and is di cult to be detected in early period. Therefore, immediate clari cation in case of uncertainty is recommended and need to be conducted. In this patient, early screen could be done in combination with serum CEA and CT scan, which would be bene cial for treatment course. Meanwhile, de nitive diagnosis can be made through the biopsy pathology. Presently, Treatment for cutaneous metastasis of gastric cancer include surgical resection, radiotherapy, chemotherapy and others [32]. Nonetheless, surgical resection of skin metastases arising from gastrointestinal cancers is mostly undertaken as palliative treatment to improve the patient's quality of life. Unfortunately, clinical effect of radiotherapy and chemotherapy is still uncertain and need to be further explored. In this patient, PD-1 treatment conspicuously remits skin metastasis and declines the level of patient's serum CEA. However, speci c reason is still out of recognition.
Immunotherapy bring new strategy for metastatic gastric cancer after the rst granted PD-1 inhibitor in 2014 year [33]. However, there have been no effective biomarkers for predicting treatment outcome of PD-1 inhibitors, although expression of PD-1, MSI-H and high TMB have been utilized in clinic treatment. In reality, ORR is still low, during 10%-40% and 30%, in PD-1 positive and TMB-H patients respectively and this data just overpasses 50%, increasing to 53%, in MSI-H cases [34]. In this case, we research both of the patient's skin metastasis samples through whole exome sequencing. Surprising to us, both samples are PD-1 negative, MSS and low TMB, which are differ to those three guideline biomarkers. These three predictors indicate poor therapeutic effect, while the left cutaneous metastasis fade away after PD-1 inhibitor treatment.
Therefore, new markers for advanced gastric cancer early diagnose and PD-1 antibody treatment are needed to improve diagnosis and cure. Analysis of gene sequencing data, gene function, tumor expression, clinical stage, survival curve, pathopoiesia gene prediction and neoantigen score from this special report.
Firstly, we analysis the reason why two skin metastasis tumors committee to different respond to PD-1 antibody. So, we compare the TMB of both metastases and use pathogenicity prediction, neoantigen, function relation, immune cells in ltration, survival curve, SCNA and immune correlation detection. The results show that the TMB of right skin metastasis is higher than the other side, although both sides' TMB are lower than the standard. On the other side, producing of neoantigen, in ltrating of CD4, CD8 and immune correlation are more possibly found in the right metastasis tumor in this patient, which complete respond to PD-1 antibody. While the left metastasis tumor has more oncogenicity. However, the effective indicators of PD-1 antibody in right side tumor are none positive, like PD-1 expression, TMB and MSI. Therefore, we search the recent research paper and gure out that these above indicators are also showing their shortages, which means that the clinical effective results of recent PD-1 indicators are lower than 50% [35][36]. The real practical effectors are still going to be researched in the future. And in our case results, SOX9 gene may be the next potential target. In summary, two skin metastases contain a part of same gene mutation and also own themselves' different genes, which are showing the progress of heterogeneity in gastric cancer.
Secondly, although the patient received surgery and postoperative chemotherapy, the metastasis of retroperitoneal lymph node came out two months later. After radiotherapy, two sides skin metastasis indicate the poor survival of this patient. In order to gure out why skin metastasis happened so quickly, we utilize the above methods to explain the mechanism between primary tumor and metastasis tumor.
The results show that metastasis tumor display higher expression, stage III or IV, poor survival, function of angiogenesis, migration and more relation with metastasis proteins. At the same time, we suggest that SOX9 mutation in the above result has statistical signi cance and deserves further clinical validation as a potential diagnostic and prognostic marker for advanced GC of PD-1 antibody treatment. So, we continue to conclude the gene, mRNA and protein levels of SOX9 in this paper.
CNVs contains genomic insertions, disrupting genes and altering genetic content, leading to phenotypic differences at cell functions [30]. Our study found that the copy number of SOX9 was majorly increased in GC tissues and the major type of SOX9 alteration was mutation (3.18%), which was associated with shorter survival. We speculated that altered SOX9 expression and SOX9 dysfunction in GC might result from alterations in chromosomal structure. Since SOX9 played several important physiological functions, its alteration might cause changes in various downstream signaling pathways. Indeed, related functional researches were involved in negative regulation of cell differentiation, cell development, cell surface signal pathway, WNT signaling, bio-rhythmic process and receptor tyrosine kinase signaling. Thus, the network of SOX9 alterations in protein level was involved in the cell dedifferentiation and WNT signaling process.
Enrichment analysis of target gene sets using mRNA sequencing data can help reveal signi cant transcriptional change and functional dysregulation. Our results suggested that the functional network of SOX9 participated primarily in the catalytic activity on RNA and ATPase activity, ncRNA processing, epidermis development and DNA replication, ribosome biogenesis in eukaryotes, cell cycle and splicesome. It is crucial to understand how alteration in a gene transcription level can lead to major dysfunction and even to cancer such as GC and respond to drug treatment.
Genomic instability and mutagenesis are fundamental characteristics of cancer cells and their associated kinases, miRNAs, proteins, signaling pathways can promote the progression of cancers as well [37]. We found that SOX9 in GC was associated with a network of kinase including CSNK2A1, CDK2 and CDK1. These kinases regulated tumor cell migration, metastasis and cell disfunction. Therefore, in GC, SOX9 may regulate tumor in ltration and metastasis via these kinases interaction.
Meanwhile, several miRNAs are associated with SOX9. These non-coding RNAs, normally involved in post-transcriptional regulation of gene expression, can contribute to human carcinogenesis [38]. The particular miRNAs, (MIR-296, MIR-214, MIR-432) in our study had been linked to tumor proliferation, apoptosis, invasion, metastasis and drug resistance [39][40][41]. In fact, hsa-mir-296 was correlated with metastatic disease in gastric cancer. In addition, mir-214 have the capacity of promoting apoptosis and participating in the drug resistance of Cisplatin in GC. So, disfunction of these miRNAs were consistent with the overexpression phenotype of SOX9 in GC.
Last but not least, pathway analysis of SOX9 showed that WNT pathway was the most possible associated with the functional network of SOX9. The WNT signaling pathway was a complex network of protein interactions, which could be commonly found in the embryonic development and cancer progression [42]. As previous research reported, SOX9 had participated into the WNT signal in many cancers, which was consistence with our results [43].

Conclusions
In conclusion, this patient of gastric cancer experiences progression of metastasis and partly responds to PD-1 antibody. Meanwhile, this special case provides evidence of divergent tumor gene transcriptome, immune cell in ltering and bioinformation analysis [36]. If this phenomenon proves generalizable, then the inter-site heterogeneity described here bespeaks a profound clinical challenge for the clinical treatment. Given the data presented in this study, it will be essential to understand not only how to therapeutically target genomic heterogeneity between and within metastases but also how to successfully mobilize an anti-tumor immune response able to control all metastases in advanced cancers. In addition, this study further on provides multi-level evidence for the gene function of SOX9 in GC and its potential as a predictor for PD-1 antibody treatment. On the other hand, this study uses online tools based on the bioinformatics theories to perform target gene analyses on tumor data from public databases (TCGA, GEO and so on). Although, these data are collected from global research, the TCGA database has limitations. One is that the TCGA GC samples contain different measurement standard.
The possible error might be found in these data. Another limitation is that the GC samples analysis just online, being short of experimental validation. Therefore, our results should be veri ed in clinical samples showing su cient coverage. These questions should be addressed in follow-up studies using molecular biology techniques.

Declarations
Ethics approval and consent to participate Gastric cancer tissue and skin tissue were collected from the surgery of one patient at Department of Gastrointestinal Surgery, the Second A liated Hospital of Zhejiang University, School of Medicine. The study was conducted in accordance with the Declaration of Helsinki. Informed consent was obtained from the patient and legally authorized representatives, and the approval of all tissues use was obtained from the Ethics Committee of the Second A liated Hospital of Zhejiang University, School of Medicine (approval No. applying).

Consent for publication
Not applicable.

Availability of data and materials
All data generated or analysed during this study are included in this published article.
Competing interests Figure 1 The pathway of this paper. This study contains three parts: rst, a special skin metastasis of gastric cancer patient report and gene result are presented; second, cancer cell and leukomonocyte are analyzed from the above result; third, hub gene SOX9 is studied from 4 aspects.     The Protein, miRNA and transcription factor-target networks of SOX9 in GC (LinkedOmics). The target factor of SOX9 in LinkedOmics.