Transcriptomic insights into skin cancer: A bioinformatics and network biology approach to biomarker identification

2.1 Data Collection

We obtained Transcriptome-Seq data for human skin cancer and adjacent normal tissues from the ENA database in raw format for further analysis (Leiter et al., 2008). These samples were sequenced using the Illumina platform. The data is associated with Project ID PRJNA546533 and comprises a total of 12 datasets. The downloaded data files were in fastq.gz format, which is a text-based format used for storing both the sequenced biological sequence and its corresponding quality scores.

2.2 Transcriptome data analysis

The initial step in the Transcriptome-Seq data analysis protocol involves assessing the quality of raw data to identify and address any poor-quality reads that might impact subsequent analyses and interpretation. We utilize the FastQC Tool, specifically optimized for Illumina reads, which evaluates key parameters such as total base quality, per-tile quality, GC content, and sequence length distribution. Trimmomatic eliminates duplicated reads and removes low-quality reads to enhance data quality (de Sena Brandine and Smith, 2019).

2.3 Mapping and alignment of Transcriptome data

Following quality control, high-quality data undergo mapping and alignment to the human reference genome GRCh38/hg38 obtained from the ENSEMBL database. BOWTIE2, a Linux-based software, efficiently aligns reads to the reference genome, utilizing indexed files for accelerated alignment (Birney et al., 2004).. Paired-end data in fastq.gz format is inputted, generating a Sequence Alignment Map (SAM) file (Langmead and Salzberg, 2012). SAM files are converted to Binary Alignment Map (BAM) format using SAMTools for efficient storage and processing (Li et al., 2009). To optimize memory usage for subsequent analyses, the BAM alignment files are sorted.

2.4 Transcript Assembly and Quantification

The process of estimating gene and transcript expression involves assembling the mapped reads to the reference genome. Read quantification is performed by averaging the raw counts of mapped reads using the featureCounts program (Liao et al., 2014) which is implemented by the SourceForge Subread package (Liao et al., 2013). The featureCounts program is specifically designed for read counting, which is essential for quantifying reads mapped to specific genomic features defined in the reference file. It operates with two inputs: sorted BAM files containing mapped reads and a GTF reference file (such as hg38.gtf from ENSEMBL) for annotation. The output consists of a tab-delimited table of read counts, where rows correspond to transcripts (genes), and columns represent individual samples. This read counts file is used in subsequent differential expression analysis to identify differentially expressed genes (DEGs) between conditions or groups.

2.5 Differential expression analysis based on skin cancer samples

In this study, we conducted differential expression analysis to identify genes with significant alterations between normal and skin cancer samples. Raw read counts were obtained using featureCounts and served as input for DESeq2, a Bioconductor package in R Studio that applies negative binomial generalized linear models to detect differentially expressed genes (DEGs) (Anders and Huber, 2010). The count data was normalized, and genes were filtered using the false discovery rate (FDR) method, with an adjusted p-value (padj) threshold of < 0.05. We used Wald tests within DESeq2 to identify significant DEGs, applying a p-value cutoff of < 0.05. Genes were categorized as upregulated or downregulated based on their log2-fold change values, with positive values indicating upregulation and negative values indicating downregulation in skin cancer samples (Love et al., 2014). This approach allowed us to uncover key molecular alterations associated with skin cancer, providing valuable insights into potential biomarkers and therapeutic targets.

2.6 Gene Set Enrichment Analysis

Understanding the biological significance underlying differential gene expression profiles can be achieved through pathway analysis using GSEA (Reimand et al., 2019). Gene Set Enrichment Analysis is an effective method widely used to study functional enrichment between two biological groups. In this study, GSEA software was utilized to predict different Gene Ontology (GO) terms present in the gene expression data. The oncogenic signature gene sets (c6.all.v7.5.1.symbols.gmt) were used, considering a significant false discovery rate (FDR) q-value threshold of < 0.02. The analysis involved providing inputs such as the expression dataset file, phenotype labels file, gene sets file in gct, cls, and gmt formats, as well as loading the Chip annotation file. Subsequently, the data files were loaded, analysis parameters were set, and the analysis was executed.

2.7 Enrichment Map and network analysis

Enrichment Map is a freely available and open-source Cytoscape plugin used for visualizing enriched pathways in a network format (Magrane and Consortium, 2011; Merico et al., 2011). It simplifies the interpretation of results obtained from Gene Set Enrichment Analysis by representing relevant pathways as nodes connected by edges that share many genes. The nodes in the network are colored based on the ES score, and the size of the edges is determined by the number of genes shared between connected pathways. To explore individual gene interactions within different pathways, we selectively identify specific genes and construct a network by annotating these genes using the String Database within Cytoscape software (Szklarczyk et al., 2021).

3.1 Quality control and alignment Results

