A 26-year-old female patient was admitted to the First Affiliated Anhui Medical Hospital in December 2016 having presented with GC-related symptoms. The patient felt uncomfortable in the upper left quadrant of her abdomen and had repeated black melena for one month. This symptom recurred and became more serious for 2 weeks. The patient felt full and uncomfortable in the upper left quadrant of her abdomen repeatedly following a meal, one month prior to the development of black melena. She also had nausea, although there was no hematemesis or hematochezia. The abdominal discomfort worsened and the patient received treatment in Luan Jinkai Hospital. Thereafter, the patient visited the First Affiliated Hospital of Anhui Medical University where the present study was conducted. The gastroscopy of the patient indicated GC with obstruction, while pathology indicated antrum considered as carcinoma mucocellulare. Pathological examination of the samples was performed using hematoxylin and eosin staining. Briefly, tissues were sliced with a dimension of 1.5×1.5×0.2-0.3 cm and soaked in 40°C warm water. The sections were deparaffinized in xylene and were gradually hydrated through graded ethanol (100, 95, 80 and 70%) each for 3 min at 25°C. The sections were stained in hematoxylin solution for 10 min at 25°C and differentiated in 1% hydrochloric alcohol. Then the slides were rinsed with tap water and distilled water until the nuclei became blue, and dehydrated in 95% ethanol. The sections were counterstained in 1% eosin solution for 5 min at 25°C, washed with 70% ethanol twice and absolute ethanol, and were cleared in two washes of xylene. The sections were mounted with neutral balsam and observed under a light microscope magnification, ×400 (×40 objective and ×10 ocular magnification; BX-42; Olympus). An abdominal computer tomography scan + enhancement showed that both armpits had multiple small lymph nodes and the gastric antrum had thickened walls along with multiple swollen lymph nodes in the surrounding area. The patient was diagnosed with a T3N3bM0 (advanced stage IIIC) GC using AJCC TNM staging system 7th edition (14).

The patient underwent a gastrostomy for tumor removal. Tumor tissue and adjacent normal tissue ~5-cm away from the tumor were resected and preserved in liquid nitrogen until use.

