Identification of Biomarker for Cutaneous Squamous Cell Carcinoma Using Microarray Data Analysis.

Cutaneous squamous cell carcinoma (CSCC) is one of the most malignant tumors worldwide. We aimed to explore the molecular mechanism of this CSCC and screen feature genes that can function as the biomarker of CSCC and thus provide a theoretical basis for the pathogenesis research and development of medicine. The method of microarray data analysis was used in this study to explore the differentially expressed genes between tissues of normal specimens and tissues of patients with CSCC. Besides, functional enrichment analysis and signal pathway were performed on these genes to screen the feature genes that are closely associated with CSCC can function as the potential biomarkers of CSCC.A total of 53 samples from two datasets, GSE45216 and GSE45164, were used in the differentially expressed analysis. And as a result, a total of 833 genes were screened out, including 465 up-regulated genes and 215 down-regulated genes. Candidate genes, including up-regulated genes like S100A12, MMP1, DEFB4B/DEFB4A, KRT16 and PI3, and down-regulated genes like EGR3, LRP4, C14orf132, PAMR1, CCL27, and KRT2 were screened out. All these genes were testified in the dataset of GSE66359. The result showed that only three genes, KRT16, PI3 and EGR3, were mostly differentially expressed and only EGR3 had the same expression pattern with both datasets, GSE45216 and GSE45164.Of note, EGR3 gene was found to be the most differentially expressed gene in cutaneous squamous cell carcinoma, which had the potential to function as the candidate genes and help in the diagnosis and prognostic treatments of CSCC.

Introduction

Cutaneous squamous cell carcinoma (SCC) is the second most common malignant tumor in the world, just after basal cell carcinoma 1. Nowadays, more and more people are affected by this disease and white skin patients with immunologic suppressionor chronic skin inflammation disease are more apt to be affected. Besides, it's found that middle-aged and elderly people are relatively more susceptible to this disease compared to younger people. The occurrence of cutaneous squamous cell carcinoma is related to multiple factors. Besides, its diversity clinically and pathologically makes it easy to escape diagnosis or be misdiagnosed 2. Currently, early diagnosis and prediction and timely treatment are still the most effective measure to improve the survival rate of patients with CSCC and prevent the disease from deteriorating. However, the research on biomarkers of squamous cell carcinoma was rare. Therefore, it is urgent to explore potential biomarker of cutaneous squamous cell carcinoma which can be quite beneficial for improving the clinical management of CSCC.

DNA microarray data is commonly used in clinical research since it can monitor expression levels of thousands of genes at the same time. And it had become a promising and prevailing method used in the identification of differentially expressed genes between normal samples and tumor specimens 3.What's more, this method provides a complete, systematic, and reliable comparison of gene expression between tissue types 4, 5. Therefore, biological processes and signaling pathway associated with the tumor may be explored by the evidence and clue provided by the DEGs. Recently, several studies indicated that the differentially expressed genes rapid and successful identification using DNA microarray data, such as in the human gliomas 6 and prostate cancer 7, etc.

The main purpose of this study is to explore the molecular mechanism of CSCC and to screen out the feature genes that can function as the biomarker of CSCC and thus provide a theoretical basis for the pathogenesis research and development of medicine. In this study, we came to the conclusion that EGR3 was closely associated with the occurrence of cutaneous squamous cell carcinoma and had the potential to function as the biomarker of CSCC, which could be quite helpful in the diversity clarification of CSCC clinically and histopathologically, therefore helpful in the improvement of prognostic and diagnostic tools and treatment of CSCC in clinical management.

Materials and Methods

Microarray data source

Microarray data was downloaded from GEO (Gene Expression Omnibus) dataset at the website of https://www.ncbi.nlm.nih.gov/geo/. Three separate datasets with the accession number of GSE45216 8, GSE45164 9, and GSE66359 10, 11 were selected for the analysis. There were 30 tumor samples in the GSE45216 dataset and we selected 10 samples with the same background and testing platform in the GSE42677 dataset as the control. GSE45164 and GSE45216 were used in the screening of critical genes. There were three normal specimens and 10 tumor samples in the GSE45164 dataset and there were five normal keratinocyte cells and 8 cutaneous squamous cell carcinoma cells. The testing platforms were Affymetrix Human Genome U133A 2.0 Array and Affymetrix Human Genome U133 Plus 2.0. GSE66359 was used in the testing of gene expression and the platform was Affymetrix Human Genome U133 Plus 2.0 Array.

Data quality control

AffyPLM package 12 was used for the data quality analysis based on the linear model at microarray level. RLE (Relative Log Expression) box figure and NUSE (Normalized Unscaled Standard Errors) figure were painted test the tendency accordance of the testing data. Besides, degradation situation of RNA was tested by AffyRNAdeg function. Finally, high-quality RNA datasets with the same tendency were selected for down-stream analysis.

Data pre-processing

