The gene expression profiles of GSE26899, GSE77346, GSE54129, and GSE65801 were obtained from the Gene Expression Omnibus (GEO; www.ncbi.nlm.nih.gov/geo). In detail, GSE26899 dataset is consisted of 96 clinical gastric tumor tissues and 12 adjacent normal tissues; GSE77346 dataset is consisted of 1 trastuzumab-sensitive cell line and 4 trastuzumab-resistant cell lines (12); GSE54129 includes 111 human gastric cancer tissues and 21 non-cancerous tissues; GSE65801 contains 32 gastric cancer tissues and 32 paired non-cancerous tissues (13).

The raw microarray data files of the datasets downloaded from the GEO website were subsequently analyzed via using the GEO2R (www.ncbi.nlm.nih.gov/geo/geo2r/), an online tool comparing two or more groups of samples in the same experimental setting (14). False Discovery Rate (FDR) of P-value adjusted (adj. P) to 0.05 and |logFC|>1 were set as the cut-off criteria.

Gene ontology (GO) analysis is a commonly used approach for functional studies with three ontologies including biological process, molecular function, and cellular component (15), while Kyoto Encyclopedia of Genes and Genomes (KEGG) is a knowledge base for the systematic study of gene functions (16). To study the functional annotations of differentially expressed genes (DEGs), we next employed Database for Annotation, Visualization and Integrated Discovery (DAVID, david.abcc.ncifcrf.gov/,) to process the GO and KEGG analyses of DEGs identified in gastric cancer samples. P<0.05 was set as the threshold.

The Search Tool for the Retrieval of Interacting Genes (STRING), an online database (string-db.org) designed to evaluate PPI information, covers 9,643,763 proteins from more than 2,000 organisms, which was used to construct the PPI. To evaluate the interactive associations of DEGs identified from GSE26899, we mapped these DEGs to the STRING (version 10.5) database. Confidence score >0.4 was selected as significant. PPI networks were constructed by STRING and visualized by Cytoscape. Subsequently, the plug-in Molecular Complex Detection (MCODE) was employed to screen the modules of PPI networks in Cytoscape with the threshold set as follows: MCODE scores >10.

To evaluate the association between COL4A1 level and its clinical outcomes, Kaplan-Meier plotter (KM plotter; www.kmplot.com), an online survival analysis tool, was performed. KM plotter is capable of assessing the effect of 54,675 genes on overall survival via using 10,188 cancer samples including 4,142 breast, 1,648 ovarian, 2,437 lung, and 1,065 gastric cancer patients (17). Patients with gastric cancer were separated into high- and low-expression groups according to the level of COL4A1, and the overall survival was then analyzed. The hazard ratio (HR) with 95% confidence intervals and log rank P-value were calculated.

GeneMANIA, an online tool (www.genemania.org/), can be used to generate hypotheses of gene function, analyze gene lists, and prioritize genes for functional assays (18). After selecting Homo sapiens from the nine optional organisms, COL4A1 was entered into the search bar and the results were then collected. Annotation of biological processes involving COL4A1 was performed by consulting the Coremine Medical online database (www.coremine.com/medical/).

To predict the miRNAs targeting the mRNA of COL4A1, miRWalk (version 2.0, zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2/), an online platform supplying information about predicted and experimentally validated miRNA-target interactions, was then employed (19). Herein, nine prediction programs (miRWalk, miRanda, miRDB, miRNAMap, Pictar2, PITA, RNA22, RNAhybrid and Targetscan) were selected. These predicted miRNAs were then overlapped by at least seven programs, and selected for further analysis. Pathway enrichment analysis of these miRNAs was performed by using the DIANA-mirPath web server (snf-515788.vm.okeanos.grnet.gr/index.php?r=mirpath) (20).

SPSS 22.0 software (IBM Corp., Armonk, NY, USA) was used to analyze data. Two tailed Student's t-test was used to compare the two groups. P<0.05 was considered to indicate a statistically significant difference.

