Screening of differentially expressed genes
Download RNA-seq data using the gene expression of edge R screening difference, based on fold change = 1 (P < 0.05) screened to 9314 differentially expressed genes DEGs were (wherein 782 genes were up-regulated as a threshold in a tcga; 1493 genes were downregulated ). The volcano map made by R software visualizes these differential genes, where logFC ≥ 1 or logFC≤-1, and the error rate FDR < 0.05 (Fig. 2). A heat map is used to show the expression of 50 genes in 50 normal tissue samples and 373 cancer tissue samples (Fig. 3)
Construction of co-expression networks
Based on the expression data of 9314 differentially expressed genes in 423 liver cancer tissues, the WGCNA algorithm was used to construct the co-expression module. The function hclust was used for sample clustering analysis (Figure. 4a), and the weighted co-expression network was constructed, and then the heat maps of different clinical features were drawn (Figure. 4b). The network construction needs to perform network topology analysis, and the standard soft threshold power value of scale-free network was determined by R calculation (Figure. 5), and then the gene modules were divided, and different modules were marked with color (Figure. 6). Network visualization, randomly select 10,000 genes for mapping, you can clearly observe the color difference between the modules, and graphically show the relationship between the modules and the traits (Figure. 7a, b).
Identification of modules
Liver cancer used clinical information about the patients downloaded from the TCGA, trait association analysis module, and found clinical traits significant correlation module (Figure.8). Pick the modules (darkgrey and navajowhite2) that are most relevant to TNM staging. Therefore, the module can calculate the correlation coefficient (MM) with the gene, and all the continuous traits can also calculate the correlation coefficient (GS) with the gene expression value. Finally, the two correlation matrices are combined, and the modules of interest are formulated for analysis. After comparing the scatter diagrams of each module, it can be found that the genes in the darkgrey and navajowhite2 gate modules are highly correlated with traits and modules at the same time (Figure. 9a, b).
Functional analysis of gene module
The GO enrichment analysis was carried out on the modules darkrey and navajowhite2, which are highly related to TNM, and the main molecular functions of the genes in the modules were observed. The results showed that the darkgyrey module genes mainly act on the reduction reaction process; The navajowhite2 modular gene mainly acts on the activation of specific DNA-binding transcription factors (Fig. 10a, b). These functions may be related to liver cancer metastasis.
Hub genes in important modules
The darkgyy module and navajowhite2 were introduced in Cytoscape, and the darkgyy module CytoHubba algorithm was used to screen the Top30 genes in BottleNect. Analyzed RP11-109J4.1, TXNRD1, C9orf106, LINC00677, CYP11B2, AC005550.3 belong to the hub gene, in the navajowhite2 module CytoHubba DMNC algorithm screening analysis of Top10 genes, and found that ABLIM2 and MIR3945HG belong to the central gene (Fig. 11a, b).
Identification and evaluation of prognostic models
Based on the clinical survival information of liver cancer samples, Cox regression was applied to two modules (darkgrey and Navajowhite2) respectively. AC005550.3, C9orf106, RP11-109J4.1, CYP11B2, LINC00677, TXNRD1, AKR1B15, SLC5A11, RP11-818F20.5, RP11-276E17.2, CTD-3037G24.4 were identified in the module darkgrey(Table 1). C9orf106, CYP11B2, LINC00677 are hub genes. Module Navajowhite2 identified MIR394HG as hub gene (Table 2). According to the expression level of these genes, the risk score of each sample was established. According to the median risk score, liver cancer samples were divided into high-risk group and low-risk group. According to the Kaplan Meier survival curve, the high-risk group had worse overall survival (OS), indicating that the prognostic characteristics of risk score were obvious (Figure. 12a, b). Heat map of gene expression profile in darkgrey in liver cancer samples (Figure. 13)
Table 1
genes of top30 based on BottleNeck, Closeness, Degree, DMNC, EcCentricity, EPC,MCC,MNC darkgrey modules.
Table 1
|
|
|
Term
|
Count
|
Genes
|
BottleNeck
|
30
|
RP11-143E21.3 SHANK1 CYP11B2 CYP1B1 MEOX2-AS1 RP11-20G13.2 SOX14 GRM5 RP11-65J3.2 TRPC5OS RP11-20G13.3 AHRR RP11-91H12.4 STC2 EGF SRXN1 RP11-818F20.5 AKR1B10P1 RP11-91H12.1 LINC01483 NTNG2 PNLIP UGT2B26P RP11-65J3.3 GRM5-AS1 LINC00677 KRT34 LHX9 MSC-AS1 PIEZO2
|
Closeness
|
30
|
RP11-143E21.3 MEOX2-AS1 CYP11B2 STC2 CYP1B1 AHRR GRM5 LINC00677 SRGAP3-AS2 RP11-20G13.3 RP11-65J3.2 SOX14 RP11-311F12.2 RP11-91H12.4 RP11-20G13.2 RP11-146F11.5 RP11-91H12.1 SHANK1 RP5-881P19.7 RP11-65J3.3 GRM5-AS1 RP11-818F20.5 TRPC5OS NTNG2 CDKL4 RP11-383J24.1 RP11-766F14.2 ALDH3A1 ANKFN1 LINC01483
|
Degree
|
30
|
RP11-143E21.3 MEOX2-AS1 STC2 GRM5 CYP11B2 AHRR CYP1B1 SRGAP3-AS2 RP11-91H12.4 RP11-65J3.2 RP11-20G13.3 LINC00677 SOX14 RP11-311F12.2 RP11-20G13.2 SHANK1 RP11-146F11.5 TRPC5OS RP11-65J3.3 ANKFN1 RP11-818F20.5 RP11-91H12.1 RP5-881P19.7 GRM5-AS1 RP11-383J24.1 CDKL4 NTNG2 ALDH3A1 LINC01483 RP11-766F14.2
|
DMNC
|
30
|
GRM5-AS1 RP11-766F14.2 LINC00677 CYP1B1 RP11-65J3.2 RP11-91H12.1 RP11-143E21.3 GRM5 RP11-65J3.3 ALDH3A1 EGF RP11-146F11.5 KRT34 PIEZO2 RP11-449J21.5 MRAS RP5-881P19.7 CDKL4 C9orf106 STC2 MEOX2-AS1 SOX14 SRGAP3-AS2 TRPC5OS SHANK1 RP11-383J24.1 RP11-818F20.5 NTNG2 RP11-20G13.3 RP11-311F12.2
|
EcCentricity
|
30
|
GRM5 RP11-65J3.2 GRM5-AS1 RP11-20G13.3 CYP11B2 LINC00677 AHRR CYP1B1 SOX14 RP11-143E21.3 STC2 EGF RP11-146F11.5 MEOX2-AS1 PNLIP RP5-1182A14.5 UGT2B26P TRPC5OS RP11-65J3.3 SHANK1 RP11-20G13.2 KRT34 PIEZO2 CLRN1-AS1 STK32B ALDH3A1 RP11-311F12.2 GLYATL1P3 RP11-91H12.4 SLC13A3
|
EPC
|
30
|
RP11-143E21.3 CYP1B1 MEOX2-AS1 STC2 CYP11B2 LINC00677 RP11-65J3.2 AHRR SRGAP3-AS2 RP11-91H12.4 GRM5 RP11-20G13.3 SOX14 RP11-20G13.2 RP11-311F12.2 RP11-91H12.1 SHANK1 RP11-146F11.5 RP5-881P19.7 TRPC5OS RP11-383J24.1 CDKL4 RP11-818F20.5 RP11-65J3.3 GRM5-AS1 NTNG2 RP11-766F14.2 ALDH3A1 MRAS ANKFN1
|
MCC
|
30
|
CYP1B1 RP11-143E21.3 MEOX2-AS1 RP11-65J3.2 LINC00677 GRM5 STC2 SRGAP3-AS2 CYP11B2 RP11-20G13.3 AHRR RP11-91H12.4 SOX14 RP11-20G13.2 RP11-311F12.2 RP11-91H12.1 SHANK1 TRPC5OS GRM5-AS1 RP11-818F20.5 RP11-146F11.5 RP11-65J3.3 RP11-383J24.1 RP11-766F14.2 NTNG2 ANKFN1 RP5-881P19.7 ALDH3A1 CDKL4 EGF
|
MNC
|
30
|
RP11-143E21.3 STC2 MEOX2-AS1 CYP11B2 CYP1B1 GRM5 AHRR RP11-91H12.4 RP11-65J3.2 RP11-20G13.3 LINC00677 SRGAP3-AS2 RP11-20G13.2 SOX14 RP11-311F12.2 TRPC5OS SHANK1 RP11-383J24.1 RP11-818F20.5 RP11-91H12.1 NTNG2 RP11-65J3.3 GRM5-AS1 ALDH3A1 EGF RP11-146F11.5 RP11-766F14.2 KRT34 PIEZO2 RP11-449J21.5
|
Table 2
genes of top30 based on BottleNeck, Closeness, Degree, DMNC, EcCentricity, EPC,MCC,MNC navajowhite modules.
Table 2
|
|
|
Term
|
Count
|
Genes
|
BottleNeck
|
30
|
ATP10B ARR3 FMO8P NDUFA5P6 UMODL1 CEACAM20 TRIM40 RP11-295G20.2 SLC25A15P3 C5orf67 RP5-965G21.4 NTS RP11-443B7.2 RP11-756A22.7 ABLIM2 DUOXA2 TRIM31-AS1 PRAC2 CTD-2091N23.1 RP4-607I7.1 GPC6-AS1 LTF CTA-363E6.6 WNT10B LINC00332 WNT1 RP11-752D24.2 ALPPL2 MIR3945HG RP11-85I21.1
|
Closeness
|
30
|
ATP10B NDUFA5P6 FMO8P UMODL1 ABLIM2 C5orf67 CEACAM20 WNT1 RP11-295G20.2 SLC25A15P3 RP4-607I7.1 CTA-363E6.6 RP11-752D24.2 LTF RP11-443B7.2 RP11-756A22.7 ARR3 RP5-965G21.4 NTS DUOXA2 TRIM40 TRIM31-AS1 WNT10B PRAC2 CTD-2091N23.1 GPC6-AS1 LINC00332 ALPPL2 MIR3945HG RP11-85I21.1
|
Degree
|
30
|
NDUFA5P6 ATP10B ABLIM2 ARR3 FMO8P C5orf67 UMODL1 RP11-295G20.2 CEACAM20 SLC25A15P3 WNT1 TRIM40 RP5-965G21.4 NTS RP11-443B7.2 RP11-756A22.7 DUOXA2 TRIM31-AS1 PRAC2 CTD-2091N23.1 RP4-607I7.1 GPC6-AS1 LTF CTA-363E6.6 WNT10B LINC00332 RP11-752D24.2 ALPPL2 MIR3945HG RP11-85I21.1
|
DMNC
|
30
|
ALPPL2 RP11-295G20.2 MIR3945HG RP11-85I21.1 RP5-965G21.4 CEACAM20 NTS RP11-443B7.2 TRIM31 RP11-756A22.7 KISS1 ABLIM2 ARR3 TMEM74 RP11-227L6.1 SLC25A15P3 SPINK1 DUOXA2 TRIM31-AS1 C19orf33 PRAC2 CTD-2091N23.1 RP4-607I7.1 GPC6-AS1 AF127577.8 LTF U95743.1 BPIFB1 RP11-66H6.4 CTA-363E6.6
|
EcCentricity
|
30
|
ATP10B CEACAM20 ARR3 FMO8P NDUFA5P6 UMODL1 RP11-295G20.2 ABLIM2 SLC25A15P3 RP4-607I7.1 CTA-363E6.6 C5orf67 WNT1 RP11-752D24.2 TRIM40 RP11-443B7.2 RP11-756A22.7 LTF RP5-965G21.4 NTS DUOXA2 PRAC2 CTD-2091N23.1 GPC6-AS1 LINC00332 TRIM31-AS1 WNT10B ALPPL2 MIR3945HG RP11-85I21.1
|
EPC
|
30
|
ATP10B UMODL1 FMO8P ABLIM2 NDUFA5P6 C5orf67 WNT1 CEACAM20 RP11-295G20.2 SLC25A15P3 RP11-752D24.2 LTF CTA-363E6.6 RP4-607I7.1 ARR3 RP11-443B7.2 RP11-756A22.7 NTS RP5-965G21.4 TRIM40 DUOXA2 WNT10B TRIM31-AS1 PRAC2 CTD-2091N23.1 GPC6-AS1 LINC00332 ALPPL2 MIR3945HG RP11-85I21.1
|
MCC
|
30
|
NDUFA5P6 ATP10B ABLIM2 ARR3 FMO8P C5orf67 UMODL1 RP11-295G20.2 CEACAM20 SLC25A15P3 WNT1 TRIM40 RP5-965G21.4 NTS RP11-443B7.2 RP11-756A22.7 DUOXA2 TRIM31-AS1 PRAC2 CTD-2091N23.1 RP4-607I7.1 GPC6-AS1 LTF CTA-363E6.6 WNT10B LINC00332 RP11-752D24.2 ALPPL2 MIR3945HG RP11-85I21.1
|
MNC
|
30
|
RP11-295G20.2 RP5-965G21.4 CEACAM20 NTS RP11-443B7.2 RP11-756A22.7 ABLIM2 ARR3 SLC25A15P3 DUOXA2 TRIM31-AS1 PRAC2 CTD-2091N23.1 RP4-607I7.1 GPC6-AS1 LTF CTA-363E6.6 FMO8P WNT10B C5orf67 LINC00332 NDUFA5P6 WNT1 ATP10B UMODL1 RP11-752D24.2 TRIM40 ALPPL2 MIR3945HG RP11-85I21.1
|
ROC curve of survival model
Univariate Cox regression identified 26 potential predictors (SLC5A11, MPZ, RP11-276E17.2, VASH2, TRIM16, RP11-552C15.1, AKR1B15, CYP11B2, RP11-109J4.1, MAFG, RP11-311F12.2, LINC00677, TRIM16L, AKR1B10P1, UNKL, RP11-65J3.2, CTD-3037G24.4, RP11-818F20.5, TXNRD1, C9orf106, GRM5, LINC00635, ALDH3A1, AKR1B10, AC005550.3, XKRX) Multivariate Cox regression analysis showed that six genes (RP11-109J4.1, TXNRD1, C9orf106, LINC00677, CYP11B2, AC005550.3) and one gene (MIR3945HG) in darkgrey were independent prognostic factors for liver cancer (Figure. 14a, b). In order to compare the sensitivity and specificity of risk score for prognosis of patients with hepatocellular carcinoma, time-dependent receptor operating characteristic (ROC) analysis was performed. The area under the ROC curve (AUC) was 0.8, 0.764, 0.76387 for the 1-year, 3-year, 5-year survival risk score of darkgrey module and 0.665 for the 1-year, 3-year, 5-year survival risk score of Navajowhite2 module (Figure. 15a, b), indicating RP11-109J4.1, TXNRD1, C9orf106, LINC00677, CYP11B2, AC005550.3 and MIR3945HG.The seven genes are highly reliable for the prognosis of liver cancer.
Prognostic analysis of key genes
AC005550.3, C9orf106, CYP11B2, LINC00677, RP11-109J4.1, TXNRD1, MIR3945HG are imported into the GEPIA database. Except MIR3945HG is not included in the database, all others have records. First, analyze the expression levels of genes such as AC005550.3, C9orf106, CYP11B2, LINC00677, RP11-109J4.1, TXNRD1. The liver cancer and threshold (log2FC and q-value) and the difference analysis method all select the default value, input a single gene one by one, analyze the pancancer of the gene, and get the expression map with liver cancer (Figure. 16a, b, c, d, e, f).Then the expression analysis of the above genes in normal tissues and tumor tissues in liver cancer is compared, and the expression box plot is obtained. (Figure. 17a, b, c, d, e, f). After that, the expression difference of the above-mentioned genes in different pathological analyses in liver cancer is analyzed to understand the expression changes of a certain gene with the change of the disease. LINC00677, CYP11B2,AC005550.3 were highly expressed in stage I, stage II, stage III.RP11-109J4.1,C9orf106 and TXNRD1 were highly expressed in all stages (Figure. 18a, b, c, d, e, f). Finally, the survival analysis of AC005550.3, C9orf106, RP11-109J4.1, TXNRD1 genes in liver cancer was analyzed (Figure. 19a, b, c, d).
Infiltrating immune correlation of key genes
The time database was used to evaluate the correlation between prognosis gene expression and immune invasion. The results showed that C9orf106, CYP11B2, TXNRD1 were related to immune invasion (Figure. 20a, b, c).