2.1 Identification of differentially expressed ATGs
316 OSCC patients and 32 normal controls RNA-seq and corresponding clinical data were downloaded from TCGA database. Then, 232 autophagy related genes expression level were extracted from the transcriptome data. Subsequently, differentially expressed analysis was performed in R software EdgeR package(Figure1A,B,C).Finally, 11 upregulated and 26 downregulated autophagy related genes were sorted out With the cut-off criteria |log2FoldChange|>1 and FDR < 0.01.
2.2 Functional Analysis of differentially expressed autophagy related genes
GO enrichment analysis indicated that these genes were mainly located in autophagosome membrane, and associated with cytokine receptor binding and regulation of apoptotic signaling pathway(Figure2A,C). In addition, KEGG pathway analyses showed that most of significant autophagy related genes were enriched in apoptosis, platinum drug resistance, ErbB signaling pathway and TNF signaling pathway(Figure2B,D).
2.3 Establishment of cox regression model
Through univariate and multivariate cox regression analysis, total of 9 variables including BID, ATG12, BAK1, SPHK1, NKX2−3, ATIC, LAMP1, ATF6, BNIP3 were enrolled in cox model(Figure3A,B). And rish_score=(0.35502*BID)+(0.69633*ATG12)+(0.22561*BAK1)- (0.24922*SPHK1)-(0.66016*NKX2-3)+(0.30945*ATIC)+(0.30416*LAMP1)+(0.50726*ATF6) +(0.26573*BNIP3). Subsequently, OSCC patients in TCGA database were divided into high risk group and low risk group according to cox formula median.Survival analysis indicated that overall survival of high risk core group patients were significantly lower than the low-risk group(Figure3C). Moreover, the expression level of protective ATGs in the low risk group was higher than that in the high risk group.On the contrary, risky genes were lower in high risk group (Figure3D). The risk score combined with survival data were visualized in R software (Figure3E,F).
2.4 Identification of cox regression model
Firstly, the ROC was plotted, and its area under the curve(AUC) is 0.75 which markedly higher than other clinical characteristics(Figure4A). Furthermore, risk score in early stage was significantly lower than terminal stage(Figure4B) indicating that risk score on basis of ATGs may realize early diagnosis in OSCC. Moreover, univariate and multivariate cox regression analysis indicating that risk score may be regarded as an independent prognostic factor(Figure4C,D). Furthmore, the survival curve of ATGs relevant to OSCC overall survival were ploted in R software. And ATG12 and BID were selected as 2 potential independent biomarkers of OSCC in the cox regression model according to univariate(Figure3A), multivariate cox regression(Figure3B) and survival analysis (Figure4E,F) .
2.5 Survival analysis
The prognostic values of the risk score for different clinicopathological parameters including age, gender, T and N in TNM system, grade and stage were further investigated. M classification in TNM system were excluded because of numerous data missing. Survival analysis demonstrated that low risk group had significantly longer overall survival than high risk group in the stratification analysis based on age, gender, T, N, grade and stage(Figure5).
2.6 Identification of potential independent prognostic biomakers
Comprehensive bioinformatics indicated that ATG12 and BID were associated with overall survival in OSCC. And univariate and multivariate cox regression showed that ATG12 and BID might be selected as potential independent prognostic biomarkers in our study. Therefore, ATG12 and BID genes expression were validated in OSCC cell lines and tissues by qRT-PCR. Our results revealed that ATG12 and BID were upregulation in OSCC cell lines and tissues than MNTs, which were similar with the results in TCGA database (Figure6A,B). Unfortunately, the correlation between ATG12,BID and clinical parameters showed no significant difference (Table1, Table2).