FastQC tool was utilized for quality control analysis on each individual sample. (de Sena Brandine and Smith, 2019) The output HTML and ZIP files from FastQC were scanned, and a comprehensive multiqc report (Ewels et al., 2016) was created to present the overall quality results for all the samples. This multiqc report allowed for easy visualization and comparative analysis of the quality metrics across all samples. Notably, the report included information on duplicate reads, which is presented in Table 1. Furthermore, the reads from the samples were aligned to a reference genome, and the overall alignment rate was calculated (Chen et al., 2022). The results of this alignment, including the alignment rates for each sample, are provided in Table 2, Table 3 and Table 4.

3.2 Functional analysis of DEGs

The differential expression analysis revealed significant alterations in gene expression between skin cancer samples and normal tissues. Two genes, C19orf42 and OST4, were found to be upregulated in skin cancer samples, with log2 fold-changes of 0.28 and 0.39, respectively (Liu et al., 2014b). While the roles of these genes in cancer are not well characterized, their upregulation suggests potential involvement in cancer cell survival, proliferation, and tumor progression processes (Hoang et al., 2017; Liu et al., 2022). C19orf42 has been implicated in some studies as a contributor to tumor cell adaptability, while OST4 plays a role in protein glycosylation, a process often dysregulated in cancer, which may influence cell adhesion, migration, and immune evasion in the tumor microenvironment (Liu et al., 2014a).

In contrast, ten key genes were significantly downregulated, including C3, LCN2, and NFKB2. C3, which is crucial for immune response, had the highest downregulation (log2 fold-change = -1.85) , suggesting immune evasion as a possible mechanism employed by cancer cells. C3 is an essential component of the complement system, which plays a critical role in recognizing and eliminating pathogens and abnormal cells. Its suppression may allow tumor cells to avoid detection and destruction by the immune system, thereby promoting cancer progression (Galvez et al., 2020a).

LCN2, involved in cancer cell migration, invasion, and metastasis, also exhibited substantial downregulation. As LCN2 has been shown to affect cellular iron transport and modulate the tumor microenvironment, its downregulation may indicate alterations in the metastatic potential of skin cancer cells (Galvez et al., 2020b). Similarly, NFKB2, a regulator of immune and inflammatory responses, was downregulated, suggesting impairments in the NF-κB signaling pathway, which is critical for orchestrating immune responses against tumors. The suppression of NFKB2 could hinder the immune system's ability to mount an effective anti-tumor response, potentially allowing the cancer to evade immune surveillance and proliferate unchecked (Love et al., 2014).

Furthermore, other downregulated genes like ATF3, a stress-response gene, and DUSP4, involved in cellular signaling pathways, suggest that the suppression of cellular defense mechanisms may contribute to the aggressiveness of the tumor. ATF3 is often upregulated in response to cellular stress, and its downregulation in skin cancer samples might reflect the cancer cells' ability to tolerate stress, promoting survival in a hostile tumor microenvironment (Anders and Huber, 2010).

While well-established melanoma biomarkers such as HMB-45, tyrosinase, and S100 are invaluable for diagnosing and prognosing melanocytic tumors, their utility is limited in non-melanoma skin cancers (NMSCs). HMB-45 and tyrosinase are specific to melanocytic cells, whereas S100, though highly sensitive, lacks specificity. Incorporating machine learning (ML) techniques, such as Support Vector Machines (SVM) or Random Forest, could enhance the identification and validation of these DEGs as biomarkers. ML models can analyze large-scale gene expression data to predict the most relevant genes involved in tumor progression, immune evasion, and metastasis. Furthermore, ML can refine diagnostic accuracy by classifying samples based on gene expression profiles, allowing for personalized therapeutic strategies (Garg et al., 2022a; Garg et al., 2022b). Targeting downregulated genes like NFKB2 and C3 to restore immune function, or inhibiting upregulated genes such as C19orf42 and OST4 to curb tumor progression, could lead to the development of novel treatments aimed at controlling skin cancer growth and metastasis (Cives et al., 2020b).

3.3 Principal Component analysis (PCA)

Principal Component Analysis (PCA) is a widely used technique that highlights variation and reveals distinct patterns within a dataset. In this study, the DESeq2 package was employed to obtain the dataset. The resulting PCA plot effectively demonstrates the grouping of samples based on two conditions (Shen et al., 2016). By evaluating the expression data, we can assess the relationship between these samples, as the plot represents the counts of one sample versus another. PCA plot is constructed by assessing the normalized differences in expression patterns, which allows us to compute a distance matrix. This matrix captures the dissimilarities between samples based on their gene expression profiles (Bhalla et al., 2019). In the plot, the X-axis and Y-axis represent a mathematical transformation of these distances, effectively visualizing the data in two dimensions—specifically, PC1 (Principal Component 1) versus PC2 (Principal Component 2). Notably, PC1 accounts for 2 % of the variance, while PC2 captures a substantial 96 % of the variance, as illustrated in Fig. 1.

