Novel significant stage-specific differentially expressed genes in hepatocellular carcinoma

Abstract

Background: Liver cancer is among top deadly cancers worldwide with a very poor prognosis, and the liver is a vulnerable site for metastases of other cancers. Early diagnosis is crucial for treatment of the predominant liver cancers, namely hepatocellular carcinoma (HCC). Here we developed a novel computational framework for the stage-specific analysis of HCC.

Methods: Using publicly available clinical and RNA-Seq data of cancer samples and controls and the AJCC staging system, we performed a linear modelling analysis of gene expression across all stages and found significant genome-wide changes in the log fold-change of gene expression in cancer samples relative to control. To identify genes that were stage-specific controlling for confounding differential expression in other stages, we developed a set of six pairwise contrasts between the stages and enforced a p-value threshold (< 0.05) for each such contrast. Genes were specific for a stage if they passed all the significance filters for that stage. The monotonicity of gene expression with cancer progression was analyzed with a linear model using the cancer stage as a numeric variable.

Results: Our analysis yielded two stage-I specific genes (CA9, WNT7B), two stage-II specific genes (APOBEC3B, FAM186A), ten stage-III specific genes including DLG5, PARI, NCAPG2, GNMT and XRCC2, and 35 stage-IV specific genes including GABRD, PGAM2, PECAM1 and CXCR2P1. Overexpression of DLG5 was found to be tumor-promoting contrary to the cancer literature on this gene. Further, GABRD was found to be signifincantly monotonically upregulated across stages. Our work has revealed 1977 genes with significant monotonic patterns of expression across cancer stages. NDUFA4L2, CRHBP and PIGU were top genes with monotonic changes of expression across cancer stages that could represent promising targets for therapy. Comparison with gene signatures from the BCLC staging system identified two genes, HSP90AB1 and ARHGAP42. Gene set enrichment analysis indicated overrepresented pathways specific to each stage, notably viral infection pathways in HCC initiation.

Conclusions: Our study identified novel significant stage-specific differentially expressed genes which could enhance our understanding of the molecular determinants of hepatocellular carcinoma progression. Our findings could serve as biomarkers that potentially underpin diagnosis as well as pinpoint therapeutic targets.

Keywords: Cancer progression; Differentially expressed genes; HCC stages; LIHC transcriptomics; Linear modelling; Metastasis; Monotonic expression; Pairwise contrasts; Significance analysis; Stage-specific biomarkers; Tumorigenesis.

Fig. 1

Major causative pathways of hepatocarcinogenesis. All pathways converge to progressive genomic alterations, leading a normal cell to acquire the hallmarks of cancer

Fig. 2

TCGA ‘Hybridization REF’ Barcode. The first 10 characters constitute an anonymized unique patient identifier and the following two characters denote whether the sample is tumor or normal

Fig. 3

Design matrices. a In the linear modeling, the control samples served as the baseline expression (intercept) of each gene against which the stage-specific expression was estimated. b the design matrix for the contrasts analysis

Fig. 4

A Venn representation of the pairwise stages contrasts. A gene could be differentially expressed in any combination of the four stages and this could be represented by a 4-bit string, one bit for each stage. For e.g., ‘1111’ at the overlap of all four stages would be assigned to genes that are differentially expressed in all four stages

Fig. 5

Boxplots of top 9 linear model genes. For each gene, notice that the trend in expression could be either overexpression or downregulation relative to the control. For e.g., GABRD, PLVAP, CXorf36, CDH13 and UBE2T are overexpressed, while CLEC4M, CLEC1B, BMP10, and CLEC4G are downregulated. It could be seen that a linear trend does not imply maximal |lfc| in stage 4, as illustrated most clearly in the case of UBE2T

Fig. 6

Boxplots illustrating stage-specificity of differentially expressed genes. Extremal expression in a stage could be either maximal expression or minimal expression relative to the control and all other stages, and could be termed maximal differential expression. Here we show genes with maximal differential expression in stage-I (WDR72; minimum expression), stage-II (GLI4, maximum expression; COLEC11, minimum expression), stage-III (CKAP2; maximum expression), and stage-IV (MAPK11; maximum expression)

Fig. 7

Principal components analysis of cancer vs control. a The first two principal components of the top 100 genes from linear modeling are plotted. It could be seen that control samples (red) clustered independent of the cancer samples (colored by stage). b The same analysis repeated with 100 random genes failed to effect a clustering of the control samples relative to the cancer samples

Fig. 8

Venn illustration of the size of each 4-bit string. The numbers of genes with each pattern of differential expression are shown

Fig. 9

Volcano plot of the 49 significant stage-specific differentially expressed genes. Stage 1 genes, red; Stage 2, blue; Stage 3, green; and Stage 4, orange. The genes are seen to orient away from the origin and the axes, indicating significance and effect size

Fig. 10

Heatmap plots of the final 24 stage-specific genes. a heatmap generated from the lfc values of all the stage-specific genes (arranged stagewise). The color gradient spans the spectrum from downregulation (blue) to overexpression (red). Log fold changes upto sixfold are seen, indicating 64 times differential expression with respect to control. b Representation of the stagewise gene expression based on clustering of differential expression profiles

Fig. 11

Boxplot of top genes with monotonic expression. These six genes (GABRD, PIGU, NDUFA4L2, CRHBP, FLVCR1, TTC13) showed monotonic trends of expression across the cancer stages, and were topranked in both the linear models given by eqns. (1) and (2)

Fig. 12

Boxplot of stage-I specific genes. It is seen that CA9 and WNT7B are both maximally downregulated in stage-I

Fig. 13

Boxplot of stage-II specific genes. It is seen that both APOBEC3B and FAM186A are maximally overexpressed in stage-II, the trend following an inverted U-shape

Fig. 14

Boxplot of stage-III specific genes. Except for GNMT, the expression of stage-III specific genes show a peak in stage-III, with the expression trend following an inverted U-shape across the stages. The expression trend is convex and reversed for the downregulated GNMT, with minimum expression in stage-III

Fig. 15

Boxplot of top 10 stage-IV specific genes. All genes, except NR1I2 and CXCR2P1, show a smooth increasing expression trend reaching peak expression in stage-IV. In the case of NR1I2 and CXC2RP1, the trend is reversed, with the expression decreasing smoothly to touch the minimum in stage-IV