Drug resistance-related genes identified from in vitro drug-induced resistant models may simply represent transcriptional changes, thus, complementary agents targeting these genes are usually failure to translate into clinical practice (21). Because drug resistance acquisition can arise before the malignant transformation stage (22), the genes playing important roles in tumorigenesis and drug resistance is thereby more likely to be critical for resistance occurrence. Therefore, we first identified 509 DEGs in gastric cancer tissues from the GSE26899 dataset using a 2-fold-change and adj. P<0.05 as the threshold cutoff. Among these DEGs, the expression of these 172 genes was significantly upregulated, while that of 337 genes was significantly downregulated in cancer tissues (Fig. 1A). To explore the potential roles of these significantly upregulated DEGs, GO and KEGG pathway analyses, which can provide valuable insights regarding protein function, were performed (23). As shown in Fig. 1B, GO results showed that these significantly upregulated DEGs were largely associated with extracellular region part, collagen, extracellular matrix, and cell or biological adhesion processes. KEGG results also consistently revealed that these DEGs were highly enriched in extracellular matrix-receptor interaction and focal adhesion pathways (Fig. 1C). Based on further DEGs analysis with the STRING database, the PPI network of these DEGs, which contained 505 nodes and 1,207 edges, was subsequently constructed. Using the MCODE plug-in in Cytoscape, we obtained the module with the highest score (Fig. 2A), and also performed a cluster analysis of these genes in the module (Fig. 2B). The module with the highest score was selected based on these DEGs identified from gastric cancer tissues via bioinformatics methods.

To further explore the potential genes contributing to trastuzumab resistance, we retrieved GSE77346 data, a microarray dataset consisted of one trastuzumab-sensitive gastric cancer cell line and four trastuzumab-resistant gastric cancer cell lines (12). After screening the significantly upregulated DEGs in trastuzumab-resistant cancer cells, COL4A1, one of hub genes in the selected module (Fig. 2A), was also found to be significantly upregulated in trastuzumab-resistant cells (Fig. 3A), suggesting COL4A1 might be important in both tumorigenesis and trastuzumab resistance. Since the gene expression is not always consistent with its protein amount (24), further validation of COL4A1 protein level in clinical gastric cancer tissues is quite necessary. By employing the Human Protein Atlas database, an online tool analyzing protein level from clinical specimens, we observed that COL4A1 was positively expressed in gastric cancer tissue, but negatively expressed in normal gastric tissue (Fig. 3B). Two other datasets, including GSE54129 and GSE65801 were also used to validate the mRNA level of COL4A1, which was indeed significantly upregulated in clinical gastric cancer samples (Fig. 3C and D). To investigate potential regulation mechanisms of COL4A1 in gastric cancer, the genomic alteration of COL4A1 in The Cancer Genome Atlas (TCGA) cohort was analyzed using the cBioPortal (www.cbioportal.org/) (25). As shown in Fig. 3E, the amplification, missense mutation, truncating mutation, and deletion of COL4A1 accounted for 4.3, 4.8, 2.0, and 0.5% of stomach adenocarcinoma cases, respectively. The prognostic value of COL4A1 in patients with gastric cancer was analyzed using the KM plotter according to the low and high expression of COL4A1. As shown in Fig. 3F, the high mRNA level of COL4A1 [HR 1.56 (1.31–1.86)] was associated with poor overall survival of gastric cancer patients. Altogether, these data suggest that COL4A1 is a potential gene candidate that promotes the development of gastric cancer and subsequent trastuzumab resistance.

