Uncovering the link between human endogenous retroviruses, inflammatory pathways, and gastric cancer development

2.1Compilation and analysis of transcriptome datasets

Transcriptome sequencing datasets were compiled from various sources, starting with 56 paired samples from gastric cancer patients, obtained from the NCBI database under the accession number PRJNA435914. To further validate the biomarkers identified, additional data from a separate gastric cancer research project (PRJNA1013242), comprising 61 normal and 48 tumor samples (totaling 109 cases), was also included. For lung cancer, datasets from 20 paired samples were retrieved under the project PRJNA344450 [20], providing a comparative analysis across different cancer types. The liver cancer dataset (PRJNA598835 [21]) included 35 tumor and 30 normal samples, offering insights into the expression patterns in hepatic malignancies. Breast cancer research data (PRJEB50623) contributed to 31 tumor and 10 normal samples, enriching the diversity of the dataset. Pancreatic cancer samples (PRJNA490335 [22]), with an equal distribution of 10 tumor and 10 normal samples, and colorectal cancer data (PRJNA288518 [23]), with 5 tumor and 5 normal samples, were also compiled. These datasets served as an external control to validate the specificity and reliability of the biomarker expression levels identified in our study.

2.2HERVs annotation and differential expression analysis

The various types of HERVs were identified on the human-reference genome (GRch38) using RepeatMasker software (V4.0.9_p2) and Dfam gene component database [24]. To ensure that the HERV position is unique, the combination of the HERV name and position was adopted as keywords in querying the annotation to display.

Figure 1.

The differential expression analysis of HEARVs. (a) Differential HERVs in the gastric cancer dataset were screened based on the criteria of |log2FC|> 1 and p< 0.05. (b) The intersection of the results obtained from the SVM and AdaBoost algorithms. (c) The expression profiles of three potential HERV markers in normal and tumor groups of the gastric cancer dataset. The expression data are presented as log-transformed TPMs. (d) The expression profiles of the three potential HERV markers in normal and tumor groups of lung adenocarcinoma samples. The expression data are presented as log-transformed TPMs.

In the transcriptome analysis, the HERVs sequences were retrieved and decompressed to fastq, and the Hisat2 (V2.2.1) [25] was adopted to query the HERVs data against the reference genome(GRch38). The Featurecounts (V2.0.1) [26] was used to quantify the aligned reads of each protein-coding and non-coding genes using EM algorithm. Based on the count matrix, the differential expression analysis was performed by DESeq2 package (V1.32.0) [27], and a combination of |log2FC|> 1 and P< 0.05 was used as cutoff to select the significantly differentially expressed genes.

2.3Identification of HERV biomarkers

In the feature selection process, a combination of SVM [28] and AdaBoost [29] algorithm was employed to develop an efficient, interpretable, and robust statistical model. This method created perturbed samples by randomly altering the sample labels and multivariate screening of variance fractions. The importance of each feature in discriminating cancer group from control groups was then determined and the mean value of each feature was subsequently calculated. By several cycles of eliminating the least significant feature and progressively reducing the features, a combination with the smallest variance was obtained. This approach enhanced model performance during the feature selection process, capitalizing the strong interpretability of linear SVM and generalization performance of AdaBoost.

Figure 2.

The WGCNA analysis of protein-coding genes related with HERVs-elevation. (a) the correlation between gene sets identified by WGCNA and LTR16A|72|451 high (LTR16A_H) and low (LTR16A_L) expression status. (b) The correlations between genes and HERVs expression status within ME3 module, the horizontal axis represents correlation coefficients between genes and modules, while the vertical axis represents the correlation coefficient between genes and LTR high expression. (c) The cluster profiles of protein coding genes in the ME3 module in gastric cancer samples.

2.4HERVs-human protein coding genes co-expression network analysis

According to the selected HERV expression intensity, the samples were divided into high and low groups, and an correlation analysis between HERV and genes was performed using the WGCNA package (V1.71) [30]. The set of genes with a high correlation with the high and low expression of HERV were compiled and the correlation network was constructed. The pathway processes were analyzed using Clusterprofiler (V4.0.5) [31] to determine the biological function of the gene sets. Finally, we utilized the chooseTopHubInEachModule function from WGCNA to identify the genes with the highest connectivity within the entire gene set. Additionally, we employed gastric cancer data from TCGA to examine the impact of these genes on survival.

Figure 3.

Genomic positions and the correlation network of protein-coding genes associated with HERVs-elevation, From the outermost ring, the first ring represents chromosomes, the second ring represents the coordinates of genes, the third ring represents the expression status, and the central position connecting lines represent the weight relationships among genes based on WGCNA results. Red lines represent correlations within the same chromosome, blue represents interrelationships between chromosomes.

3.1Differential expression of HERVs in gastric cancer and screening of biomarkers

We compiled the transcriptome datasets of gastric cancers and processed the raw expression matrix using DESeq2. The significantly expressed genes were identified between cancer and normal groups (Fig. 1a), in which 22 HERVs were upregulated and 18 HERVs were downregulated.

