Preprocessing of sequencing reads
AGS cells grouped into untreated and treated conditions with PRU for 24h (Fig. S2a) were sequenced using Ilumina novaseq6000 platform. Each condition were sequenced in triplicates. The sum of all the raw sequences were found to be 324,191,686 and the genome mapped sequences were 289,528,324, respectively. Samples from each conditions are mapped with average total reads 89.32% and average clean reads 95.98% as shown in Fig. S2b. The obtained cleaned reads were subjected to normalization based on log2 transcripts per kilo base million (TPM) (Fig. S2c). The samples were clustered based on Euclidean distance method and represented by dendogram as shown in Fig. S2d. The correlation of the sequenced samples is represented by heat map (Fig. S2e) based on read counts ≥5 and TPM ≥0.3, respectively. The overall read mapping statistics of all the samples are depicted in Table S1.
Screening of differentially expressed genes
The differentially expressed genes among PRU treated and untreated control AGS cells were identified using the edgeR package in R. The expression profile showed a total of 1,118 DEGs among the control and treated conditions which are clustered using complete linkage method as shown in Fig. 1a. Among the DEGs, the up-regulation and down-regulation category of genes were distinguished based on a logFC>=2 and false discovery rate (FDR) < 0.05 represented in a volcano plot (Fig. 1b). Upon identification, 463 genes were up-regulated and 655 were down-regulated as depicted in a bar graph in Fig. 1c.
Functional enrichment of DEGs
The DEGs were subjected to functional enrichment analysis with the up-regulated and down-regulated cluster of genes using ClusterProfiler package in R. The gene ontology functions predicted within the q-value 2.0 were taken into consideration. The enriched biological processes of 463 up-regulated genes are found to be involved in cell matrix adhesion, cellular modified amino acids metabolic process followed by positive regulation of chromosome organization enriched in cell-matrix adhesion shown in Fig. 2a. Similarly, the 655 down-regulated genes were enriched in the categories of negative regulation of hydrolase activity, small molecule catabolic process, protein targeting followed by regulation of response to DNA damage stimulus shown in Fig. 2b.
Pathway enrichment of DEGs
To obtain detailed information on the pathways executed by the DEGs, pathway enrichment was performed using KEGG database44. The top 10 significantly enriched pathways in DEGs with the number of genes along with FDR values are identified and listed in Table 1. Among the 10 pathways, maximum number of genes with higher FDR values were found to be involved in four distinct mechanisms: necroptosis, TNF signaling, MAPK signaling and Ubiquitin proteolysis represented in Fig. 3.
Expression of necroptosis related genes
With the significant pathway information involved in necroptosis signaling, the available necroptosis related genes were identified among the DEGs. Interestingly around 40 genes were found to be related with necroptosis as shown in Fig. 4. The enriched expression of genes involved in the necroptosis signaling network were mapped using KEGG pathway maps using Pathview R and visualized in Fig. 5.
Functional annotation of the necroptosis genes
The significant necropsis DEGs were subjected to functional analysis on GeneCodis. The top 10 GO terms of the necroptosis DEGs were listed and provided in Table S2–S4. Based on the constructed network clusters, the biological process of the genes were significantly enriched in identical protein binding (GO:0042802), tumor necrosis factor receptor binding (GO:0005164), protein binding (GO:0005515), protein-containing complex binding (GO:0044877), death receptor binding (GO:0005123), ubiquitin protein ligase binding (GO:0031625), death effector domain binding (GO:0035877), death domain binding (GO:0070513), JUN kinase kinase kinase activity (GO:0004706) and thioesterase binding (GO:0031996) shown in Fig. 6a.
Regarding MF, Regulation of tumor necrosis factor-mediated signaling pathway (GO:0010803), Positive regulation of I-kappaB kinase/NF-kappaB signaling (GO:0043123), Regulation of necroptosis process (GO:0060544), Necroptotic process (GO:0070266), Death-inducing signaling complex assembly (GO:0071550), TRIF-dependent toll-like receptor signaling pathway (GO:0035666), Regulation of extrinsic apoptotic signaling pathway via death domain receptors (GO:1902041), I-kappaB kinase/NF-kappaB signaling (GO:0007249), Apoptotic signaling pathway (GO:0097190) and Cellular response to mechanical stimulus (GO:0071260) shown in Fig. 6b.
Followed by, in terms of cellular component the genes were enriched in death-inducing signaling complex (GO:0031264), cytoplasm (GO:0005737), membrane raft (GO:0045121), ripoptosome (GO:0097342), cytosol (GO:0005829), CD40 receptor complex (GO:0035631), protein-containing complex (GO:0032991), tumor necrosis factor receptor superfamily complex (GO:0002947), cytoplasmic side of plasma membrane (GO:0009898) and plasma membrane signaling receptor complex (GO:0098802) shown in Fig. 6c.
Identification of target gene and its associated prognostic markers
To determine the hub target gene among the significant necroptosis genes, a protein-protein interaction network was constructed using the STRING database with the threshold confident score of 0.9 and all the unconnected nodes were removed. Subsequently, on the observation of the PPI network, RIPK was found to be the core hub gene interconnected with all the necroptosis proteins (Fig. 7a). Among the total DEGs, 16 candidate genes were identified to be in close association with the target RIPK family of genes. In which, 9 genes (RIPK1, EXD2, BBS7, LRRC75-AA51, GPR107, TUBA4A, KDM1B, TYMP, and MATR3) were found to be up-regulated and 8 genes (NRP1, MNX1AS1, SSRP1, PRDX2, PLRG1, LGALS4, SNX5 and FXYD3) were found to be down-regulated in PRU treated conditions presented in Fig. 7b, which can be considered as possible candidate biomarkers. The List of candidate genes significantly associated with RIPK family with their functional description is represented in Table S5.
Candidate gene validation through GEPIA
Consistent with the GEO analysis, GEPIA analysis showed among the 16 candidate genes, 8 genes (NRP1, MNX1-AS1, SSRP1, PRDX2, PLRG1, LGALS4, SNX5, and FXYD3) were overexpressed in stomach cancer samples compared with normal tissues. The box plots of gene expression by pathological stages based on the TCGA clinical annotation revealed their high expression levels significantly associated with advanced adenocarcinoma stage (P-value<0.05) are shown in Fig. S3. In addition, the overall survival rate for the eight candidate genes were predicted with their HR and Logrank P ratios as follows: HR=0.68, P=0.017 (NRP1), HR=1.8, P=0.00019 (MNX1), HR=0.84, P=0.28 (SSRP1), HR=0.88, P= 0.42 (PRDX2), HR=1.1, P=0.55 (PLRG1), HR=0.82, P= 0.21 (LGALS4), HR=1.0, P=0.9 (SNX5) and HR=1.2, P=0.28 (FXYD3). The survival plots generated using GEPIA are provided in Fig. S4.