Skip to main content

ORIGINAL RESEARCH article

Front. Immunol., 18 June 2021
Sec. Autoimmune and Autoinflammatory Disorders

A Computational Systems Analyses to Identify Biomarkers and Mechanistic Link in Psoriasis and Cutaneous Squamous Cell Carcinoma

  • 1Research Center for Modeling and Simulation, National University of Sciences and Technology, Islamabad, Pakistan
  • 2Jackson Laboratory for Genomic Medicine, Farmington, CT, United States
  • 3National Institute of Health (Pakistan), Islamabad, Pakistan

Psoriasis is the most common and chronic skin disease that affects individuals from every age group. The rate of psoriasis is increasing over the time in both developed and developing countries. Studies have revealed the possibility of association of psoriasis with skin cancers, particularly non-melanoma skin cancers (NMSC), which, include basal cell carcinoma and cutaneous squamous cell carcinoma (cSCC). There is a need to analyze the disease at molecular level to propose potential biomarkers and therapeutic targets in comparison to cSCC. Therefore, the second analyzed disease of this study is cSCC. It is the second most common prevalent skin cancer all over the world with the potential to metastasize and recur. There is an urge to validate the proposed biomarkers and discover new potential biomarkers as well. In order to achieve the goals and objectives of the study, microarray and RNA-sequencing data analyses were performed followed by network analysis. Afterwards, quantitative systems biology was implemented to analyze the results at a holistic level. The aim was to predict the molecular patterns that can lead psoriasis to cancer. The current study proposed potential biomarkers and therapeutic targets for psoriasis and cSCC. IL-17 signaling pathway is also identified as significant pathway in both diseases. Moreover, the current study proposed that autoimmune pathology, neutrophil recruitment, and immunity to extracellular pathogens are sensitive towards MAPKs (MAPK13 and MAPK14) and genes for AP-1 (FOSL1 and FOS). Therefore, these genes should be further studied in gene knock down based studies as they may play significant role in leading psoriasis towards cancer.

Introduction

Psoriasis is one of the most common and chronic skin diseases (1). It is characterized by rough scales on the skin. It could appear on any part of the body and it affects both children and adults (1). About 2-3% individuals are affected by psoriasis worldwide, that means more than 125 million people suffer from psoriasis (2). Males and females are equally targeted by psoriasis i.e., male to female ratio is 1:1. The usual onset of the disease is reported to be between 18 and 39 years, or between 50 and 69 years (2). The onset age of the disease is associated with genetic and epigenetic factors (2). The incidence rate of psoriasis is also observed to increase with time (3). The pathogenesis of psoriasis is still unclear and there are two hypotheses regarding it (4). The first hypothesis assumes that in psoriasis the immune system of the body appears to be hyperactive to an extent that it attacks itself (4); which is the characteristic of the autoimmune disease. Therefore, psoriasis shows characteristics of an autoimmune disorder (4, 5). According to the second hypothesis, the life cycle of the skin cells speeds up in psoriasis from 28 days to 3-6 days, which results in excessive cell growth (4, 6). It eventually causes the formation of dead skin cells. Despite of being shed off, these dead cells accumulate at the surface of the skin and form lesions and patches (6). Frequently, red patches or dry scales appear on elbows, knees, and scalps (7). In an individual’s lifetime psoriasis often occur repeatedly (7). The exact cause and cure of psoriasis has not yet been discovered. However, multiple factors including environmental & genetic factors and immune system of the body are considered for it (8). It is predicted that environmental factors like cold and dry weather could trigger the symptoms of psoriasis (9). Moreover, the symptoms of psoriasis could appear after certain infections like strep throat and even in stress (9). Both genetic and epigenetic factors play the substantial role in the formation and exacerbation of the symptoms of the disease (10). Genetic risk factors are HLA-Cw6, mutations in CARD14 (10). CARD14 is caspase recruitment domain family member 14 gene. Obesity, stress, smoking, alcohol consumption, HIV infection, and certain medications serve as the epigenetic factors of the disease (11).

Association of psoriasis with other comorbidities, like cardiovascular diseases, diabetes, and metabolic syndrome, has been observed (12). An interesting study (13) was conducted for analyzing the implication of cutaneous squamous cell carcinoma (cSCC) antigen with psoriasis. This study was based on the fact that serum (fluid after blood coagulation) examination of psoriasis patients always revealed the presence of cSCC antigen. The study confirmed the concentration and over-expression of cSCC antigen in psoriatic lesions. A transcriptome level analysis was performed to reveal common patterns in cSCC and psoriasis (14). In this experiment, expression levels in cSCC were compared with multiple cutaneous disorders by performing microarray analysis. The comparison of cSCC with psoriasis (psoriasis vulgaris) revealed some genes depicting same expression behavior by up regulation in both disease states of cSCC and psoriasis. Subsequently, RT (real time) reverse transcription PCR (polymerase chain reaction) along with immunohistochemistry were also performed to validate the differential expression of genes attained after comparison in both disease states of cSCC and psoriasis. The common over-expressed genes found in both cSCC and psoriasis were: DEFB4 (defensin B4), SERPINB3, STAT1, K16 (keratin 16), WNT5A, and CEACAM5. In another study (15), histopathological examination of a biopsy of a psoriasis patient revealed the development of cSCC. That patient showed negative family history regarding psoriasis and only received homeopathic treatment for psoriasis.

Studies have revealed the possibility of association of psoriasis with skin cancers, particularly non melanoma skin cancer (NMSC), which include basal cell carcinoma (BCC) and squamous cell carcinoma (SCC) (12, 16). A cohort study (1988-2011) (16) was conducted for estimating malignancy rates in psoriasis patients treated with systemic therapy. This study reported that there exists a higher rate of developing cSCC among psoriasis patients who were treated with biologics as a treatment. The speculation regarding the association of psoriasis with NMSC was built on the fact that immunosuppressors and ultraviolet phototherapy are used for treating psoriasis (12, 16). Moreover, the long-lasting inflammatory nature of psoriasis might also escalate the threat of NMSC in psoriasis patients (12). That is why, the second studied disease in this study is cSCC.

Most commonly, skin cancer is characterized by the abnormal growths of the skin (17). A sore that does not heal or takes time to heal could also be the indicator of skin cancer. Symptoms of skin cancer vary differently for its different subtypes. For instance, it could appear as small, smooth lumps. However, in some cases skin cancer appears as scaly, hard lumps that often bleed (17). Skin cancer is divided into melanoma and non-melanoma categories. Basal Cell Carcinoma (BCC) and Squamous Cell Carcinoma (SCC) set up non-melanoma skin cancer (17, 18). Squamous Cell Carcinoma is also known as cutaneous squamous cell carcinoma (cSCC). cSCC is the second most prevalent non-melanoma skin cancer (19, 20) and has the potential to recur/revert (19). It is mostly reported to appear in the Caucasian population (21). Furthermore, cSCC also accounts for about 20% of all the cutaneous cancers (19, 20, 22, 23). Regardless of the low occurrence of cutaneous SCC, it is responsible for the majority of the metastasis in disease, eventually causing deaths (24). It is stated to be counted in the list of topmost expensive cancers in the US (United States) (24). The prevalence rate of the cSCC is also vastly growing annually (25). Cutaneous SCC is caused by malignant multiplication and spread of epidermal keratinocytes (26). Often it appears on the sun exposed areas of the body; like face, head, hands, dorsal side of forearms and neck regions (17), and is more prevalent in elderly males (26). It has also been observed and reported that the particular rough patches resulted from sun exposure (aka Actinic Keratosis) also transform or lead to cSCC (17, 26, 27). Light skin tone, ultraviolet radiations, arsenic compounds, and immunosuppression are the reported risk factors and are found to increase the likelihood of developing cSCC (28). Some of the common pathological and clinical factors for the examination of the disease are age, gender, ethnicity, tumor size, histology, organ transplantation, psoriasis, leukemia, history of previously injured skin & nodal metastasis (29). Along with patient’s characteristics, tumor site & size, tumor depth, and histologic features accounts for recurrence and metastasis risks of cSCC (26). Genes which are reported to be frequently mutated in cSCC are TP53, CDKN2A, Ras, and NOTCH1 (20).