To identify potential HERV markers for gastric cancer, we carried out a multivariate screening process using the SVM model and AdaBoost algorithm. The LTR16A|72|451, LTR27D|87|222, and LTR91|636|874 emerged as the intersections of the two algorithms (Fig. 1b). The original results and explanations of the two algorithms are provided in Supplementary Material 1. The LTR16A|72|451 and LTR27D|87|222 were significant upregulated in tumor subgroups, while LTR91|636|874 was significantly downregulated (Fig. 1c). Conversely, the expression of LTR16A|72|451 and LTR91|636|874 were not significantly different between lung cancer and normal tissues, and LTR27D|87|222 was significantly downregulated in lung cancer groups. This situation is demonstrated in Fig. 1d. Additionally, we validated the expression levels of three potential biomarkers using raw data collected from various research projects. In the gastric cancer validation set, the expression patterns were completely consistent with those observed in this study. However, in tumors with higher mortality rates, such as breast cancer, liver cancer, colon cancer, and pancreatic cancer, we did not observe any significant correlation. Due to space limitations in the article, we have presented specific expression matrices and box plots in Supplementary Material 2 and 3.

Figure 4.

Functional analysis of HERV-elevation related genes based on (a) KEGG pathway (b) GO molecular function domain (c) GO biological process domain (d) GO cellular component domain.

Figure 5.

GSEA of HERV-elevation related genes. (a) Profile of the running ES score and the gene set members on the rank ordered list. (b) The functional relationships among the HERV-elevation related pathways.

3.2Co-expression network analysis of HERV marks and protein-coding genes

To further investigate the relationship between HERV biomarkers and gene expression, we used the median expression levels of these three biomarkers as a standard for high and low grouping. Combined with the previous expression level box plots, we found that the biomarker LTR16A|72|451 exhibits significant differences between normal and tumor groups, with a noticeably higher overall expression level than the other two HERVs, making it more suitable for observing correlations with protein-expressing genes. Therefore, we used this biomarker for high and low grouping to study the correlation between its expression changes and the expression of genes encoding proteins.

Weighted genes co-expression network analysis (WGCNA) analysis was performed to determine the relationship between HERV expression and gene expression in the high and low expression groups, respectively. The Principal Component Analysis (PCA) of the gene expression data within a specific module generated the Module Eigengene (ME) and represented the overall expression characteristics of the genes. In brief, genes with similar expression patterns were grouped into modules, PCA was then performed on the gene expression data for each module. The first principal component captured the largest portion of variation in the gene expression data and was used as the Module Eigengene. The ME carried the overall expression pattern within the module and was utilized to evaluate the correlation between the module and phenotypic traits, in this way, the gene modules significantly associated with traits could be identified.

We displayed the correlation between the gene sets after clustering and traits in Fig. 2a, where a high correlation between the gene set represented by ME3 and the high expression of LTR16A is clearly visible. Subsequently, the correlation between all genes within the ME3 module and the condition of high expression of LTR16A|72|451 is shown in Fig. 2b, revealing a strong positive correlation between them. According to the clustering heatmap of gene expression, genes with the strongest correlation with high HERV expression displayed a notable distinction between normal and tumor groups (Fig. 2c).

According to the WGCNA results, we retrieved the interaction relationships between genes within the corresponding gene set (Fig. 3). We selected the weight greater than 0.4 as cutoff and constructed the gene correlation network. The network degrees were computed, identifying RAB31 as a central hub. This designation of RAB31 as the hub gene was achieved using the ‘chooseTopHubInEachModule’ function. This function operates by selecting genes within distinct modules and assessing the correlation among these genes based on their expression data. These correlations are subsequently converted into measures of connectivity strength, which are used to construct an adjacency matrix. Through this process, RAB31 was pinpointed as the gene exhibiting the highest connectivity within its module, establishing its role as a key hub.

3.3Functional analysis of protein-coding genes correlated with HERV expression elevation

To further explore the biological functions and pathways related by the genes with high correlation to HERV elevation, we performed KEGG and GO analysis on this gene set. The KEGG analysis revealed that the genes associated with the three potential HERV markers were mainly enriched in the PI3K-Akt pathway, which is critical to tumorigenesis (Fig. 4). We also performed Gene Set Enrichment Analysis (GSEA) to score the identified pathways using the foldchange information. The top five pathways with significant enrichment showed an overall upregulation trend. Among them, the IL-17 signaling pathway showed the most significantly enriched and followed by the JAK-STAT and TNF pathway (Fig. 5a).

Subsequently, we investigated the interrelationships of several major pathways in Fig. 5b. The genes associated with the gene set were mainly enriched in the IL-17 pathway, and several major pathways belong to typical inflammatory response pathways [32]. We can clearly see the process starting from the IL17 pathway, affecting TNF through genes such as MMP3 and CXCL1, and then extending to NF-kappa B through genes like CXCL1 and VCAM1. Meanwhile, IL17 also connects to the PI3K-Akt pathway through IL6 and AKT3, broadly affecting more related genes. In this process, we can observe several hallmark tumor pathways that influence and interact with each other, closely linked.

3.4Survival analysis of key genes RAB31 which correlated with HERV expression elevation

To investigate the functional impact of HERV-related genes on the survival of gastric cancer patients, we identified RAB31 as the hub gene within our selected gene set, proposing it as a candidate prognostic marker. The gene RAB31 in gastric cancer can also be a potential therapeutic target [33, 34]. We obtained a gastric cancer dataset from TCGA and used the median expression level of RAB31 as the basis for dividing into high and low expression groups. The results of the survival analysis based on these groups were then presented in Fig. 6.

Figure 6.

Survival analysis of the RAB31 gene in TCGA gastric cancer dataset. The patients were classified into RAB31-high and RAB31-low groups according to with median tpm of RAB31 across samples.