Gcrma package 13 was used in the normalization and background correction of microarray data to ensure the integrity and comparability of the dataset. The error in the microarray and among was eliminated and average value was calculated for the genes that were tested for more than once and the value was used for the down-stream analysis.

Relevance analysis of gene expression level among samples was an essential indicator in the testing of experimental reliability and rationality in sample selection. Therefore, the global and principle component analysis were performed on the testing samples, and the Pearson correlation coefficient was calculated and correlation and distribution figure were painted.

Identification of differentially expressed genes

Limma packages 14 were used in the identification of the differentially expressed genes between samples in the control group and processed samples. Finally, differentially expressed genes with the log 2 value (fold change) larger than 1 and the p-value less than 0.05 were screened out. then, ggplot2, VennDiagram and heatmap in R language were used to paint the volcano, Venn figure and heatmap for the visualization of differentially expressed situation.

Functional enrichment analysis of differentially expressed genes

Pathways and functional enrichment analysis were performed on these differentially expressed genes. Functional Annotation Tool of DAVID (The Database for Annotation, Visualization, and Integrated Discovery) 15 were used in functional annotation and GO/KEGG enrichment analysis. P value before and after correction (Benjamini correction or FDR correction) was calculated by DAVID and p-value less than 0.05 was used as the threshold.

Results

Quality control of data

Regression calculation was used on the raw data by affyPLM in R language. Relative logarithm (RLE) expression box figure (Supplementary Figure ) and standard deviation chart were painted to test the homogeneity among microarray data. The result of RLE chart showed that the gene expression value of most samples in GSE45216 and GSE45164 datasets were enriched at zero with high accordance, indicating their feasibility for the use of down-stream analysis.

Normalization of data

Gcrma package was used for the normalization of samples after the quality screening. Theresult was shown in Figure . The result of expression density curve and box plot showed that the expression value in the two groups ranged from 0 to 15, with light variation, which corresponded with reality situation. The expression value of GSE45216 and GSE45164 after normalization was focused at 3. The expression trends of the two sets of microarray datasets are similar.

Cor function in R language was used on the expression data after normalization and logarithmic transformation to calculate the pearson correlation coefficient among samples and the correlation coefficient figure was painted subsequently (Figure ). The minimum value of correlation coefficient in GSE45216 dataset was 0.726 and normal samples and tumor specimens were divided into a different group. The minimum value of correlation coefficient in GSE45216 was 0.832, higher than that of GSE45216 dataset. There was a significant difference in normal tissue samples but the tumor samples can be divided into the same type. Besides, the result of principal component analysis (PCA) had the consistent results (Figure ).

Differentially expressed genes identification of microarray data after normalization was performed by Limma package. A total of 8032 differentially expressed genes were identified by the comparing the 30 tumor tissues and 10 normal tissues in the GSE45216 dataset, including 3474 up-regulated genes and 4560 down-regulated genes. The number of genes with the annotation information was 3118 (Figure ). A total of 1750 differentially expressed genes were identified by the comparing the 10 tumor tissues and 3 normal tissues in the GSE45216 dataset, including 956 up-regulated genes and 794 down-regulated genes. The number of genes with the annotation information was 1678 (Figure ). The number of Common genes in both GSE45216 and GSE45164 was 833, and genes with the same tendency were 680, including 465 up-regulated genes and 215 down-regulated genes (Figure ).

Functional enrichment analysis

Functional enrichment analysis was performed on these common genes by DAVID software (Figure ). There were 207 mostly enriched terms in the GO functional enrichment analysis. The maximum number of terms enriched in biological processes was 130 and the minimum number of terms enriched in molecular functions was 35. The top three functional terms enriched in molecular functions were 2'-5'-oligoadenylate synthetase activity, CXCR3 chemokine receptor binding, and cAMP-dependent protein kinase regulator activity. The mostly enriched terms enriched in biological processes were the maintenance of centrosome location, negative regulation of stress-activated MAPK cascade, positive regulation of ATP biosynthetic process and response to interferon-beta. The top three functional terms enriched in cell components were condensin complex, cell projection membrane, and meiotic spindle. A total of 38 signaling pathways were mostly enriched in the KEGG enrichment analysis, in which Leishmaniasis was the mostly enriched one.

Validation analysis of independent samples

The gene sets with fold change ranging the top 20 in the up and down-regulated genes in the two datasets were selected. As a result, up-regulated genes of S100A12, MMP1, DEFB4B///DEFB4A, KRT16 and PI3, and the down-regulated genes EGR3, LRP4, C14orf132, PAMR1, CCL27, and KRT2 were obtained and selected as the candidate genes for further analysis.The GSE66359 dataset, including 13 samples (8 tumor tissues and 5 normal tissues), at the platform of Affymetrix Human Genome U133 Plus 2.0, was used as the testing sample. The result of independence testing analysis showed that genes of KRT16, PI3 and EGR3 had the significant reduction of gene expression in CSCC tumor cells, but KRT16 and PI3 genes had the contrary to expected results (Figure ).

