The GSE62254 cDNA microarray dataset was downloaded from the Gene Expression Omnibus website (https://www.ncbi.nlm.nih.gov/geo/). Corresponding patient information, including Lauren classification and survival data, was obtained from the supplementary materials of the original article (25). The Robust Multichip Average algorithm (26) was applied for background correction, and qspline was applied for normalization (27). Data were perfect match-corrected and summarized using the Li-Wong method (28). All probes were mapped to Ensembl Gene Symbols in the R package ‘mygene’ (29).

Position weight matrix data of TF-binding motifs in vertebrates were obtained from the JASPAR database (30). The methods and parameters of binding site scanning used were previously described (31). PPI data were obtained from a publicly available dataset (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2836267/ bin/NIHMS177825-supplement-03.xls) (32).

Networks of DF, MX and IT GC subtypes were constructed by combining the corresponding gene expression, TF motif and PPI data with an update parameter of α=0.25. Confident TF-target edges were identified by a False Discovery Rate (FDR) of <0.05.

AnaPANDA software (23) was used to further identify TFs specifically activated in a certain subtype of GC, and the probability cutoff was set to 0.8 to build sub-networks. The hypergeometric distribution model was utilized to evaluate the overlap between genes co-targeted by each two given TFs.

Overall survival (OS) was the primary endpoint in the present analysis, which was defined as the time from tumor resection to death or last follow-up. The median mRNA expression level of a given gene was chosen as the cutoff to divide patients into two subgroups. Log-rank tests and Kaplan-Meier plots were used to evaluate the difference in OS between subgroups. Cox proportional hazard model was applied for multiple-variants analysis, in which ‘backward LR’ stepwise logistic regression was used for variable selection.

Human GC cell lines [MGC803 (the cell line used has been authenticated by STR profiling) and SGC-7901] were purchased from the Cancer Institute and Hospital, Chinese Academy of Medical Sciences (Beijing, China). All cells were maintained in Dulbecco’s modified Eagle’s medium (cat. no. 10-013-CVR; Corning, Inc., Corning, NY, USA) supplemented with 10% fetal bovine serum (cat. no. 10270-106; Gibco; Thermo Fisher Scientific, Inc.) at 37°C in an incubator containing 5% CO₂. Approximately 5×10⁵ cells/well were cultured in a 6-well plate. After 24 h, cells were transfected with small interfering RNAs (siRNAs) (sequences: si-NC, 5′-UUCUCCGAACGUGUCACGUTT-3′; si-NFYA, 5′-CAAACAAUACCACCGUAUUTT-3′) using RNA-mate (BioChain Institute, Inc., Newark, CA, USA) (5 µg siRNA + 10 µl RNA-mate) for 6 h at 37°C according to the manufacturer’s protocol. siRNAs were purchased from Shanghai GenePharma Co., Ltd. (Shanghai, China).

Nuclear transcription factor Y subunit α (NFYA) antibody (cat. no. 12981-1-ap, 1:1,000) was purchased from Wuhan Sanying Biotechnology (Wuhan, China). GAPDH antibody (cat. no. A01020, 1:1,000) was purchased from ARP American Research Products, Inc. (Waltham, MA, USA). Proteins were extracted from the cells using lysis buffer (50 mM Tris-HCl, pH 7.4; 10 mM EDTA; 0.5% NP-40; 1% Triton X-100) and quantified using the bicinchoninic acid method, after which they were separated (50 µg) by 10% SDS-PAGE and transferred to PVDF membranes. The membranes were blocked with 5% bovine serum albumin (cat. no. A1933; Sigma-Aldrich; Merck KGaA, Darmstadt, Germany) for 1 h at 37°C, and were then incubated with primary antibodies for 2 h at 37°C, washed with PBS-1% Tween (PBST) five times (5 min/wash), and then incubated with a secondary antibody (cat. no. 7074, 1:2,000; Cell Signaling Technology, Inc.) for 30 min at 37°C. Subsequently, the membranes were washed a further three times with PBST (5 min/wash) and proteins were detected using an enhanced chemiluminescence kit (cat. no. 34076; Thermo Fisher Scientific, Inc.).

For the cell proliferation analysis, a total of 12 h post-transfection, MGC803 and SGC-7901 cells were cultured in 96-well plates at a density of 2×10³ cells/well. Cell viability was measured using a Cell Counting kit-8 (CCK-8; Dojindo Molecular Technologies, Inc., Kumamoto, Japan), according to the manufacturer’s protocol. For colony formation assay, MGC803 and SGC-7901 cells were cultured in 6-well plates at a density of 1×10³ cells/well. After 14 days, colonies were fixed in methanol, stained with 0.25% crystal violet for 10 min at room temperature and counted. All assays were repeated three times independently.

Categorical, baseline data were compared with Pearson’s χ² test or Fisher’s exact test, and continuous baseline data were compared with one-way analysis of variance followed by Student Nerman Keuls test in Table I. For comparisons of continuous data, an unpaired Student’s t-test was conducted. Wald test was performed to evaluate the overall multivariate Cox model. For all statistical analyses, P<0.05 was considered to indicate a statistically significant difference, and a cutoff value of FDR<0.05 was used for multiple comparison corrections. All statistical analyses were two-sided and performed using R Software 3.3.1 (www.r-project.org). R packages ‘VennDiagram’ and ‘ggplot2’ were used for data visualization; Mygene was used for gene symbol mapping; MASS and survival were used for survival analysis; q value was used for FDR analysis.

A total of 300 patients were included in the present analysis. Patient characteristics are presented in Table I. In total, there were 134 patients in the DF group, 146 in the IT group and 20 in the MX group. The age in each group was 58.44±12.53, 64.41±9.61 and 67.35±7.90 years, respectively. A larger proportion of patients were male (66.33%). Patients in the DF, IT and MX groups were similar with regards to tumor location, mutL homolog 1 expression and recurrence.

Expression data were extracted from 134 DF, 20 MX and 146 IT GC samples from the GSE62254 dataset. Combining TF motif and PPI data, TF-target regulatory networks for these three subtypes of GC were generated using PANDA software (Fig. 1). For each TF-target edge in each subtype, a Z-score was calculated based on the confidence level of the potential regulatory relationship. Edge Z-score distributions of the various subgroups of GC are presented in Fig. 2A. Different subgroups were assigned different colors.

All edges with an FDR-adjusted P-value of <0.05 were considered confident. The overlap of confident edges among the three subtypes of GC was displayed as a Venn diagram (Fig. 2B); >85% of TF-target edges were shared among all three subtypes, indicating that the TF-target relationship was strongly conserved. According to the definition of DF-specific edges (ModuleDF), IT-specific edges (ModuleIT) and commonly conserved edges, different edges were marked with different colors in a 3D scatter plot, in which each axis represented a Z-score of each subtype of GC (Fig. 2C). Fig. 2D exhibited the projection of this 3D plot through each axis. The overlap of differentially expressed genes between DF and IT GC, as well as the target genes of ModuleDF and ModuleIT were illustrated in a Venn diagram (Fig. 2E). By applying the hypergeometric distribution model to the target genes of each TF, it was revealed that most TFs with a high activity in ModuleDF also had a high activity in ModuleIT (Fig. 2F).

Using the MX subtype as a control, the AnaPANDA algorithm was applied to further identify TFs, which were specifically activated in each subtype of GC. A total of 13 TFs were activated in DF GC (Fig. 3A), and 8 TFs were activated in IT GC (Fig. 3B). Additionally, Fisher’s exact test was applied to evaluate the overlap between target genes shared by different pairs of TFs. In DF GC, RELA proto-oncogene (RELA) and forkhead box L1 (FOXL1), sex determining region Y (SRY) and NK3 homeobox (NKX3)-2, NFYA, paired box 2 and cAMP responsive element binding protein 1 (CREB1) were identified as co-targeted, which suggested that those TFs had very similar target profiles (Fig. 3C). In IT GC, NK2 homeobox 5 and NKX3-1, transcription factor AP-2α (TFAP2A), early growth response 1 (EGR1) and Sp1 transcription factor (SP1) were also identified as co-targeted (Fig. 3D).

Among the TFs specifically activated in DF or IT subtypes, TFAP2A, FOXL1 and SP1 were identified as potential prognostic biomarkers in all GC (Fig. 4A-C). NFYA and SRY were identified to be potential prognostic factors in DF GC (Fig. 5A-D), whereas T-cell leukemia homeobox 1 (TLX1) identified to have potential prognostic value in IT GC (Fig. 5E and F). Meanwhile, heart and neural crest derivatives expressed 1 (HAND1) and CREB1 were also identified to be potential prognostic factors in DF GC (Fig. 6A-D), whereas EGR1 was identified to have potential prognostic value in IT GC (Fig. 6E and F).

Cox proportional hazards model was also applied and respectively implemented for the aforementioned genes. NFYA [hazard ratio (HR) (95% confidence interval, CI)=0.560 (0.349, 0.900), P=0.017] and SRY [HR (95% CI)=0.603 (0.375, 0.969), P=0.037] were identified as independent prognostic factors in DF GC (Tables II and III), whereas TLX1 [HR (95% CI)=0.547 (0.321, 0.9325), P=0.027] was identified as an independent prognostic factor in IT GC (Table IV). Conversely, EGR1 was not associated with prognosis in either DF or IT GC (Table V).

To confirm the biological function of NFYA in DF and IT GC, NFYA expression was knocked down by siRNA in DF GC-derived MGC803 cells and IT GC-derived SGC-7901 cells (Fig. 7A). CCK-8 assays indicated that knockdown of NFYA expression markedly inhibited the rate of cell growth in DF GC-derived MGC803 cells, and partially inhibited the rate of cell growth in IT GC-derived SGC-7901 cells (Fig. 7B). Colony formation assays also demonstrated that the colony formation abilities of DF GC-derived MGC803 cells were nearly eliminated by NFYA siRNA (Fig. 7C). However, the colony formation ability of IT GC-derived SGC-7901 cells was only partially inhibited under the same conditions (Fig. 7C).