As a user-friendly web interface for functional prediction of genes, GeneMANIA has been widely used as an effective tool to predict and explore drug resistance-related genes (18). Accordingly, COL4A1 showed interactions with 20 proteins/genes; among these 6 genes were involved in conferring drug resistance in cancer (Fig. 4A). The mRNA expression of these 6 genes in gastric cancer cells and trastuzumab-resistant gastric cancer cells were also evaluated using the microarray data of GSE77346. As shown in Fig. 4B, the mRNA level of angiopoietin 2 (ANGPT2), COL3A1, COL15A1, and secreted protein acidic and cysteine rich (SPARC) were upregulated in trastuzumab-resistant gastric cancer cells. Among these genes, COL4A1 co-expressed and co-localized with ANGPT2, collagen type VI α 3 chain (COL6A3), and SPARC, respectively. It has been reported bevacizumab resistance in glioblastoma and doxorubin resistance in liver cancer are largely mediated by an enforcement of ANGPT2/TIE2 signaling (26,27). Additionally, overexpression of COL6A3, one of the most highly upregulated genes in oxaliplatin and cisplatin resistant ovarian cancer cells, has also been recognized to confer cisplatin resistance in sensitive cancer cells (28). Moreover, the accumulation of intracellular SPARC can drive imatinib resistance in chronic myelogenous leukemia cells (29), and also modulate cisplatin resistance via modulating the Let-7f-1 miRNA/HMGB1 signaling in medulloblastoma cells (30). COL3A1 and COL15A1, two types of fibrillar collagen highly expressed in chemotherapeutic drugs resistant ovarian cancer cells (31), were also exhibited to co-express, share protein domains, and co-localize with COL4A1 (Fig. 4A). Additionally, COL4A1 also was co-expressed, shared pathways, and had physical interactions with platelet derived growth factor subunit B (PDGFB). Akt/PDGF-B signaling can regulate Akt, and thus confer the hypoxia-induced cisplatin resistance in liver cancer cells (32). PDGFB may contribute to the resistant phenotype and sustain signaling through MAPK and Akt in breast cancer cells (33). The Coremine Medical database is a freely available online tool for obtaining information on health, medicine, and biology (33); thus, we used this database to annotate the biological process of COL4A1. As shown in Fig. 4C, 12 biological processes significantly associated with COL4A1, gastric cancer, and drug resistance (P<0.01) were annotated. Considering that the close relationships of COL4A1 with these processes and the close relationships of the 12 processes with gastric cancer and drug resistance, COL4A1 might be involved in the development of drug resistance in gastric cancer via its effect on these biological processes. In detail, cell growth related biological processes (including 4 cell growth, cell proliferation, and apoptotic process), gene expression regulation-related (including gene expression, RNA interference, and reverse transcription), and especially, the epithelial-mesenchymal transition (EMT) biological process were identified to be closely associated with the development of COL4A1 in the gastric cancer drug resistance. Therefore, these results suggest that COL4A1 may confer drug resistance in gastric cancer via regulating cell growth, gene expression, and in particular, epithelial mesenchymal transition (EMT) processes.

Post-transcriptional gene expression regulated by miRNAs is important for multiple cellular processes during development and pathogenesis (34). Amplification and overexpression of tumor-promoting miRNAs or genetic loss of tumor-suppressing miRNAs are tightly correlated with the development of cancer and chemotherapeutic resistance (35). Therefore, studying target genes of miRNAs is a focus of interest due to the diagnostic and therapeutic relevance, and the gene functions can also be predicted according to functionality of these miRNAs targeting the gene (36). To obtain the miRNAs targeting COL4A1 mRNA, the miRNA-mRNA interaction analysis was performed by using the miRWalk. We identified 553 miRNAs predicted to transcriptionally target COL4A1 mRNA, indicating the regulation of COL4A1 by miRNAs. These miRNAs predicted by at least 8 of 9 prediction tools were then submitted to pathway enrichment using DIANA miRPath (20). As shown in Table I, the top 5 highly associated pathway were ECM-receptor interaction, TGF-β signaling pathway, viral carcinogenesis, proteoglycans in cancer, and Hippo signaling pathway, which are almost reported to be related with drug resistance (37–41). It has been established TGFβ signaling can activate autophagy process, and thereby lead to the oxaliplatin resistance in colorectal cancer (38). In addition, the activation of Hippo signaling also contributes to the drug resistance and cancer relapse (42). YAP, an effector of hippo signaling, has also been reported to alter clinical response of EGFR-tyrosine kinase inhibitors in lung cancer patients (43). Furthermore, it has also been reported that inactivation of Hippo pathway can restore gemcitabine sensitivity among a variety of cancers (44). In the present study, as shown in Table II, 9 of top 10 miRNAs that target COL4A1 mRNA were tightly associated with drug resistance in cancers (45–51). For instance, loss of intracellular miR-29b promotes cisplatin resistance in gastric cancer. Consistently, ectopic overexpression of miR-29b in cholangiocarcinoma cells can confer gemcitabine sensitivity to HuH28 cells (45). MiR-506 overexpression can confer hydroxycamptothecin resistance in colon cancer cells by inhibiting PPARα expression (50). Collectively, these data provide further supports for the drug resistance function of COL4A1 in gastric cancer.