A total of two GC cohorts with whole-exome sequencing (WES) data were used to construct the FAM signature, the one from The Cancer Genome Atlas (TCGA) as the discovery set, and the other from Affiliated Hospital of Jiangsu University (AHJU) as the validation set. The enrollment criteria for patients included the following: 1) matched transcriptome and WES data for FAM signature construction in the TCGA discovery set; 2) available WES data for FAMscore calculation; 3) pathological diagnosis of GC; and 4) no prior history of anti-cancer therapies (including neoadjuvant therapy). Clinical and clinicopathological classifications and staging were based on the American Joint Committee on Cancer criteria. The research protocol of the AHJU cohort was approved by the Ethics Committee of AHJU, and all patients provided written informed consent.

Multidimensional data of TCGA, including mRNA expression profiles, somatic mutations, and microsatellite instability (MSI), were downloaded from the cBioPortal database (http://www.cbioportal.org/). Tumor neoantigen burden (TNB) was obtained from the immune landscape of TCGA pan-cancer (Thorsson et al., 2018). Transcriptome data of AHJU were stored in the dataset-EGAD00001004164 in the European Genome-phenome Archive. Genome data of AHJU were stored in the dataset-HRA001647 in the Genome Sequence Archive for Human.

A total of four published immunotherapy cohorts of solid tumors with available data on WES and response evaluation and/or survival were downloaded and analyzed to determine the role of FAMscore, including 1 GC cohort wherein the participants had received PD-1 inhibition with pembrolizumab (Kim et al., 2018), one melanoma cohort wherein participants were treated with pembrolizumab or nivolumab (Hugo et al., 2017), one non-small-cell lung cancer (NSCLC) cohort (NSCLC-1) wherein participants were treated with combined PD-1 (nivolumab) and CTLA-4 (ipilimumab) blockade (Hellmann et al., 2018), and another NSCLC cohort (NSCLC-2) wherein participants were treated with pembrolizumab alone (Rizvi et al., 2015).

A total of 34 patients with GC underwent RNA-seq in the AHJU cohort. The RNA-seq process has been reported previously (Duan et al., 2020). In brief, the library was prepared using the KAPA Stranded RNA-seq Kit with RiboErase (KAPA Biosystems, United States), after extraction of total RNA from fresh samples and depletion of ribosomal RNA. Then, the library was sequenced on Illumina HiSeq4000 NGS platforms (Illumina, United States), after determination of the library concentration and assessment of library quality. The sequence reads were generated by base calling performed on bcl2fastq v2.16.0.10 (Illumina, United States) in the FASTQ format (Illumina 1.8 + encoding). Following quality control and isoform and gene-level quantification, the transcriptome was mapped using STAR (version 2.5.3a).

A total of 74 patients with GC underwent WES in the AHJU cohort. As previously described (Wang et al., 2021), after genomic DNA extraction, the whole genome library was prepared using the KAPA Hyper Prep Kit (KAPA Biosystems), and exome capture was performed using the Illumina Rapid Capture Extended Exome Kit (Illumina Inc.). Then, enriched libraries were sequenced on the Illumina HiSeq4000 NGS platform (Illumina, United States) as paired 150-bp reads. The data processing has been described elsewhere (Xu et al., 2018).

We selected the xCell algorithm to quantify the abundance of infiltrating cells in the TME based on transcriptome data. It uses a novel technique for reducing associations between closely related cell types and validates signatures using both in silico simulations and cytometry immunophenotyping (Aran et al., 2017).

We included somatic nonsynonymous variants (SNVs) in coding regions for analysis, including frameshift, missense, nonsense, and splice-site mutations, which may be functional (Devarakonda et al., 2018). The tumor mutation burden (TMB) was calculated as the number of SNVs.

To facilitate detection, genes were asked to have a mutation frequency of ≥5%. Mutations (no matter the number) that occurred in one gene of the same sample were evaluated as 1 or 0. Genetic mutations associated with overall survival (OS) were selected using univariate Cox regression models. Then, CAF abundances between the genetically mutated type and wild-type were compared, and genes with mutations that significantly impacted CAF abundance were selected. Furthermore, genes with the most useful prognostic mutations were screened using the least absolute shrinkage and selection operator (LASSO) Cox regression model, and a FAMscore model was subsequently built based on the fraction of selected genes using Cox regression coefficients. The formula is established as follows: FAMscore = sum (each gene’s mutation × corresponding coefficient).

After genomic DNA extraction from both GC and normal tissues, the MSI status was determined by single fluorescent multiplex PCR based on five well-known mononucleotide repeats, BAT-25, BAT-26, NR-21, NR-24, and NR-27 (Cristescu et al., 2015). MSI was defined as allelic size variations in at least two microsatellites in samples or microsatellite stability (MSS).

Based on the criteria of adjusted p-value<0.05 and log2 (fold change) > 1, differentially expressed genes (DEGs) were determined between the high and low FAMscore subgroups. NetworkAnalyst 3.0 (https://www.networkanalyst.ca/) was used for GSEA based on the Kyoto Encyclopedia of Genes and Genomes (KEGG).

Antibodies such as anti-CD8 (CST70306, Cell Signaling Technology, United States), anti-CD56 (CST3576), anti-panCK (CST4545), anti-S100 (ab52642, Abcam, UK), anti-HLA-DR (ab92511), and anti-CD68 (BX50031, Biolynx, China) were used in mIF using the PANO 7-plex IHC kit (Panovue, Beijing, China), according to the manufacturer’s instructions. The Mantra System (PerkinElmer, Waltham, MA, United States) was used to scan the stained slides and reconstruct the images of the sections. inForm image analysis software (PerkinElmer, Waltham, MA, United States) was used to quantify the cells in the images.

The χ² test, Student’s t-test, Fisher’s exact probability test, and Mann–Whitney U test were used for comparisons between groups. The Pearson correlation coefficient R was used to evaluate the correlations between groups. The Kaplan–Meier method and log-rank test were used for survival analysis. The optimal cut-off values to define high and low FAMscore in TCGA and AHJU cohorts were determined based on the association of FAMscore with OS using the survminer R package. Univariate and multivariate Cox proportional hazard models were used to evaluate prognostic factors, and hazard ratios (HRs) along with their 95% confidence intervals (CIs) were determined. The receiver operating characteristic curve (ROC) and the areas under the ROC curves (AUCs) were used to evaluate the predictive power of FAMscore for the objective response of patients to immunotherapy. Statistical significance was set at p < 0.05. R (version 4.0.1), R Bioconductor packages, and SPSS (version 19.0, Chicago, IL) were used in the analyses.

In the study, we included 383 patients from TCGA and 74 patients from AHJU. An overview of the patient characteristics is provided in Table 1, and the details are in Supplementary Table S1. No significant differences were found in patient characteristics between the TCGA and AHJU cohorts.

In the TCGA GC cohort, 828 genes presented a mutation frequency of ≥5% (Supplementary Table S2). Of these, mutations in 47 genes were significantly associated with OS (Figure 1A), and the CAF abundances calculated by the xCell algorithm (Supplementary Table S3) were significantly different between the mutated type and wild-type in 31 of these genes (Figure 1B). After LASSO Cox regression analysis (Figure 1C), 16 genes were selected to construct the FAMscore of OS, as follows: FAMscore = ANK1 mutation*-0.211 + AR mutation*-0.612 + CARD11 mutation*-0.077 + CDH23 mutation*-0.313 + CDHR2 mutation*-0.421 + COL5A2 mutation*-0.088 + FAT2 mutation*-0.033 + GRIN3A mutation*-0.134 + ITPR3 mutation*-0.287 + PARD3B mutation*-0.154 + PCDH20 mutation*-0.402 + PTPRT mutation*-0.019 + THSD1 mutation*-0.803 + UTRN mutation*-0.288 + YLPM1 mutation*-0.104 + ZNF462 mutation*-0.151.

The FAMscore in TCGA and AHJU is provided in Supplementary Table S1. In TCGA, FAMscore was significantly higher in histological grade III than in grade I/II (p = 0.047; Supplementary Figure S1A). A similar result was observed in AHJU, but the significance was limited by the sample size (p = 0.22; Supplementary Figure S1B). The FAMscore was not significantly associated with the other patient characteristics including age, sex, and stage.

In TCGA, high FAMscore was significantly associated with poor OS (p < 0.0001; Figure 1D), which was further validated in AHJU (p = 0.038; Figure 1E). In multivariate analysis, high FAMscore was an independent predictor of poor OS in both TCGA (HR = 4.36, 95% CI: 2.61–7.30, p < 0.001; Figure 1F) and AHJU (HR = 3.75, 95% CI: 1.01–13.95, p = 0.048; Figure 1G).

MSI GC had a significantly lower FAMscore than MSS GC in both TCGA (p = 9.11e-24) and AHJU (p = 0.023) (Figure 2A). TMB (Figure 2B) and TNB (Figure 2C) were negatively correlated with FAMscore in both TCGA (R = -0.69, p < 2.2e-16 and R = -0.6, p = 3.1e-10, respectively) and AHJU (R = -0.58, p = 6.9e-08 and R = -0.54, p = 7.4e-07, respectively). Moreover, the immune subtype C2 (IFN-gama dominant) in TCGA, which may override an evolving type I immune response (Thorsson et al., 2018), had the lowest FAMscore among immune subtypes (Figure 2D).

A total of 180 DEGs were identified between the high and low FAMscore subgroups in the AHJU cohort (Figure 3A and Supplementary Table S4), and the top 50 DEGs are shown in Figure 3B. GSEA showed that genes in KEGG pathways involving DNA repair, cell death, cell cycle, and aerobic metabolism were significantly enriched in the low FAMscore subgroup (Figure 3C), while those involving autoimmune diseases, infection, and signaling associated with malignant biology were significantly enriched in the high FAMscore subgroup (Figure 3D). These results suggest that GC with a high FAMscore is more aggressive and characterized by abnormal immunity.

mIF for selected immune cells was performed in samples from 47 patients in the AHJU cohort (Figure 4A and Supplementary Table S5). We found that the density of M1 macrophages in the tumor parenchyma was significantly higher in the low FAMscore subgroup than that in the high FAMscore subgroup (p = 0.008; Figure 4B). The low FAMscore subgroup also had a nonsignificantly increased density of NK cells (including CD56bright and CD56dim subtypes) in the tumor parenchyma. In TCGA (Figure 4C and Supplementary Table S6), cell abundances calculated by the xCell algorithm further verified that low FAMscore was significantly associated with increased infiltration of M1 macrophages (p = 0.003) and NK cells (p = 1.12e-06). These results indicate that a low FAMscore may reflect a favorable TME characterized by dynamic immune infiltration.

In four immunotherapy cohorts (Supplementary Table S7), patients were divided into partial response/complete response (PR/CR) and stable disease/progressive disease (SD/PD) subgroups. The PR/CR subgroups generally had lower FAMscores than the SD/PD subgroups (p = 0.014 for the GC cohort, p = 7.9e-05 for the NSCLC-1 cohort, p = 0.023 for the NSCLC-2 cohort, and p = 0.14 for the melanoma cohort; Figure 5A). The ROC analysis showed that the AUC for the prediction of FAMscore to immunotherapy responses ranged from 0.639 to 0.742 (Figure 5B). Based on the optimal threshold of FAMscore for the maximum ROC curve values (−0.315, −0.002, −0.3, and −0.225 for the GC, NSCLC-1, NSCLC-2, and melanoma cohorts, respectively), patients were dichotomized into high and low FAMscore subsets. Furthermore, the objective response rates (ORRs) between the low and high FAMscore subgroups were 66.7 vs. 14.6% (8/12 vs. 6/41, p < 0.001) for the GC cohort, 68.2 vs. 19.6% (15/22 vs. 9/46, p < 0.001) for the NSCLC-1 cohort, 75 vs. 23.1% (6/8 vs. 6/26, p = 0.007) for the NSCLC-2 cohort, and 75 vs. 40.9% (12/16 vs. 9/22, p = 0.037) for the melanoma cohort (Figure 5C). These findings clearly indicate that the FAMscore was correlated with immunotherapy responses.

Progression-free survival (PFS) was available in the NSCLC-1 and NSCLC-2 cohorts, and OS was available for the melanoma cohort. A high FAMscore was associated with poor PFS in both the NSCLC-1 (p = 0.019) and NSCLC-2 (p = 0.01) cohorts (Figure 6A,B). In multivariate analysis including PD-L1 expression, smoking history, and other clinical features, high FAMscore was the only independent predictor of poor PFS in both the NSCLC-1 (HR = 2.55, 95% CI: 1.16–5.62, p = 0.02; Figure 6C) and NSCLC-2 (HR = 5.0, 95% CI: 1.13–22.19, p = 0.034; Figure 6D) cohorts. Moreover, a high FAMscore was also associated with poor OS of immunotherapy in the melanoma cohort (p = 0.011; Figure 6E). In multivariate analysis including BRAF mutation, therapy history of MAPK inhibitors, and other clinical features, a high FAMscore was still the only independent predictor of poor OS (HR = 3.48, 95% CI: 1.27–9.55, p = 0.015; Figure 6F). These results revealed that the FAMscore was a strong predictor of immunotherapy-associated survival.