With the increasing annual incidence of psoriasis and cSCC, there is a great urge to comprehend and compare the diseases at the molecular level by interpreting the genomic data. There is a need to validate known biomarkers of psoriasis and cSCC through high throughput data. Moreover, common therapeutic target for both the diseases is not yet discovered. It has also been discussed that prolonged psoriasis could lead to cSCC (15). According to our understanding, the association of both the diseases has not been comparatively analyzed at molecular level in detail. In the current study, microarray and RNA-seq data analyses have been performed on the publicly available datasets of both psoriasis and cSCC to reveal the common and uncommon patterns in both the diseases.

Materials and Methods

The workflow to identify potential biomarkers, therapeutic targets, and to comprehend the association between cSCC and psoriasis, is given in Figure 1 and described below in detail.

FIGURE 1
www.frontiersin.org

Figure 1 General Workflow of the Project. In order to achieve the aims and objectives of this study, integrated data analysis of microarray and RNA-seq were performed for cSCC and psoriasis. From differential expression analysis, common DEGs between psoriasis and cSCC were identified that are proposed as potential biomarkers and therapeutic targets for both the diseases. Pathway analysis was then performed to analyze the disrupted pathways and their possible relation with the disease. The significant pathways (at p-value <= 0.05) from every dataset were then compared for common pathway. As a result of this comparison, one common pathway, named IL-17 signaling pathway, was identified. Afterwards, quantitative systems biology was implemented to analyze the results at a holistic level. A sub-section of complete IL-17 pathway was modeled by using SimBiology application in Matlab (30). The average expressions values of DEGs that were present in the pathway were calculated for microarray datasets of psoriasis and cSCC. The model was then simulated by adding the average expression values. After simulation of the model, sensitivity analysis was performed to analyze the response of the biological processes to model quantities (DEGs and genes). Sensitivity analysis of the biological processes of both psoriasis and cSCC revealed significant genes that might lead psoriasis to cancer and could also be proposed as potential biomarkers and therapeutic targets of both the diseases.

Datasets Inclusion Criteria