Discussion

Squamous Cell carcinoma (SCC), also known as epidermoid carcinoma, is a kind of malignant tumor that often occurs in epidermal cells or cells of appendages 16. The cancer cells are featured by different degree of keratosis. Organs covered by squamous epithelium, such as skin, mouth, lips, esophagus, cervix, vagina, etc, are more feasible to be affected. Some other organs, like bronchus, bladder and renal pelvis, though without the cover of squamous epithelium, can also transform to SCC by squamous metaplasia phenomenon. Currently, the incidence of SCC has been increased dramatically worldwide 17. Therefore, early diagnosis and treatment are still the most efficient measures to avoid the occurrence and development of the tumor. Nowadays, the biomarker research of SCC was rare and the identification of biomarker of this disease was urgent, which may quite helpful in the predictive and prognostic treatment of cancer in clinical research.

In this study, a total of 833 co-expression genes were identified by the differentially expressed genes analysis on the total of 53 genes from three datasets, that is, GSE45216, GSE42677, and GSE45164. As a result, a total of 680 genes with similar expression trend were figured out, including 465 up-regulated genes and 215 down-regulated genes. Furthermore, the candidate genes were verified in the GSE66359 dataset, which made the result more reliable. Finally, three genes, KRT16, PI3, and EGR3 were the most differentially expressed in cutaneous squamous cell carcinoma, but KRT16 and PI3 genes had the opposite expression pattern with those in GSE42677 and GSE45164 dataset. Unlike previous study that identify biomarker of CSCC in two contexts, in vitro and in vivo, in this study, three datasets, GSE45216, GSE42677, GSE45164, were selected for the research and the results were verified in GSE66359 dataset, precluding the overlapping of the context and made the results more reliable. A total of 53 samples, including 30 tumor samples in GSE45216 dataset, combined with 10 normal tissues in GSE42677 dataset after normalization and 10 tumor tissues and 3 normal tissues in GSE45164 dataset, were selected for the differentially expressed analysis, aiming to explore the molecular mechanism of squamous cell carcinoma. As a result, a total of 833 co-expression genes as KRT16 and PI3 etc, were screened out. The number of genes with the same tendency was 680, including 465 up-regulated genes and 215 down-regulated genes. Up-regulated genes of S100A12, MMP1, DEFB4B///DEFB4A, KRT16 and PI3, and down-regulated genes as EGR3, LRP4, C14orf132, PAMR1, CCL27, and KRT2 were identified as the potential target genes. Besides, the result was testified in GSE66359 dataset. The result showed that three genes, KRT16, PI3 and EGR were the most significant differentially expressed but the expression pattern of KRT16 and PI genes were contrary to that in GSE42677 and GSE45164 datasets.

EGR3 (Early Growth Response 3) gene belongs to EGR family and there were 4 genes in this gene family, that is EGR1~EGR4. There was a highly conserved DNA-binding domain encoding zinc finger proteins in the EGR family. This kind of gene had one reaction component that can bind to the zinc finger structure. Proteins encoded by EGR3 gene were the transcription factor of Cys2His2 zinc finger structure, which was also one type of genes responding to early growth process[18, 19]stimulated by karyomitosis. It was reported in early research that EGR2 and EGR3 genes stimulated NFκB and MAPK signaling pathway in the upstream 20 with the help of breast adipose fibroblasts TNFα. It played essential roles in fibrotic response and up-regulated in scleroderma disease 21. Genes in EGR family were down-regulated or usually absent in oral squamous cell carcinoma tissues and some squamous cell lines. However, the expression and significance of EGR3 in cutaneous squamous cell carcinoma haven't been made clear. Liao's study showed that the expression of EGR3 gene played a critical role in the differentiation, proliferation, metastasis and progression of gastric cancer cells 22 and Inoue's study addressed the role of Egr3 as an intracellular mediator of the estrogen-signaling pathway in breast cancer. All this study demonstrated the close association between EGR3 and cancer 23. In this study, EGR3 was found to be highly enriched and had exactly the same expression pattern with GSE45216 and GSE45164. All these results showed that EGR3 may be closely associated with the occurrence of squamous cell carcinoma and may function as the potential biomarker of this disease, suggesting a potential application in the improvement of prognostic tools and treatments of this disease.

Unlike conventional method, in this study, the candidate genes identified were further verified further in the GSE66359 dataset, which made the result more reliable. Though three candidate genes, that is, KRT16, PI3, and EGR3 were found to be mostly enriched, only EGR3 was found to have similar expression pattern with that in the GSE4267 and GSE45164 dataset. What's more, the association of EGR family genes with squamous cell carcinoma had been supported by a series of early studies. Therefore, we had reason to believe that EGR3 may be closely related to the occurrence and development of cutaneous squamous cell carcinoma and may function as the biomarkers of this disease.

Supplementary figures.

Click here for additional data file.