Genomic DNA was isolated from the resected tissues using a genomic DNA isolation kit (Bioo Scientific; http://www.biooscientific.com) according to the manufacturer's recommendations. Extracted DNA was quantified using a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies; Thermo Fisher Scientific, Inc.) and the integrity was assessed using 1% agarose gel electrophoresis and visualized in gel imager (Tanon; http://www.biotanon.com). RNA was extracted from the tumor and adjacent normal tissue using TRIzol^® (Invitrogen; Thermo Fisher Scientific, Inc.), and eluted in RNAse-free water. RNA quantity and quality were assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies Inc.). DNA and RNA libraries were prepared according to standard protocols according to Bioo Scientific.

For whole transcriptome sequencing, raw sequencing reads were produced using an Illumina X10 sequencer (Illumina Inc.) with 150-bp paired-end reads according to the manufacturer's instructions. Raw RNA-Sequencing (RNA-Seq) data was trimmed by removing the adaptors and low quality reads using trim_galore v0.4.4 software (15). The data was subsequently aligned to the human reference genome GRCh37.p13 using STAR v2.6.1a software (16) in a 2-pass mapping. Read count in the tumor and adjacent normal tissues was calculated using htseq-count v0.10.0 software (17), and differentially expressed genes were detected using edgeR v3.24.3 software (18). Genes with P≤0.001 and log₂ fold change ≥1 were considered as differentially expressed.

Whole genome sequencing of the prepared libraries was performed using the Illumina X10 sequencer (Illumina Inc.) with 150-bp paired-end reads according to the manufacturer's instructions. The raw fastq data was mapped to the human reference genome (b37) provided by the Broad Institute with the Burrows-Wheeler Aligner v.7.0.12 (19). Read duplications were marked using the Picard tool v.2.10.10; local realignment and base quality score recalibration was performed using Genome Analysis ToolKit (GATK) v. 3.8-0 (20) and the final bam files were used for variant calling.

Germline variants were called using the GATK HaplotyperCaller joint v3.8 software (21) and then filtered using GATK variant quality score recalibration (VQSR). Variants that passed using the VQSR module with coverage ≥6X and supporting reads ≥0 were retained, while the variants showing genotype ‘0/0’ or ‘./.’ in the tumor sample were filtered out. For somatic variant calling, three different types of software were used: GATK HaplotypeCaller joint v3.8, Mutect v1.1.4 and MuTect2 v3.8 (22). In Mutect and MuTect2, all variants that were flagged as ‘PASS’ were kept. Other variants were retained with a coverage ≥5X and supporting ≥3 reads, if they were flagged among a subset of ‘alternative (23) allele in normal’, ‘triallelic site’, ‘possible contamination’, ‘clustered read position’ in Mutect, or ‘germline risk’, ‘alt allele in normal’, ‘t lod fstar’ (tumor does not meet likelihood threshold), ‘str contraction’ (site filtered due to contraction of short tandem repeat region), ‘triallelic site’ (site filtered because > two alternate alleles pass tumor likelihood threshold) in MuTect2. Germline variants were removed from both the Mutect and MuTect2 output data. In all three softwares, variant allele frequency in the tumor sample had to be 3 times higher compared with that in the normal sample, and >5% of the variant reads in the tumor sample was required. Functional annotations for both germline and somatic variants were added to each mutation using the ANNOVAR software v.2018Apr16 (17) and several publicly available databases, such as 1,000 Genome Project (24), the Exome Aggregation Consortium (25), the Genome Aggregation Database (26) and the compiled scores prediction system dbnsfp33a (27), dbSNP (28) and COSMIC (29).

The transition to transversion ratio for germline variants was calculated by dividing total transitions by total transversions. Transitions refer to variations that involve a change from purine to purine, while transversions refer to variations that involve a change from purine to pyrimidine or vice versa (30). The transition to transversion ratio between homologous strands of DNA is generally ~2, and for human, the ratio is ~2.1 (31).

Copy number variants (CNV), tumor purity and ploidy was analyzed using cnv_facets software (https://github.com/dariober/cnv_facets) based on FACETS v0.5.14 (32), and run with the ‘-cval 25 400’ command. Duplicated and deleted segments were retained and annotated using an ENCODE gene symbol in R v3.5.1 (33). iCallSV (https://github.com/rhshah/iCallSV), which applied Delly v0.7.5 prediction method (34) was used to detect structural variations and fusions using default parameters. Inter- and intra-chromosomal fusions >50 kb were retained and only in-frame deletions, duplications and fusions were selected for subsequent analysis.

The 1,000 Genome Project (24), the Exome Aggregation Consortium (25) and the Genome Aggregation Database (26) provided alternative allele frequency data in different populations that ANNOVAR software applied to the present variations. Annotations using a set of algorithm scoring system from dbNSFP v.33a, including FATHMM, LRT, Mutation Assessor, Mutation Taster, M-CAP, Polyphen2, PROVEAN and SIFT (35), were also applied to the variations. Then the common variants were filtered out, which are sites with a minor allele frequency ≥1%, to get the most informative ones. Mutations predicted to be deleterious by at least two algorithms of dbNSFP were considered to be potential pathogenic sites of these informative variants.

To identify genes with a high possibility of involvement in cancer development, the Cancer Gene Census (CGC) (https://cancer.sanger.ac.uk/census) (36) and Cancer Predisposition Genes (CPG) (37) databases were used. In total, 850 candidate genes from cancer-associated databases (723 genes from CGC and 114 genes from CPG with 87 overlapping genes), were investigated. The gene lists were downloaded directly from the two databases with no specific selection criteria.

Gene enrichment was performed using the R package clusterProfiler v3.10.1 software (38), which implements a hypergeometric model to test for gene set overrepresentation relative to a background gene set. Enrichment was analyzed using the functions: ‘Enriched Kyoto Encyclopedia of Genes and Genomes (KEGG)’, ‘gseKEGG’ (Gene Set Enrichment Analysis (GSEA) (39) of KEGG, and a function of clusterProfiler), P<0.05, OrgDb=org.Hs.eg.db (a database for Genome Wide Annotation for Human, referenced by clusterProfiler).

The differentially expressed genes were compared with mutated genes identified by whole genome sequencing. To achieve this, genes that were mutated and with differentially expressed transcripts were merged in R v3.5.1 (33).

The RNA-Seq data from The Cancer Genome Atlas stomach adenocarcinoma (TCGA STAD) was downloaded from TCGAbiolinks online tool (https://github.com/BioinformaticsFMRP/TCGAbiolinks). To identify the most significantly differentially expressed genes in GC, the results from three RNA-Seq differential expression analysis tools, including limma v3.42.2 (40), DEseq2 v1.26.0 (41) and edgeR v3.28.1 (18), were combined and the overlapping genes were separated for further analysis. For all three tools, P<0.05 and the average log₂ fold change ± standard deviation were used as the cut-off values. Gene Ontology (42) enrichment analysis was performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.7 functional annotation clustering tool to determine the potential functions of the differentially expressed genes. The KEGG pathway analysis was performed to determine the involvement of differentially expressed genes in different biological pathways. The-log10 (P) denotes the enrichment score, indicating the significance of the pathway associations.

KM-plotter (https://kmplot.com/analysis/) is an online database and tool which can be used for the discovery and validation of survival biomarkers (43). It contains cancer gene expression data and survival information from Gene Expression Omnibus (44), European Genome-phenome Archive (EGA) (https://www.ebi.ac.uk/ega/home), and TCGA (https://www.cancer.gov/tcga), and uses each percentile of mRNA expression between the lower (25%) and upper quartiles (75%) of expression as a cut-off point to divide patients into high and low expression groups. The KM-plotter software, which included available transcriptome and survival data (both overall and disease-free survival times), was used to analyze a total of 1,440 patients with GC, based on default parameters.

The proportion and the ploidy of tumor cells in the sample were ~22% and 1.74, respectively. Using whole genome sequencing, the total number of paired reads were 8.13 and 8.34 billion, the mean read coverage was 37.16 and 35.48, while the mean library insert size was 319 and 323 bp for normal and tumor tissue samples, respectively. For the consensus coding sequence region, ~90% of the bases had over 30X coverage. Using whole transcriptome sequencing, 1.04 and 1.11 billion clean reads were obtained with an average input read length of 253 and 261 bp for normal and tumor tissue, respectively.

The RNA-Seq pipeline was used (Fig. S1A) for the whole transcriptome sequencing. Genes with P≤0.001 and absolute log₂ fold change value ≥1 were considered as significantly differentially expressed genes. In total, 766 down- and 765 upregulated genes were obtained (Fig. S1B). The RNA-Seq result was subsequently matched to the whole genome sequencing data for further analysis.

Germline mutation analysis revealed a total of 4,228,339 single nucleotide polymorphisms (SNPs) and Indels with 95.03% of the sites previously reported in the dbSNP database. The transition to transversion ratio was 2.05, indicating that the germline analysis had high confidence, as the empirical value for whole genome analysis is ~2.1 (31). Among the germline variants, 21,405 SNPs/Indels were located in the exonic regions, of which 10,136 were predicted to alter the protein, including 9,504 non-synonymous single nucleotide variations (SNVs), 112 frameshift deletions, 90 frameshift insertions, 173 non-frameshift deletions, 169 non-frameshift insertions, 78 stop-gain, and 10 stop-loss. More than 1% minor allele frequency variants from publicly available databases were removed, as those sites are common in the population. Of the candidate genes from the CGC and CPG databases, 13 SNPs and two Indels were identified in 15 genes (Table I). A deleterious score evaluation of the 15 variants revealed that mutations in DOCK8 (p.D63N), HLF (p.V59A), NKX2-1 (p.G322S), NRG1 (p.S17C), PRDM2 (p.S823T), WNK2 (p.C652R) were deleterious, while mutations in KMT2D (p.V401M), RAD51D (p.A181V) and ROS1 (p.E1605Q) were possibly damaging.

A total of 3 softwares from the Broad Institute were used to identify the somatic variants and 114 variants in the exonic regions were found. This included 45 non-synonymous SNVs, 6 frameshift deletions, 11 frameshift insertions, 17 non-frameshift deletions, 3 non-frameshift insertions, 2 stop-gain and 30 synonymous SNVs. Similar to the germline mutation analysis, common variants were filtered out with publicly available datasets to select the most informative somatic variants. In total, 29 variants were obtained, including 13 mutations which were reported in the Catalogue of Somatic Mutations in Cancer database (Table SI). Of these, 3 protein-truncating variants (BPNT1, p.S144I; FRG1, p.N153D; and TAS2R31, p.L59F) were predicted to be deleterious and one variant (KRTAP5-3, p.C185S) was predicted to be possibly damaging. However, for transcriptional regulation, 3 genes (MUC2, MUC4 and SLC8A2) were found to be differentially expressed.

To identify large-scale variations in the patient, the copy number alterations and complex structural events in the tumor sample were investigated. The structural variant discovery tool, DELLY, reported 701 structural events; however, only in-frame events or fusions over 50 kb between two genes were retained. Ultimately, 20 significant structural variations were reported (Table SII). Another sensitive tool, FACETS, detected 58 large segments, of which 62 were deletions and 112 were duplications found in the coding regions (Table SIII). Combining the results from the two tools, nine genes (CDKN1B, FBXW7, MUC4, CCND3, ETV6, FGFR2, FGFR3, TFEB and KMT2C) were reported in either the CGC or CPG databases. CDKN1B, FBXW7 and MUC4 showed duplications, and there were deletions in CCND3, ETV6, FGFR2, FGFR3 and TFEB. Although KMT2C rearranged with BAGE2, the latter was not reported in the 2 databases. The mRNA expression levels of the genes with structural variations was also investigated and 21 were found to be differentially expressed in the tumor sample (Fig. 1).

RNA-Seq and mutation data were analyzed to characterize the genomic alterations in known signaling pathways. The GO enrichment for somatic mutated genes revealed terms, which included 3 mucin family genes, MUC2, MUC4 and MUC6, and were enriched in ‘maintenance of gastrointestinal epithelium’, ‘O-glycan processing’, and digestive system process, while KRTAP5-3, TCHH, RPTN, POU3F1 were enriched in epidermal cell differentiation pathway (Table SIV). The KEGG enrichment analysis revealed SLC8A2 and SLC25A5 were enriched in calcium signaling and cGMP-PKG signaling pathways (Table SIV). From the RNA-Seq data, GSEA showed that the upregulated genes were enriched in ‘cGMP-PKG signaling’, ‘calcium signaling’, and ‘oxytocin signaling’ pathways, while the downregulated genes were enriched in ‘cell cycle’ pathway, ‘oxidative phosphorylation’ pathway, ‘transcriptional mis-regulation in cancer’ pathway and ‘GC pathway’ (Fig. 2).

To ensure the present study was representative and valuable, the genomic features of the patient was compared with data from TCGA STAD dataset, which was downloaded from TCGAbiolinks, and 329 (136 up- and 193 downregulated genes) differentially expressed genes were identified, which were also found in the RNA-Seq analysis. A total of three tools were used to identify the most significantly differentially expressed genes, and 136 up- and 193 downregulated genes were identified (Fig. S2A), and were subsequently investigated in the tumor tissue sample obtained from the patient with GC and compared with that in the adjacent normal tissue. In total, 27 genes were identified, including five that were up- and 22 that were downregulated (Fig. 3A). GO and KEGG pathway enrichment analysis revealed one (GO:MF) and three (KEGG) significant terms, respectively, including digestion and gastric acid secretion (Fig. 3B).

To determine the associated functional genomic alterations of these 27 genes, the previous studies that reported them were analyzed. In the present study, there were no somatic mutations in the 27 genes, however, 12 of the genes, including gastrokine 2 (GKN2), prosaposin-like 1 (PSAPL1), ethanolamine-phosphate phospho-lyase (ETNPPL), cytochrome P450 family 2 subfamily W member 1 (CYP2W1), SFRP4, collagen type XXII α 1 chain, prostate stem cell antigen, lipase F, gastric type, solute carrier family 15 member 1 (SLC15A1), ATPase H+/K+ transporting subunit β (ATP4B), aldehyde dehydrogenase 3 family member A1 (ALDH3A1) and cystatin SN (CST1) were mutated in germline. Notably, two (SFRP4 and SLC15A1) of the 12 germline mutated genes have been reported in previous studies. Germline mutations in SFRP4 (P320T and R340K) were reported to be strongly associated with human cancer types (45) while mutations in SLC15A1 were linked with bowel diseases (46).

The KM plotter database was subsequently used with the log-rank test to investigate the association between the expression level of the 27 genes and clinical outcomes (both disease-free and overall survival times). In the database, patients with GC are divided into high and low expression groups, with each percentile of mRNA expression between the lower (25%) and upper quartiles (75%) used as a cut-off point. A total of 10 genes (ALDH3A1, ETNPPL, activity regulated cytoskeleton associated protein (ARC), ATP4B, CST1, CYP2W1, GKN2, PSAPL1, SFRP4 and SLC15A1) with differential expression were significantly associated (P<0.05) with survival times. In all the 10 genes, those with high expression (red curve) were correlated with a poorer survival while genes with low expression (black curve) correlated with higher survival in patients with GC (Fig. S2B).