Publicly available datasets of human origin, free from drugs, cell lines, and mutations, were selected from EMBL-EBI (https://www.ebi.ac.uk/) and GEO (https://www.ncbi.nlm.nih.gov/geo/) databases against query words like, psoriasis and cutaneous squamous cell carcinoma (cSCC). Additionally, to make the statistical analysis significant, those datasets with enough number of samples were selected. The Table 1 shows the details of analyzed datasets of both the diseases.

TABLE 1
www.frontiersin.org

Table 1 Microarray and RNA-seq datasets.

A total of 3 datasets were analyzed for microarray data analysis of cSCC, as given in Table 1. Two of these datasets, with accession numbers E-GEOD-32868 and E-GEOD-34536, are of miRNA type, containing 4 and 14 samples, respectively. The third dataset, with accession number E-GEOD-42677 is of mRNA type, containing 20 samples.

For microarray data analysis of psoriasis, a total of 5 datasets were analyzed, as given in Table 1. All these datasets are of mRNA type, with accession numbers E-GEOD-75343, E-GEOD-50790, E-GEOD-41664, E-GEOD-34248 and E-GEOD-13355, containing 45, 8, 76, 28 and 180 samples, respectively.

Datasets for RNA-seq analysis of both the diseases are given in Table 1. All datasets are of mRNA type and from illumina platform. Two datasets were analyzed for cSCC, with accession numbers E-MTAB-5678 and GSE84293, containing 9 and 15 samples, respectively. Two datasets were analyzed for psoriasis, with accession numbers E-GEOD-41745 and E-GEOD-67785, containing 6 and 28 samples, respectively.

Differential Expression Analysis

Microarray Data Analysis

Microarray data analysis was performed in R (version of 3.5.1). A published maEndToEnd pipeline (31) for Affymetrix platform, was used for microarray datasets of mRNA type. The quality control (QC) of the data was performed by data visualization and filtering techniques which include box plot, principal component analysis (PCA) plot, relative log expression (RLE) plot, and Heatmap clustering analysis. PCA and heatmap clustering analysis were performed and plotted to visualize the logarithmic transformations by using Bioconductor R packages (32, 33). Afterwards, the probes having low intensity along with the low variance were removed to avoid noise in data. Then to normalize gene expression intensities quantile normalization method (34, 35) was used. The differentially expressed genes (DEGs) were identified by fitting linear model in disease samples relative to control/normal (defined groups) by using the R package limma (36). The R package limma has linear model features and is specifically designed for comparative analysis between various experimental groups. Empirical bayes variance moderation method was implemented to the data by using eBayes function to generate the statistical significance values for DEGs (32).

The two selected microarray datasets of cSCC (E-GEOD-32868 and E-GEOD-34536) are of miRNA type. The dataset with accession number E-GEOD-32868 is of Affymetrix platform. The analysis of that dataset was performed by modifying the maendtoend pipeline (31) by omitting the ‘Annotation step’, as annotation is not required for microarray data analysis of miRNAs. The dataset with accession number E-GEOD-34536 is from Agilent platform. The analysis of that dataset was performed by AgiMicroRna package (37) in R. The AgiMicroRna package integrates the RMA (robust multiarray average) algorithm for preprocessing and differential expression computation of Agilent miRNA data. The RMA algorithm implements the linear model from the R package limma (36). After obtaining miRNAs, their target genes were extracted from miRWalk (38). miRWalk is a publicly available database of miRNAs targets, which integrates multiple prediction programs in a user-friendly interface. For extraction of target genes, p-value was set to < 0.005 and all the prediction programs were selected.

RNA-Sequencing Data Analysis

A published ‘new Tuxedo’ protocol (39) for RNA-seq data analysis was followed in this study. A Unix based workstation was used for this analysis. To address the issues regarding raw data reads of RNA-seq data, like adapter contamination, over-represented or inaccurately represented sequences and base correction etc., fastp (40) preprocessor was used. After preprocessing, the pre-processed reads were aligned to the genome by using Hisat2 (41) aligner. Afterwards, StringTie (42) was used to assemble complex aligned reads. Transcripts were assembled individually for each sample of the dataset being analyzed. After successful assembling, the collective merging of transcripts from all samples was also performed by using StringTie. In order to match the output from StringTie with reference file, a specialized program named gffcompare (http://ccb.jhu.edu/software/stringtie/gff.shtml) was used. The aligned and assembled reads were then subjected to quantification by using StringTie. Quantification refers to the computation of abundances of mapped assembled transcripts. Finally, Ballgown package (43) of R was used for computing differential expression.

Parameters for Identifying DEGs

In this study, a gene is identified as differentially expressed when it crossed both the defined thresholds of p-value and fold change. Stringent thresholds of p-value (p-value < 0.005) and fold change (-1 > Log2FC | Log2FC > +1) were used between groups of normal and disease phenotypes. Less stringent thresholds of p-value (p-value < 0.05) and fold change (-0.5 > Log2FC | Log2FC > +0.5) were used when the comparisons were made between subtypes of a disease, and between non-lesional and normal phenotypes. Our reason for this is that the non-lesional skin of the psoriasis patients behave normally, and the differential expression between subtypes of a disease would be smaller as compared to the healthy counterparts. The DEGs were visually analyzed by using the R package EnhancedVolcano (44).

Pathway Analysis

Pathway analysis was performed to analyze the disrupted pathways and their possible relation with both diseases, based on the number of DEGs. The enriched pathways of DEGs were identified using Enrichr (45) database and KEGG based pathways were identified. The significant pathways (at p-value <= 0.05) from every dataset were then compared for common pathway. As a result of this comparison, one common pathway, named IL-17 signaling pathway, was identified (given in Table 4).

Quantitative Systems Biology

IL-17 signaling pathway was selected for further analysis, as it was common in both psoriasis and cSCC datasets (The KEGG IL-17 signaling pathway - Homo sapiens is provided in Supplementary File 1). The average expression values of the DEGs present in the IL-17 pathway were calculated for both psoriasis and cSCC microarray datasets, as given in Table 5. Quantitative systems biology was implemented by using ‘SimBiology’ application in Matlab (30). A sub-section of complete IL-17 signaling pathway was modeled in Matlab for analyzing the sensitivities of autoimmune pathology (AP), neutrophil recruitment (NR), and immunity to extracellular pathogens (IEP) towards DEGs, by using block diagram editor in SimBiology. The species were used as biological entities (genes/biological processes) in the model. The activation and inhibition reactions were then added among species in accordance with the IL-17 signaling pathway. The model was then simulated by adding the average expression values to the species (genes) as initial state values. The genes that were not differentially expressed in the pathway were set to 0.5 as initial state value. The reason we have proposed that value to be at 0.5 because it resides between the values of 0 and 1, which are the minimum and maximum values that can be settled for these genes. The biological processes being analyzed in the model were autoimmune pathology, neutrophil recruitment, and immunity to extracellular pathogens. The initial state values of these biological processes were set at zero because we wanted to analyze the response of these biological processes to model quantities (DEGs and other). Once all the initial state values were set, the model was simulated separately for DEGs in psoriasis and cSCC by using ‘SimBiology’ application in Matlab (30). As a result of simulations, graphs were generated (given in Figure 5).

After simulating the model, sensitivity analysis (46) was performed. We wanted to analyze the assumptions about the influences of model entities and reactions on overall model response. Sensitivity of biological processes (autoimmune pathology, neutrophil recruitment, and immunity to extracellular pathogens) in a model were analyzed against genes. All the genes were taken as input genes, whereas biological processes were taken as output. The sensitivity graphs were then generated for further analysis (given in Figures 6 and 7).

Results

Differential Expression Analysis

Microarray Data Analysis of Psoriasis

The microarray data analysis was performed on the selected psoriasis datasets using the methodology described comprehensively in the section “Methods”. In the current study, five mRNA type microarray data sets were selected for psoriasis as given in Table 1. For two datasets of psoriasis (E-GEOD-13355 & E-GEOD-75343), three phenotypic subgroups were available, which were used to draw comparison between three groups including lesional & non-lesional, lesional & normal, and non-lesional & normal.

The common genes among all the subgroups based comparisons were identified. Supplementary Figures S1–S3 and S4–S6 depict the volcano plots displaying DEGs between the defined subgroups of datasets E-GEOD-75343 and E-GEOD-13355, respectively. The red dots in the volcano plots are the identified DEGs. Supplementary Figures S7–S9 depict the volcano plots displaying DEGs of E-GEOD-50790, E-GEOD-41664 and E-GEOD-34248 datasets, respectively. The number of identified DEGs from all the analyzed microarray datasets of psoriasis are given in Table 2. The datasets of psoriasis were compared for common genes by comparing the identified DEGs from all five analyzed datasets. As a result of this comparison, one common DEG, named s100a8, was identified. S100a8 is a protein coding gene that mainly functions in inflammatory and immune responses (47).

TABLE 2
www.frontiersin.org

Table 2 Number of identified DEGs/differentially expressed miRNAs of microarray datasets.

RNA-Seq Data Analysis of Psoriasis

Differentially expressed genes were identified between lesional & non-lesional phenotypes in RNA-seq data analysis of psoriasis. Two datasets of psoriasis were analyzed for RNA-seq analysis as shown in Table 1. The number of identified DEGs from these two datasets is given in Table 3. From datasets E-GEOD-41745 and E-GEOD-67785, 153 and 833 DEGs were identified, respectively. Supplementary Figures S10 and S11 depict the volcano plots displaying DEGs between psoriatic lesional & non-lesional phenotypes of datasets E-GEOD-41745 and E-GEOD-67785, respectively. The datasets were then compared for common genes by comparing the identified DEGs from both the analyzed datasets. As a result of this comparison, 16 common DEGs were identified. The common DEGs are given in Table 3. This comparison also revealed that datasets E-GEOD-41745 and E-GEOD-67785 have 137 and 817 unique identified DEGs, as shown in Figure 2.

TABLE 3
www.frontiersin.org

Table 3 Number of identified and symbols of common DEGs of RNA-seq datasets.

FIGURE 2
www.frontiersin.org

Figure 2 (A) Comparative analysis of RNA-seq datasets of psoriasis. The analysis revealed 16 common DEGs between psoriasis datasets, and, 137 and 817 unique identified DEGs from datasets E-GEOD-41745 and E-GEOD-67785, respectively. (B) Comparative analysis of RNA-seq datasets of cSCC. The analysis revealed 26 common DEGs between cSCC datasets, and, 1667 and 323 unique identified DEGs from datasets E-MTAB-5678 and GSE84293, respectively.

Microarray Data Analysis of cSCC

Studies have proposed the association of psoriasis with NMSC, which include BCC and cSCC (12, 16). The chronic inflammatory nature of psoriasis (12) and use of biologics (16) might also increase the risk of cSCC in psoriasis patients. Therefore, cSCC is analyzed in the current study to interpret its association with psoriasis at the molecular level. A summary of datasets analyzed for cSCC is provided in Table 1. Two miRNA (E-GEOD-32868 and E-GEOD-34536) and one mRNA datasets (E-GEOD-42677) of cSCC are analyzed in this study. A total of six and 12 differentially expressed miRNAs were identified from E-GEOD-32868 and E-GEOD-34536, respectively (provided in Supplementary Tables S12 and S13). The number of target genes of these differentially expressed miRNAs is given in Table 2. The number of identified DEGs from E-GEOD-42677 dataset is given in Table 2. The datasets E-GEOD-32868 and E-GEOD-34536 are composed of two phenotypes i.e., normal and cSCC. However, the dataset E-GEOD-42677 is composed of three phenotypes i.e., normal, cSCC-invasive, and cSCC-in situ. Therefore, the comparisons for finding differential expression between all the possible combinations of phenotypes within the dataset were performed for this dataset. The comparisons were between: cSCC-invasive & normal, cSCC-insitu & normal, and cSCC-invasive & cSCC-insitu groups. Supplementary Figures S14–S16 show the volcano plots of DEGs identified in the defined comparisons of dataset E-GEOD-42677. Supplementary Figures S12 and S13 depict volcano plots displaying differentially expressed miRNAs between cSCC and normal phenotypes, for datasets E-GEOD-32868 and E-GEOD-34536, respectively. The datasets were compared for common genes by comparing the identified target genes of differentially expressed miRNAs with that of mRNAs. For two miRNA datasets (E-GEOD-32868 & E-GEOD-34536) 200 common DEGs were identified. However, no common DEGs were identified between two miRNA datasets and one mRNA dataset of cSCC.

RNA-Seq Data Analysis of cSCC

DEGs were identified between normal & cSCC phenotypes in RNA-seq data analysis of cSCC. Two datasets of cSCC were analyzed for RNA-seq analysis as shown in Table 1. The Table 3 shows the number of identified DEGs from RNA-seq datasets of cSCC. From datasets E-MTAB-5678 and GSE84293, 1693 and 349 DEGs were identified, respectively. Supplementary Figures S17 and S18 depict the volcano plots displaying DEGs between cSCC and normal phenotypes of E-MTAB-5678 and GSE-84293, respectively. The datasets were then compared for common genes by comparing the identified DEGs from both the analyzed datasets. As a result of this comparison, 26 common DEGs were obtained. The common DEGs are given in Table 3. This comparison also revealed that E-MTAB-5678 has 1667 and GSE84293 has 323 unique identified DEGs, as shown in Figure 2.

Pathway Analysis

Pathway analysis of DEGs of psoriasis and cSCC datasets revealed the significant enrichment of these DEGs in various pathways. The significant pathways (at p-value <= 0.05) of all the analyzed datasets are provided in Supplementary Table 1. Out of these significant pathways, IL-17 signaling pathway was found common among the analyzed datasets of both psoriasis and cSCC, as given in Table 4.

TABLE 4
www.frontiersin.org

Table 4 Common pathway of both the diseases. The DEGs of both the datasets of psoriasis and cSCC were enriched in IL-17 signaling pathway.

Quantitative Systems Biology

Simulation of the Model

The DEGs of psoriasis and cSCC were significantly enriched in the IL-17 signaling pathway. A sub-section of complete IL-17 signaling pathway was modeled in Matlab for analyzing the sensitivity of autoimmune pathology, neutrophil recruitment, and immunity to extracellular pathogens towards DEGs. These models for psoriasis and cSCC are shown in Figures 3 and 4, respectively. The models were simulated by inserting the average expression values (given in Table 5) of DEGs in psoriasis and cSCC enriched in IL-17 signaling pathway. The simulation of model unveiled the alterations in the concentrations of species over time. The Figures 5A, B depict simulated processes of autoimmune pathology, neutrophil recruitment, and immunity to extracellular pathogens in psoriasis and cSCC, respectively. It is cleared from Figures 5A, B that the rate of these biological processes increases over time in both psoriasis and cSCC. However, the rate is relatively higher in psoriasis. The peak value of these biological processes in psoriasis is around 28 in 9.5 time units, whereas it is observed to be around 21 in 10 time units in cSCC.

FIGURE 3
www.frontiersin.org

Figure 3 Model of IL-17 signaling pathway for psoriasis. The model is built as network of species (oval) and reactions (circle), which are connected by edges. Species represent biological entities of the pathway, and are color coded. Dark blue color shows processes of autoimmune pathology, neutrophil recruitment, and immunity to extracellular pathogens, yellow color depicts DEGs, green color shows DEG that are also acting as inhibitor, red color shows inhibitors, and sky-blue color represents genes that are part of the pathway, however, not differentially expressed. The species MAPKs includes MAPK13 and MAPK14. Chemokines includes CXCL1, CXCL2, CXCL5, CXCL8, CXCL10, CCL2, CCL7 and CCL20. The species Anti-microbial includes S100A7, S100A8, S100A9, and LCN2 genes. Tissue remodeling includes MMP1 and MMP9 genes. C1, C2, C3 and AP-1 are complexes of genes. C1 includes Hsp90 and Act1. C2 includes TAB2, TAB3 and TAK1 genes. C3 incorporates ANAPC5 and A20 genes, and AP-1 includes FOSL1 and FOS genes.

FIGURE 4
www.frontiersin.org

Figure 4 Model of IL-17 signaling pathway for cSCC. The model is built as network of species (oval) and reactions (circle), which are connected by edges. Species represent biological entities of the pathway, and are color coded. Dark blue color shows processes of autoimmune pathology, neutrophil recruitment, and immunity to extracellular pathogens, yellow color depicts DEGs, red color shows inhibitors, and sky-blue color represents genes that are part of the pathway, however, not differentially expressed. The species MAPKs includes MAPK3, MAPK9, MAPK10, MAPK11, and MAPK14. Chemokines includes CXCL1, CXCL2, CXCL8, and CCL20. The species Anti-microbial includes DEFB4A, DEFB4B, S100A7, S100A8, S100A9, and LCN2 genes. Tissue remodeling includes MMP1, MMP3, MMP9, and MMP13 genes. C1, C2, C3 and AP-1 are complexes of genes. C1 includes Hsp90 and Act1. C2 includes TAB2, TAB3 and TAK1 genes. AP-1 includes FOSL1 and JUN genes.

TABLE 5
www.frontiersin.org

Table 5 DEGs of psoriasis and cSCC significantly enriched in IL-17 signaling pathway.

FIGURE 5
www.frontiersin.org

Figure 5 Simulated graphs of autoimmune pathology, neutrophil recruitment, and immunity to extracellular pathogens. The x-axis represents time, whereas y-axis depicts concentration. (A) Simulated graph of autoimmune pathology, neutrophil recruitment, and immunity to extracellular pathogens in psoriasis. In psoriasis, activation of autoimmune pathology, neutrophil recruitment, and immunity to extracellular pathogens start immediately and reach their peak value around 28 in 9.5 time units and become constant after it. (B) Simulated graph of autoimmune pathology, neutrophil recruitment, and immunity to extracellular pathogens in cSCC. Autoimmune pathology, neutrophil recruitment, and immunity to extracellular pathogens activate immediately in cSCC and reach their peak value around 21 in 10 time units.

Sensitivity Analysis

Sensitivity analysis was performed to study the impact of the species and reactions of the model on biological processes of autoimmune pathology, neutrophil recruitment, and immunity to extracellular pathogens. The Figures 6 and 7 depict the sensitivity analysis of these biological processes in psoriasis and cSCC, respectively. It is cleared from Figure 6 that the biological processes are highly sensitive towards MAPKs, chemokines, tissue remodeling genes, and AP-1 (FOSL1 and FOS genes) in psoriasis, having sensitivity values greater than 1.5. Whereas in cSCC, these biological processes are highly sensitive towards MAPKs, AP-1 (FOSL1 and JUN genes), and TNF genes, having sensitivity values greater than 1.5 as shown in Figure 7.

FIGURE 6
www.frontiersin.org

Figure 6 Sensitivity analysis of Autoimmune pathology (AP), neutrophil recruitment (NR), and immunity to extracellular pathogens (IEP) in psoriasis. The x-axis represents input species (genes), whereas y-axis depicts sensitivity of AP, NR and IEP. The high influence of MAPKs, chemokines, tissue remodeling genes, and AP-1 is observed in AP, NR and IEP. It indicates that AP, NR and IEP are highly sensitive towards these genes, having sensitivity values greater than 1.5.

FIGURE 7
www.frontiersin.org

Figure 7 Sensitivity analysis of Autoimmune pathology (AP), neutrophil recruitment (NR), and immunity to extracellular pathogens (IEP) in cSCC. The x-axis represents input species (genes), whereas y-axis depicts sensitivity of AP, NR and IEP. The high influence of MAPKs, AP-1, and TNF is observed in AP, NR and IEP. It indicates that AP, NR and IEP are highly sensitive towards these genes, having sensitivity values greater than 1.5.

Discussion

Psoriasis is one of the most common and chronic skin diseases (1). It is characterized by rough scales on the skin. It could appear on any part of the body and it affects both children and adults (1). About 2-3% individuals are affected by psoriasis worldwide (2), and the incidence rate of psoriasis is also observed to increase with time (3). Association of psoriasis with other comorbidities, like cardiovascular diseases, diabetes, and metabolic syndrome, has been observed (12). Studies have also revealed the possibility of association of psoriasis with skin cancers, particularly non-melanoma skin cancers (NMSC), which include basal cell carcinoma (BCC) and cutaneous Squamous Cell Carcinoma (cSCC) (12, 16). Hence, the second analyzed disease in this study is cSCC. It is the second most prevalent nonmelanoma skin cancer (19, 20) with potential to metastasize and recur in some cases (19). Moreover, the incidence rate of cSCC is also increasing annually (19, 20). Therefore, there exists a dire need to analyze both the diseases at the molecular level to propose potential biomarkers and common therapeutic targets.

This study was designed to interpret and compare the high-throughput data of psoriasis and cSCC, aiming to reveal potential biomarkers, common therapeutic targets, and the association between the two diseases. Microarray and RNA-seq data analyses were performed on the datasets of psoriasis and cSCC (given in Table 1) to reveal the common and uncommon patterns. After performing differential expression analysis, pathway analysis was performed to analyze the disrupted pathways and their possible relation with the disease. The common pathway between the diseases was selected for further analysis. Afterwards, quantitative systems biology was implemented, by using SimBiology application in Matlab (30), to analyze the results at a holistic level. The differential expression analysis of both the analyzed diseases resulted in the identification of DEGs. Top 10 DEGs (along with p-value and Log2FC) from all the analyzed datasets are provided in Supplementary Tables S1–S18. The identified number of DEGs from microarray data analysis are given in Table 2. The comparison of microarray datasets of psoriasis resulted in one common DEG, named S100A8. The significant role of S100A8 in prognosis of psoriasis is also illustrated in the literature (48). The overexpression of S100A8 causes keratinocyte proliferation and angiogenesis in psoriatic skin (47, 49). The comparison of microarray datasets of cSCC resulted in no common gene. The number of identified DEGs from RNA-seq data analysis is given in Table 3. The comparison of RNA-seq datasets resulted in 16 and 26 common DEGs for psoriasis and cSCC, respectively (given in Table 3). These common DEGs could be proposed as potential biomarkers of the diseases. The overall comparative analysis of identified DEGs from all the analyzed datasets of both the diseases resulted in common DEGs between the two diseases, as given in Table 6. These DEGs may represent the potential biomarkers or therapeutic targets of psoriasis and cSCC.

TABLE 6
www.frontiersin.org

Table 6 Potential biomarkers and therapeutic targets for psoriasis and cSCC.

Pathway analysis of DEGs of psoriasis and cSCC revealed the enrichment of these DEGs in various pathways (provided in Supplementary Table 1). The current study has proposed one common significant pathway, IL-17 signaling pathway, based on the enrichment of DEGs of both the diseases in that pathway (given in Table 4). Genome-wide association and other genetic studies have clearly linked multiple IL-17–related genes to psoriasis pathogenesis. Psoriasis is one of the most typical IL-17A-driven human diseases (50). The expression of various IL-17 signaling pathway genes elevated in psoriatic lesional tissue compared with non-lesional tissue (51). The significance of IL-17 signaling pathway in the prognosis of cSCC is also illustrated in the literature, as the disrupted IL-17 pathway promotes tumor development by activating cytokines production (52).

IL-17 signaling pathway was selected for further analysis, as it was common in both psoriasis and cSCC datasets. IL-17 signaling play a critical role in tumor microenvironment by facilitating tumor formation and metastasis (53). Quantitative systems biology was implemented to analyze the results at a holistic level. A sub-section of complete IL-17 pathway was modeled by using SimBiology application in Matlab (30). The average expressions values of DEGs that were present in the pathway were calculated for microarray datasets of psoriasis and cSCC, as given in Table 5. The model was then simulated by adding the average expression values to the species (genes) as initial state values. The biological processes being analyzed in the model were autoimmune pathology, neutrophil recruitment, and immunity to extracellular pathogens. Autoimmunity is defined as when the immune system of the body becomes hyperactive and acts against self-antigen, thus, resulting in tissue damage (54). Cancer is often associated with autoimmunity. Autoimmune disorders like scleroderma and myositis might increase the risk of cancer (55). Neutrophils are the type of leukocytes that protect the body in state of infections and inflammation (56). Clinical evidence has revealed that tumor associated neutrophils (TANs) might play a critical role in tumor progression. Adaptive immune system of the host becomes activated by exhibiting particular responses to extracellular pathogens (57). The clinical evidence has suggested that abnormal adaptive immunity lead to tumor cell proliferation by inducing immunosuppression (58).

The Figures 5A, B depict simulated processes of autoimmune pathology, neutrophil recruitment, and immunity to extracellular pathogens in psoriasis and cSCC, respectively. It is cleared from Figures 5A, B that the rate of these biological processes increases over time in both psoriasis and cSCC. However, the rate is relatively higher in psoriasis. The peak value of these biological processes in psoriasis is around 28 in 9.5 time units, whereas it is observed to be around 21 in 10 time units in cSCC.

After simulation of the model, sensitivity analysis was performed to analyze the response of the biological processes to model quantities (DEGs and genes). The sensitivity analysis revealed that autoimmune pathology, neutrophil recruitment, and immunity to extracellular pathogens are highly sensitive towards MAPKs, chemokines, tissue remodeling genes, and AP-1 (FOSL1 and FOS genes) in psoriasis, having sensitivity values greater than 1.5 Figure 6. It indicates that the slight changes in these genes profoundly affect the biological processes. These results are also consistent with the literature. Disrupted signaling pathways are associated with the complex pathogenesis of psoriasis, particularly, the possible pathogenic role of disrupted MAPK signaling in the development of psoriasis (59, 60). Chemokines and chemokine receptors are involved in the pathogenesis of psoriasis by facilitating the trafficking of T-cells in psoriasis (61). Therefore, targeting chemokines in psoriasis would result in therapeutic effect and improved treatment for psoriasis (61, 62). Tissue remodeling genes (MMPs) contribute to the pathogenesis of psoriasis by stimulating angiogenesis and intrusion of immune cells in the dermal tissue (63). AP-1 (the transcription factor activator protein 1) functions in cell proliferation. However, the binding ability of AP-1 is highly damaged in psoriatic lesional skin (64), thus genes for AP-1 are significant in pathogenesis of psoriasis and proposed as candidate target genes (65).

Whereas the sensitivity analysis of cSCC revealed that the biological processes are highly sensitive towards MAPKs, AP-1 (FOSL1 and JUN genes), and TNF genes, having sensitivity values greater than 1.5 as shown in Figure 7. The slight changes in these genes profoundly affect the biological processes. These results are also consistent with the literature. Dysregulated MAPK signaling by high levels of HMG-1 protein stimulates metastasis in cSCC (66). Dysregulation in the expression of the genes for AP-1 (the transcription factor activator protein 1) promote cancer development as AP-1 factors are involved in keratinocyte proliferation of skin (67). TNF(Tumor necrosis factor) is a proinflammatory cytokine that mainly functions in cell proliferation, differentiation, and death, hence plays an important role in carcinogenesis (68).

Sensitivity analysis of the biological processes of both psoriasis and cSCC revealed MAPKs and genes for AP-1 as significant genes. Slight change in these genes would greatly impact autoimmune pathology, neutrophil recruitment, and immunity to extracellular pathogens. Therefore, consistent with the discussion we propose that MAPKs (MAPK13 and MAPK14) and genes for AP-1 (FOSL1 and FOS) identified from sensitivity analysis of psoriasis are significant genes that might lead psoriasis towards cancer.

Conclusion

In the current study, the integrated data analysis of microarray and RNA-seq was performed to reveal common patterns between psoriasis and cSCC. The differential expression analysis of both the diseases resulted in common DEGs, which are proposed as potential biomarkers and therapeutic targets for psoriasis and cSCC. IL-17 signaling pathway is identified as common significant pathway as DEGs of both the diseases were enriched in that pathway. IL-17 was further studied, and sensitivity analysis was performed that resulted in the identification of significant genes in psoriasis, that might lead psoriasis towards cancer. These significant genes are MAPKs (MAPK13 and MAPK14) and genes for AP-1 (FOSL1 and FOS). These genes should be further investigated at experimental level for validating their roles in the transfer of psoriatic condition to cSCC.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author Contributions

SA and RP conceived and designed the study. SA conducted the analysis and wrote the manuscript. All authors contributed to the article and approved the submitted version.

Funding

National University of Sciences and Technology, Islamabad paid the publication fee for this article.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We are grateful for the supercomputing research and education center (ScREC) at Research Center for Modelling and Simulation (RCMS) NUST, Pakistan, for providing the facility of supercomputer to carry out this research.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2021.662528/full#supplementary-material

References

1. Kim WB, Jerome D, Yeung J. Diagnosis and Management of Psoriasis. Can Family Physician (2017) 63:278–85. doi: 10.1016/j.jaad.2016.10.006

CrossRef Full Text | Google Scholar

2. Parisi R, Symmons DP, Griffiths CE, Ashcroft DM. Global Epidemiology of Psoriasis: A Systematic Review of Incidence and Prevalence. J Invest Dermatol (2013) 133:377–85. doi: 10.1038/jid.2012.339

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Icen M, Crowson CS, McEvoy MT, Dann FJ, Gabriel SE, Kremers HM. Trends in Incidence of Adult-Onset Psoriasis Over Three Decades: A Population-Based Study. J Am Acad Dermatol (2009) 60:394–401. doi: 10.1016/j.jaad.2008.10.062

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Xie S, Chen Z, Wang Q, Song X, Zhang L. Comparisons of Gene Expression in Normal, Lesional, and Non-Lesional Psoriatic Skin Using Dna Microarray Techniques. Int J Dermatol (2014) 53:1213–20. doi: 10.1111/ijd.12476

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Rendon A, Schäkel K. Psoriasis Pathogenesis and Treatment. Int J Mol Sci (2019) 20:1475. doi: 10.3390/ijms20061475

CrossRef Full Text | Google Scholar

6. Dreyer LN, Brown GC. Oral Manifestations of Psoriasis. New Y State Dental J (2012) 33(4):14–8. doi: 10.1111/j.0303-6987.2006.00516.x

CrossRef Full Text | Google Scholar

7. Langley R, Krueger G, Griffiths C. Psoriasis: Epidemiology, Clinical Features, and Quality of Life. Ann Rheumatic Dis (2005) 64:ii18–23. doi: 10.1136/ard.2004.033217

CrossRef Full Text | Google Scholar

8. Pariser DM, Bagel J, Gelfand JM, Korman NJ, Ritchlin CT, Strober BE, et al. National Psoriasis Foundation Clinical Consensus on Disease Severity. Arch Dermatol (2007) 143:239–42. doi: 10.1001/archderm.143.2.239

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Raut AS, Prabhu RH, Patravale VB. Psoriasis Clinical Implications and Treatment: A Review. Crit Rev Ther Drug Carrier Syst (2013) 30:697. doi: 10.1615/CritRevTherDrugCarrierSyst.2013005268

CrossRef Full Text | Google Scholar

10. Lee EB, Wu KK, Lee MP, Bhutani T, Wu JJ. Psoriasis Risk Factors and Triggers. Cutis (2018) 102:18–20. doi: 10.1155/2013/537028

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Griffiths CE, van de Kerkhof P, Czarnecka-Operacz M. Psoriasis and Atopic Dermatitis. Dermatol Ther (2017) 7:31–41. doi: 10.1007/s13555-016-0167-9

CrossRef Full Text | Google Scholar

12. Yin L, Hu YY, Jameel AAB, Xu JL, Yin ZQ. Psoriasis and Skin Cancer. Biomed Res (2017) 28:2111–3. doi: 10.1189/jlb.1109767

CrossRef Full Text | Google Scholar

13. Takeda A, Higuchi D, Takahashi T, Ogo M, Baciu P, Goetinck PF, et al. Overexpression of Serpin Squamous Cell Carcinoma Antigens in Psoriatic Skin. J Invest Dermatol (2002) 118:147–54. doi: 10.1046/j.0022-202x.2001.01610.x

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Haider AS, Peters SB, Kaporis H, Cardinale I, Fei J, Ott J, et al. Genomic Analysis Defines a Cancer-Specific Gene Expression Signature for Human Squamous Cell Carcinoma and Distinguishes Malignant Hyperproliferation From Benign Hyperplasia. J Invest Dermatol (2006) 126:869–81. doi: 10.1038/sj.jid.5700157

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Gupta M, Das J, Gangopadhyay A. Multicentric Squamous Cell Carcinoma Arising on Psoriatic Plaque. Indian J Dermatol (2013) 58:151. doi: 10.4103/0019-5154.108065

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Asgari MM, Ray GT, Geier JL, Quesenberry CP. Malignancy Rates in a Large Cohort of Patients With Systemically Treated Psoriasis in a Managed Care Population. J Am Acad Dermatol (2017) 76:632–8. doi: 10.1016/j.jaad.2016.10.006

PubMed Abstract | CrossRef Full Text | Google Scholar

17. National Cancer Institute N. (1988). What You Need to Know About, Skin Cancer. NIH publication (U.S. Department of Health and Human Services, Public Health Service, National Institutes of Health, National Cancer Institute).

Google Scholar

18. Linares MA, Zakaria A, Nizran P. Skin Cancer. Primary Care (2015) 42:645–59. doi: 10.1016/j.pop.2015.07.006

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Burton KA, Ashack KA, Khachemoune A. Cutaneous Squamous Cell Carcinoma: A Review of High-Risk and Metastatic Disease. Am J Clin Dermatol (2016) 17:491–508. doi: 10.1007/s40257-016-0207-3

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Que SKT, Zwald FO, Schmults CD. Cutaneous Squamous Cell Carcinoma: Incidence, Risk Factors, Diagnosis, and Staging. J Am Acad Dermatol (2018) 78:237–47. doi: 10.1016/j.jaad.2017.08.059

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Bernstein SC, Lim KK, Brodland DG, Heidelberg KA. The Many Faces of Squamous Cell Carcinoma. Dermatol Surg (1996) 22:243–54. doi: 10.1111/j.1524-4725.1996.tb00315.x

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Johnson TM, Rowe DE, Nelson BR, Swanson NA. Squamous Cell Carcinoma of the Skin (Excluding Lip and Oral Mucosa). J Am Acad Dermatol (1992) 26:467–84. doi: 10.1016/0190-9622(92)70074-P

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Trakatelli M, Ulrich C, Del Marmol V, Euvrard S, Stockfleth E, Abeni D. Epidemiology of Nonmelanoma Skin Cancer (Nmsc) in Europe: Accurate and Comparable Data are Needed for Effective Public Health Monitoring and Interventions. Br J Dermatol (2007) 156:1–7. doi: 10.1111/j.1365-2133.2007.07861.x

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Housman TS, Feldman SR, Williford PM, Fleischer AB Jr., Goldman ND, Acostamadiedo JM, et al. Skin Cancer is Among the Most Costly of All Cancers to Treat for the Medicare Population. J Am Acad Dermatol (2003) 48:425–9. doi: 10.1067/mjd.2003.186

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Gray DT, Suman VJ, Su WD, Clay RP, Harmsen WS, Roenigk RK. Trends in the Population-Based Incidence of Squamous Cell Carcinoma of the Skin First Diagnosed Between 1984 and 1992. Arch Dermatol (1997) 133:735–40. doi: 10.1001/archderm.133.6.735

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Cassarino DS, DeRienzo DP, Barr RJ. Cutaneous Squamous Cell Carcinoma: A Comprehensive Clinicopathologic Classification: Part Two. J Cutaneous Pathol (2006) 33:261–79. doi: 10.1111/j.0303-6987.2006.00516.x

CrossRef Full Text | Google Scholar

27. Kwa RE, Campana K, Moy RL. Biology of Cutaneous Squamous Cell Carcinoma. J Am Acad Dermatol (1992) 26:1–26. doi: 10.1016/0190-9622(92)70001-V

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Hemminki K, Zhang H, Czene K. Familial Invasive and in Situ Squamous Cell Carcinoma of the Skin. Br J Cancer (2003) 88:1375. doi: 10.1038/sj.bjc.6600909

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Mullen JT, Feng L, Xing Y, Mansfield PF, Gershenwald JE, Lee JE, et al. Invasive Squamous Cell Carcinoma of the Skin: Defining a High-Risk Group. Ann Surg Oncol (2006) 13:902–9. doi: 10.1245/ASO.2006.07.022

PubMed Abstract | CrossRef Full Text | Google Scholar

30. MATLAB Optimization Toolbox. Matlab Optimization Toolbox. (2019).

Google Scholar

31. Klaus B, Reisenauer S. An End to End Workflow for Differential Gene Expression Using Affymetrix Microarrays. F1000Res (2018) 63(4):278–85. doi: 10.12688/f1000research.8967.2

CrossRef Full Text | Google Scholar

32. Klaus B, Reisenauer S An End to End Workflow for Differential Gene Expression Using Affymetrix Microarrays. F1000Research (2016) 5:1384. doi: 10.12688/f1000research.8967.1

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Gu Z, Eils R, Schlesner M. Complex Heatmaps Reveal Patterns and Correlations in Multidimensional Genomic Data. Bioinformatics (2016) 32:2847–9. doi: 10.1093/bioinformatics/btw313

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Hicks SC, Okrah K, Paulson JN, Quackenbush J, Irizarry RA, Bravo HC. Smooth Quantile Normalization. Biostatistics (2017) 19:185–98. doi: 10.1093/biostatistics/kxx028

CrossRef Full Text | Google Scholar

35. Hicks SC, Irizarry RA. When to Use Quantile Normalization? BioRxiv (2014) 012203:ii18–23. doi: 10.1101/012203

CrossRef Full Text | Google Scholar

36. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma Powers Differential Expression Analyses for Rna-Sequencing and Microarray Studies. Nucleic Acids Res (2015) 43:e47–7. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Lopez-Romero PA. Processing and Differential Expression Analysis of Agilent Microrna Chips. R Package Version (2016) 2:336–6. doi: 10.1049/iet-syb.2011.0015

CrossRef Full Text | Google Scholar

38. Dweep H, Gretz N. mirwalk2. 0: A Comprehensive Atlas of Microrna-Target Interactions. Nat Methods (2015) 12:697. doi: 10.1038/nmeth.3485

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Pertea M, Kim D, Pertea GM, Leek JT, Salzberg SL. Transcript-Level Expression Analysis of Rna-Seq Experiments With Hisat, Stringtie and Ballgown. Nat Protoc (2016) 11:1650. doi: 10.1038/nprot.2016.095

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Chen S, Zhou Y, Chen Y, Gu J. Fastp: An Ultra-Fast All-in-One Fastq Preprocessor. Bioinformatics (2018) 34:i884–90. doi: 10.1093/bioinformatics/bty560

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Kim D, Langmead B, Salzberg SL. Hisat: A Fast Spliced Aligner With Low Memory Requirements. Nat Methods (2015) 12:357. doi: 10.1038/nmeth.3317

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. Stringtie Enables Improved Reconstruction of a Transcriptome From Rna-Seq Reads. Nat Biotechnol (2015) 33:290. doi: 10.1038/nbt.3122

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Frazee AC, Pertea G, Jaffe AE, Langmead B, Salzberg SL, Leek JT. Flexible Isoform-Level Differential Expression Analysis With Ballgown. Nature Biol (2014) 33(3):243–6. doi: 10.1155/2013/569751

CrossRef Full Text | Google Scholar

44. Blighe K, Rana S, Lewis M. Enhancedvolcano: Publication-ready Volcano Plots With 667 Enhanced Colouring and Labeling (2020). Available at: https://github.com/kevinblighe.

Google Scholar

45. Xie Z, Bailey A, Kuleshov MV, Clarke DJ, Evangelista JE, Jenkins SL, et al. Gene Set Knowledge Discovery With Enrichr. Curr Protoc (2021) 1:e90. doi: 10.1002/cpz1.90

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Zi Z. Sensitivity Analysis Approaches Applied to Systems Biology Models. IET Syst Biol (2011) 5:336–46. doi: 10.1049/iet-syb.2011.0015

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Wang S, Song R, Wang Z, Jing Z, Wang S, Ma J. S100a8/a9 in Inflammation. Front Immunol (2018) 9:1298. doi: 10.3389/fimmu.2018.01298

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Schonthaler HB, Guinea-Viniegra J, Wculek SK, Ruppen I, Ximénez-Embún P, Guío-Carrión A, et al. S100a8-s100a9 Protein Complex Mediates Psoriasis by Regulating the Expression of Complement Factor C3. Immunity (2013) 39:1171–81. doi: 10.1016/j.immuni.2013.11.011

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Lee Y, Jang S, Min JK, Lee K, Sohn KC, Lim JS, et al. S100a8 and s100a9 are Messengers in the Crosstalk Between Epidermis and Dermis Modulating a Psoriatic Milieu in Human Skin. Biochem Biophys Res Commun (2012) 423:647–53. doi: 10.1016/j.bbrc.2012.05.162

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Hawkes JE, Yan BY, Chan TC, Krueger JG. Discovery of the il-23/il-17 Signaling Pathway and the Treatment of Psoriasis. J Immunol (2018) 201:1605–13. doi: 10.4049/jimmunol.1800013

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Brembilla NC, Senra L, Boehncke WH. The il-17 Family of Cytokines in Psoriasis: Il-17a and Beyond. Front Immunol (2018) 9:1682. doi: 10.3389/fimmu.2018.01682

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Nardinocchi L, Sonego G, Passarelli F, Avitabile S, Scarponi C, Failla CM, et al. Interleukin-17 and interleukin-22 Promote Tumor Progression in Human Nonmelanoma Skin Cancer. Eur J Immunol (2015) 45:922–31. doi: 10.1002/eji.201445052

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Gorczynski R. Il-17 Signaling in the Tumor Microenvironment. Tumor Microenviron (2020) 78(2):47–58. doi: 10.1007/978-3-030-38315-2_4

CrossRef Full Text | Google Scholar

54. Rahat MA, Shakya J. Parallel Aspects of the Microenvironment in Cancer and Autoimmune Disease. Mediators Inflamm (2016) 2016:17. doi: 10.1155/2016/4375120

CrossRef Full Text | Google Scholar

55. Giat E, Ehrenfeld M, Shoenfeld Y. Cancer and Autoimmune Diseases. Mosaic Autoimmun (2019) 30(3):453–65. doi: 10.1016/B978-0-12-814307-0.00041-4

CrossRef Full Text | Google Scholar

56. Uribe-Querol E, Rosales C. Neutrophils in Cancer: Two Sides of the Same Coin. J Immunol Res (2015) 2015:28. doi: 10.1155/2015/983698

CrossRef Full Text | Google Scholar

57. Gatfield J, Ferrari G, Pieters J. Immunity Against Extracellular Pathogens. Protoplasma (2000) 210:99–107. doi: 10.1007/BF01276849

CrossRef Full Text | Google Scholar

58. Gonzalez H, Hagerling C, Werb Z. Roles of the Immune System in Cancer: From Tumor Initiation to Metastatic Progression. Genes Dev (2018) 32:1267–84. doi: 10.1101/gad.314617.118

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Mavropoulos A, Rigopoulou EI, Liaskos C, Bogdanos DP, Sakkas LI. The Role of p38 Mapk in the Aetiopathogenesis of Psoriasis and Psoriatic Arthritis. Clin Dev Immunol (2013) 2013:1405–22. doi: 10.1155/2013/569751

CrossRef Full Text | Google Scholar

60. Haase I, Hobbs RM, Romero MR, Broad S, Watt FM. A Role for Mitogen-Activated Protein Kinase Activation by Integrins in the Pathogenesis of Psoriasis. J Clin Invest (2001) 108:527–36. doi: 10.1172/JCI12153

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Lee CH, Hwang STY. Pathophysiology of Chemokines and Chemokine Receptors in Dermatological Science: A Focus on Psoriasis and Cutaneous T-Cell Lymphoma. Dermatol Sin (2012) 30:128–35. doi: 10.1016/j.dsi.2012.08.004

CrossRef Full Text | Google Scholar

62. Singh TP, Lee CH, Farber JM. Chemokine Receptors in Psoriasis. Expert Opin Ther Targets (2013) 17:1405–22. doi: 10.1517/14728222.2013.838220

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Mezentsev A, Nikolaev A, Bruskin S. Matrix Metalloproteinases and Their Role in Psoriasis. Gene (2014) 540:1–10. doi: 10.1016/j.gene.2014.01.068

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Johansen C, Kragballe K, Rasmussen M, Dam T, Iversen L. Activator Protein 1 Dna Binding Activity is Decreased in Lesional Psoriatic Skin Compared With Nonlesional Psoriatic Skin. Br J Dermatol (2004) 151:600–7. doi: 10.1111/j.1365-2133.2004.06088.x

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Piruzian E, Nikolskaya T, Abdeev R, Brouskin S. Transcription Factor Ap-1 Components as Psoriasis Candidate Genes. Mol Biol (2007) 41:974–85. doi: 10.1134/S0026893307060143

CrossRef Full Text | Google Scholar

66. Sun Y, Tu Y, He L, Ji C, Cheng B. High Mobility Group Box 1 Regulates Tumor Metastasis in Cutaneous Squamous Cell Carcinoma Via the pi3k/akt and Mapk Signaling Pathways. Oncol Lett (2016) 11:59–62. doi: 10.3892/ol.2015.3843

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Eckert RL, Adhikary G, Young CA, Jans R, Crish JF, Xu W, et al. Ap1 Transcription Factors in Epidermal Differentiation and Skin Cancer. J Skin Cancer (2013) 2013. doi: 10.1155/2013/537028

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Wang X, Lin Y. Tumor Necrosis Factor and Cancer, Buddies or Foes? 1. Acta Pharmacol Sin (2008) 29:1275–88. doi: 10.1111/j.1745-7254.2008.00889.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: microarray data analysis, RNA-seq data analysis, psoriasis, cutaneous squamous cell carcinoma, pathway analysis, biomarkers, therapeutic target, differentially expressed genes

Citation: Adil S, Paracha RZ, Tariq S, Nisar M, Ijaz S, Siddiqa A, Hussain Z and Amir A (2021) A Computational Systems Analyses to Identify Biomarkers and Mechanistic Link in Psoriasis and Cutaneous Squamous Cell Carcinoma. Front. Immunol. 12:662528. doi: 10.3389/fimmu.2021.662528

Received: 02 February 2021; Accepted: 19 May 2021;
Published: 18 June 2021.

Edited by:

Eddie A. James, Benaroya Research Institute, United States

Reviewed by:

Jenan Almatouq, Mohammed Al Mana College for Health Sciences (MACHS), Saudi Arabia
Eric Munger, National Institutes of Health Clinical Center (NIH), United States

Copyright © 2021 Adil, Paracha, Tariq, Nisar, Ijaz, Siddiqa, Hussain and Amir. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Rehan Zafar Paracha, cmVoYW5AcmNtcy5udXN0LmVkdS5waw==

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.