The data on SCNAs in STAD patients, overall survival of STAD patients, gene expression in STAD patients, and locations of genes were obtained from TCGA (https://tcga-data-secure.nci.nih.gov/tcgafiles/tcga4yeo/tumor/prad/cgcc/broad.mit.edu/genome_wide_snp_6/snp/). The data was available in appendices. The version of the genome was hg19 (https://genome.ucsc.edu/cgi-bin/hgTracks?db=hg19&lastVirtModeType=default&lastVirtModeExtraState=&virtModeType=default&virtMode=0&nonVirtPosition=&position=chr2%3A25383722%2D25391559&hgsid=1653852260_7skZ6eFcyA94KTIXSNz9IEKYP60J). The data of SCNAs in STAD cell lines and expressions of relevant genes were obtained from CCLE (https://depmap.org/portal/download/all/). The two sets of data were analyzed separately.

The data was imported into RStudio 4.2.2 and classified into two sets: large-scale SCNAs with a width >100 kb and SCNAs with a width≤100 kb (McCarroll et al., 2008). They were analyzed separately to better understand their characteristics and implications. For ease of reference, SCNAs with a width >100 kb and SCNAs with a width≤100 kb were called Set1 and Set2, respectively (Table 1). See supporting material.

We hypothesized that the rate of DSB follows a Poisson Distribution with the parameter λ. To obtain λ based on the number of DSB events in each locus, a log-likelihood function was utilized to estimate λ values, taking into account the deviation from a perfect Poisson Distribution observed in the data (van Opheusden et al., 2020).

The study analyzed two datasets to determine the burden of SCNAs in specific genomic regions. The SCNA burden, defined as the number of SCNAs, was calculated for each patient within 1 Mb windows along the chromosomes.

To determine whether an SCNA fell within a window, the following criteria were applied: if the start point of the SCNA was smaller than the endpoint of the window, and the endpoint was larger than the start point of the window, then the SCNA was counted within the window.

Using RStudio, λ was estimated for each window. The parameter λ, representing the average occurrence of SCNAs within a window, was calculated using log-likelihood functions. It was observed that the distribution of λ within the windows did not perfectly match the Standard Poisson Distribution models.

To visually represent the distribution of λ, histograms of log2(λ) were plotted. Additionally, Manhattan plots were generated, using a threshold value approximately equal to the minimum value of λ observed in the hotspots, to illustrate the hotspots.

Hotspots were determined based on the top 5% of the λ values, as this threshold was considered statistically significant. The information regarding the hotspots in both sets of experiments can be found in the appendices, which include details such as the location of the hotspot and the corresponding λ value. We also identified the relevant genes located within these hotspots.

Before doing the GO analysis, hyper-polymorphic genes such as genes in OR, TAS, IGH, IGK, IGL, HLA family, and genes in the list of HLA (See appendix) were excluded from the analysis. Odor genes and taste receptor (TAS) (Keller et al., 2007; Carrai et al., 2017), HLA molecules (Dendrou et al., 2018), IGK rearrangements (Jackson et al., 2012), IgH class switching (Stavnezer and Schrader, 2014), Immunoglobulin light (IgL) chains (Edholm et al., 2011) were proven to be highly variable in all individuals. Their variability could potentially interfere with the GO analysis. By excluding these genes, we aimed to focus on identifying less variable but significant genes associated with the progression of stomach adenocarcinoma (STAD) and determining corresponding significant functional terms.

After modification, GO analysis was performed on GProfiler (Raudvere et al., 2019). The website was available at: https://biit.cs.ut.ee/gprofiler/gost. We tested the statistical enrichment of expression genes related to hotspots in the KEGG path. The ratio of intersection size to term size was calculated. Terms with the top 10 ratio value were selected to be graphed in a bubble graph. The adjusted p-value of all terms was smaller than 0.05, which illustrated that the terms were significantly enriched by expressed genes.

Violin plots visualized the relationships between the expression of significant genes identified in GO analysis and SCNA burden. The patients were divided into two groups, namely, “High” and “Low”, based on their median SCNA burdens. Both consisted of an equal number of patients. In addition, T-tests demonstrated the relationship between the APOBEC family copy number statuses and SCNA burden. If the start point of the SCNA was smaller than the endpoint of the analyzed APOBEC gene, and the endpoint of the SCNA exceeded the start of the APOBEC gene, patients with this SCNA were classified as “Mutant”, while the other patients were classified as ‘Wild’. T-tests with p-values less than 0.05 were considered statistically significant, and corresponding genes were deemed significant.

To explore the relationship between SCNA burden and overall survival, we utilized the R packages survival (version 3.4-0) and survminer (version 0.4.9) to generate Kaplan-Meier curves. Furthermore, patients were categorized based on the expression of APOBEC genes for survival analysis. Each survival analysis comprised three groups: Low, Median, and High. Results with p-values smaller than 0.05 were deemed statistically significant.

According to data from TCGA and CCLE, the majority of the log2(λ) values in Set1 (Width >100 kb) ranged from 0 to 0.5, while most of the log2(λ) values in Set2 (Width≤100 kb) were close to 0 (Figures 2A–D). Estimation revealed that the minimum value of λ in hotspots of TCGA Set1 was around 1.4, while the CCLE Set1 was around 1.6. In TCGA Set2 and CCLE Set2, the minimum value may be 1. Since the minimum value of Set2 was close to 1, which represents the overall minimum value of λ, the Manhattan plot of Set2 won’t be displayed. The Manhattan plot of Set1 with a threshold value of approximately 1.4 and 1.6, respectively, was plotted (Figures 2E, F).

A total of 145 hotspots were identified from the analysis of 2,897 windows in Set1 and Set2. In TCGA Set1, 1,437 genes were found within these hotspots, and after modification, 1,230 genes remained. Similarly, in TCGA Set2, 452 genes were located within the hotspots, and 395 genes remained after modification. In CCLE-Set1, 1,462 genes were found within these hotspots, and after filtering for genes whose variability could potentially interfere with the GO analysis, 1,265 genes remained. In CCLE Set2, 1,458 genes were found within these hotspots, and 1,336 genes were left after the filtering (Table 2). The complete list of genes before and after modification for both sets of experiments can be found in the appendices, along with additional information on the hotspots.

The GO analysis revealed that DNA Cytosine Deamination appeared among the top 10 terms in Set1. However, the top 10 terms in Set2 did not show any significant functions, leading to the exclusion of Set2 from further analysis. Interestingly, the GO analysis also highlighted the presence of genes belonging to the APOBEC family in both TCGA-Set1 and CCLE-Set1. These APOBEC genes were found to be associated with hotspots that exhibited a relatively low value of λ compared to other hotspots. The analysis revealed that these genes were predominantly located in hotspots with λ values around 1.4 and 1.6 in two sets of data, respectively (Figures 2G, H).

In TCGA data, the results indicated a significant association between the expression levels of APOBEC3C, APOBEC3D, APOBEC3F, APOBEC3G, APOBEC3H, and SCNA burden (Figures 3A–E). The findings also revealed a significant association between DSB points affecting APOBEC3D and SCNA burdens (Figure 3F). In CCLE data, the results showed a significant association between the expression levels of APOBEC3D, APOBEC3F, APOBEC3G, APOBEC3H, and SCNA burdens (Figures 3G–J).

The cohort of 332 patients from TCGA was divided into three groups based on their SCNA burden. The survival analysis demonstrated that patients with higher SCNA burden tended to better overall survival compared to those with lower SCNA burden (Figures 4A, B). Additionally, we investigated the correlation between the expression of specific APOBEC genes and the overall survival of STAD patients. Similarly, the cohort of 332 patients was divided into three groups based on SCNA burden. The results indicated that patients with median expression of APOBEC3C tended to better overall survival, as opposed to those with high or low expression (Figures 4C, D).