Close

Fig. 1

The PCA plot illustrates the clustering of 12 samples, with six samples categorized as Normal and six samples categorized as Epidermoid Carcinoma. The samples are distributed in a two-dimensional (2D) plane. In the plot, the color scheme allows us to distinguish between the groups, with light blue indicating the Normal samples and red indicating the Epidermoid Carcinoma samples.

Export to PPT

3.4 Heatmap

The heatmap (Fig. 2) visually illustrates the differential expression of key genes between epidermoid carcinoma and normal tissue samples. Hierarchical clustering effectively segregated the samples into distinct groups based on their expression profiles, with epidermoid carcinoma samples clustering separately from normal tissues. Several genes, including C1orf42, KRT23, LCN2, and PI3, exhibited markedly higher expression in carcinoma samples, as indicated by the darker shades, while genes such as PIGR and MYO1E showed lower expression levels in carcinoma compared to normal samples. The heatmap’s color gradient, ranging from yellow (low expression) to dark brown (high expression), clearly distinguishes upregulated and downregulated genes across the sample groups. This visualization underscores the significant variation in gene expression between carcinoma and normal tissues, highlighting potential genes involved in epidermoid carcinoma pathogenesis. The clustering of genes such as TNFAIP3, CFB, and LCN2 further suggests their involvement in shared pathways or processes, which may contribute to the disease's progression.

Close

Fig. 2

Heatmap of Differentially Expressed Genes in Epidermoid Carcinoma and Normal Samples. The heatmap displays the expression levels of 19 differentially expressed genes across epidermoid carcinoma (red) and normal (blue) tissue samples. Color gradients range from yellow (low expression) to dark brown (high expression), with hierarchical clustering showing distinct gene expression patterns between the two groups. Notably, genes such as C1orf42, KRT23, and LCN2 are upregulated in carcinoma samples, while PIGR and MYO1E are downregulated.

Export to PPT

3.5 GSEA and Enrichment Map

GSEA (Gene Set Enrichment Analysis) is a powerful computational method employed to identify specific Gene Ontology (GO) pathways that are significantly associated with a particular disease or biological condition. By analysing high-throughput gene expression data, GSEA helps uncover underlying molecular mechanisms and provides valuable insights into the functional roles of genes within these pathways. The genes linked to the disease are listed in order of rank, with the enrichment score (ES) provided. This analysis helps to pinpoint specific genes that play a significant role in the disease progression, including metastasis, and can aid in early-stage identification. An enrichment map is a visualization tool organizing related gene sets into a network. Gene sets are nodes, while member genes are edges. The map arranges nodes based on similarities among gene sets (Ding et al., 2015). Node color represents gene expression levels in pathways: red for upregulated pathways and blue for downregulated pathways. (as shown in Fig. 3).

Close

Fig. 3

Enrichment Map constructed for skin Cancer using Enrichment Map plugin.

Export to PPT

3.6 Gene network analysis for enriched genes using Cytoscape

The StringApp is a Cytoscape plugin used for visualizing the interconnections between genes. It retrieves functionally enriched genes involved in different pathways, which are presented in Table 5. After importing the genes into this plugin, it initiates a search within the STRING database and proceeds to annotate the targeted gene network (Newman et al., 2019). The nodes in the network represent genes, and the edges indicate similar functional activities between the genes (Kashyap et al., 2022). A String network was constructed using 40 common genes that were predominantly involved in 5–8 pathways (Fig. 4). IL6 plays a pivotal role in the IL6_JAK_STAT3 signaling pathway, which is crucial in promoting tumor growth, inflammation, and immune evasion. The persistent activation of this pathway contributes to uncontrolled cell proliferation and survival in cancer cells, and its involvement in skin cancer suggests a mechanism through which inflammation-driven tumor progression occurs. CCND2 is a key regulator of the cell cycle, and its upregulation promotes cell cycle progression and uncontrolled proliferation. CCND2 activation in skin cancer may facilitate the rapid growth of tumor cells, contributing to tumor aggressiveness and progression. PLAUR is involved in the regulation of extracellular matrix degradation and cell migration through the urokinase plasminogen activator system. Its role in angiogenesis and metastasis is particularly relevant to skin cancer, as it promotes tumor cell invasion and metastasis, further contributing to the cancer’s spread and progression. CD44 is a cell surface glycoprotein involved in cell adhesion, migration, and interactions with the extracellular matrix. Its role in the APICAL_SURFACE pathway indicates its contribution to maintaining the structure and invasiveness of cancer cells, particularly in the context of skin cancer metastasis.

Close

Fig. 4

A string network highlighting the interaction of the most common genes, IL6, CCND2, and PLAUR, was constructed using Cytoscape. The network consisted of 38 nodes and 184 edges and was involved in multiple pathways.

Export to PPT