Cell lines and culture
The human GC cell line MGC-803 was purchased from the Cell Bank of the Chinese Academy of Science (Shanghai, China) and cultured in RPMI-1640 supplemented with 10% fetal bovine serum (FBS, Hyclone, Logan, UT, USA) at 37℃ with an atmosphere of 5% CO2. Cisplatin-resistant MGC-803 cells (MGC-803/DDP) were established by our laboratory. Briefly, MGC-803 cells were initially exposed to 0.03 μM of DDP for 2 weeks. The drug concentration was intermittently increased up to 3.33 μM over one year. To maintain the chemoresistant phenotype, MGC-803/DDP cells were cultured in complete culture medium with 1.67 μmol/L DDP. These two GC cells were stored in liquid nitrogen for further experiments. The complete experiment workflow is shown in Fig. 1.
The half-maximal inhibitory concentration (IC50) analysis by CCK8 assay
For cisplatin IC50 analysis, MGC-803 and MGC-803/DDP cells were seeded into 96-well plates and treated with the indicated concentration of cisplatin for 48 hours. 10 μl Cell Counting Kit-8 (CCK-8) solution (Dojindo, Tokyo, Japan) was added to each well. After 2 hours of incubation at 37 ℃, the optical density (OD) values were measured at a wavelength of 450 nm by Thermo Scientific Varioskan Flash (Thermo Fisher Scientific, USA). Cell growth inhibition rates were described as cell inhibiting curves and The IC50 parameters were calculated by GraphPad Prism V8 (GraphPad Prism, Inc., La Jolla, CA, USA).
RNA extraction, quality control, library construction and microarray analysis
Total RNA from MGC-803 and MGC-803/DDP cells was isolated using TRIzol reagent (Life Technologies, CA, USA) according to the manufacture’s instruction. The RNA concentration and quality were assessed with the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA). All RNA samples had an RNA Integrity Number (RIN) > 8.0. RNA-seq was assessed to identify the differentially expressed (DE) lncRNAs, mRNAs and miRNAs in MGC-803 and MGC-803/DDP cells. NEBNext® UltraTM II RNA Library Prep Kit for Illumina® (NEB, USA) was used with 2 μg of total RNA for the construction of sequencing libraries on the Illumina NovaSeq 6000-sequencing platform (Tsingke Biotechnology Co. Ltd, Beijing, China). The raw reads were filtered and cleaned by removing the adaptor reads and low-quality tags. Then clean paired-end reads were aligned to the reference genome using HISAT2 (v2.0.4) for transcripts and Bowite (v2.2.5) for miRNA as described [18, 19]. Subsequently, raw counts were generated by the StringTie (v.2.1.2) for transcripts and the FeatureCounts (v1.6.1) for miRNAs [20, 21].
Identification of DERNAs
To screen the DE lncRNAs, miRNAs and mRNAs in MGC-803/DDP in comparison with MGC-803 cells, the expression profiles were analyzed using the edgeR package [22] installed in R (version 4.0.5, www. r-project. org). Fold change (log2 absolute) ≥ 2 were considered as the thresholds for identifying DERNAs. Heatmaps were drawn using the “pheatmap” R packages to show DERNAs.
Functional annotation and pathway enrichment of the DEmRNAs
To investigate the biological functions implicated with the DEmRNAs, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis and Gene Ontology (GO) functional analysis were performed using the commonly used online database DAVID [23] (version 6.8; https://davidd. ncifcrf.gov/). The Gene Ontology Biological Process (GO-BP) term and pathway with the enriched gene count ≥ 2 and the significance threshold P < 0.05 were considered significant. The top 20 significant biological processes or pathways were visualized by Bubble maps created by “ggplot2” package in R software.
Construction of protein-protein interaction (PPI) network of DEmRNAs
The PPI network of DEmRNAs was constructed by the Search Tool for the Retrieval of Interacting Genes (STRING, version 11.0; https://string-db.org/cgi/input.pl). The input gene set was set as all DEmRNAs, and the species was set as homo sapiens. The minimum required interaction score was set as the highest confidence (0.900). The PPI network was visualized by Cytoscape software (version 3.8.2; https://cytoscape.org/). The CytoNCA plug-in in Cytoscape software was implemented to determine the network topology properties of the nodes, and the parameters were set as “without weight”. We used the MCODE plug-in in Cytoscape software to analyze functionally related and highly interconnected modules from the PPI network with the threshold score ≥ 12. Then hub genes from the significant clustering modules were subjected to KEGG pathway analysis with the enriched gene count ≥ 2 and the significance threshold P < 0.05 as the cut-off criterion.
Prediction of DEmiRNAs regulatory relationship
The miRNA-lncRNA regulatory relationship of DEmRNAs was predicted using the miRcode database (http://www.mircode.org/). The predicted miRNA-lncRNA regulatory relationship was integrated with DElcnRNAs and DEmiRNAs to obtain the DEmiRNA-DElncRNA relationship. In the light of this regulatory relationship, these screened DEmiRNAs were further used to predict the miRNA-mRNA regulatory relationship via four publicly profile datasets (TargetScan, miRDB, miRTarbase and ENCORI). The screening criteria were the number of supporting datasets ≥ 3. The predicted miRNA-mRNA regulatory relationship was integrated with DEmRNAs to obtain the DEmiRNA-DEmRNA relationship.
KEGG pathway analysis of DEmiRNAs and DElncRNAs
Based on the obtained DEmiRNA-DElncRNA and DEmiRNA-DEmRNA regulatory relationships and the co-expression relationship between DEmRNAs and DElncRNAs, the compareCluster function in the clusterProfiler package (http://bioconductor.org/ packages/release/bioc/html/clusterProfiler.html) in R software was performed to calculate and compare KEGG pathway enrichment analysis of screened DEmiRNAs and DElncRNAs. The pathway with the threshold P < 0.05 was considered significant. Bubble maps are implemented to indicate pathways in which DEmiRNAs and DElncRNAs may be implicated.
Construction of the ceRNA regulatory network
According to the ceRNA hypothesis, lncRNAs can function as decoys for various miRNAs and modulate target genes’ stability or translation [23]. We used the obtained DEmiRNA-DElncRNA and DEmiRNA-DEmRNA regulatory relationships and the correlation coexpression of DElncRNA-mRNA pairs to construct a lncRNA-miRNA-mRNA regulatory network. Because DEmiRNAs expression was negatively associated with DElncRNAs or DEmRNAs based on the ceRNA hypothesis, the positive correlation expression of DEmiRNA-DElncRNA pairs and DEmiRNA-DEmRNA pairs were excluded from our ceRNA network. We conducted the Cytoscape software (version 3.8.2) to visualize the ceRNA network, and the cytoHubba plug-in was implemented to rank nodes by their network features and select the top 50 genes from maximal clique centrality (MCC) as the hub genes [24]. To screen potential ceRNA regulatory axes, R package “ggalluvial” and “ggplot2” were performed to identify the ceRNA axes.
GDC GC data validation and identification of prognostic predictors
For data validation, the RNA-seq and miRNA expression profiles and clinical information of gastric cancer were downloaded from the Genomic Data Commons (GDC) data portal (https://gdc.xenahubs.net). The lncRNAs and mRNAs data contained 343 GC tissues and 30 matched non-cancerous tissues, and the miRNAs data included 410 GC samples and 42 adjacent cancer normal samples. The differentially expressed lncRNAs, mRNAs and miRNAs between GC and normal tissues were screened out using the edgeR package (http://bioconductor.org/packages/edgeR/) in R software [22]. The survival package was applied to explore the prognostic value of DElncRNAs, DEmiRNAs and DEmRNAs in the ceRNA network, with clinical data of GDC combined. This test is based on the Kaplan-Meier method and P value < 0.05 was considered statistically significant. Regression analysis was performed to detect the correlation between DElncRNAs and DEmRNAs expression from GDC GC data with P < 0.01 and r > 0.3 as the criteria for significance.
Statistical analysis
All statistical analyses were conducted using SPSS 24.0 (SPSS, Chicago, IL, USA), GraphPad Prism V8 (GraphPad Prism, Inc., La Jolla, CA, USA) and R software (version 4.0.5). Data are presented as the mean ± SEM. The two-tailed Student’s t-test and one-way ANOVA analysis were used to compare the different expression levels between two groups. Pearson correlation analysis was performed to evaluate the ceRNA network construction. The Kaplan-Meier method and log-rank test were employed to generate the survival curve and compare differences between survival curves, respectively. A P value < 0.05 was considered statistically significant.