A Transcriptomic Meta-Analysis Identifies Differentially Expressed Genes across Different Types of Cancer and Their Association with Checkpoint Inhibitor Immunotherapy

Background Checkpoint inhibitor, anti-Programmed cell death protein1/Ligand (PD1 / PDL1) immunotherapy achieved great success in modulating the immune system to reinvigorate targeted immune fight against several types of cancer. However, patients’ clinical response to this cancer immunotherapy varies widely. Biomarkers that can predict the greatest response to the anti-PD1 immunotherapy could help to personalize the maximum benefit for the right patient.Results This study explored transcriptomic biomarkers that were associated with the response (or no-response) to the anti-PD1/PDL1 immunotherapy in seven types of cancer Malignant Pleural Mesothelioma (n=8), Head and Neck cancer (n=3), Lung Non-Squamous cancer (n=15), Squamous Lung cancer (n=8), Melanoma (n=23), Melanoma Skin Metastasis (n=10), and Renal Cell Carcinoma (n=11). The most common differentially expressed genes (105 genes) were between melanoma skin metastasis and renal cell carcinoma; followed by malignant pleural mesothelioma and renal cell carcinoma, which had 79 common differentially expressed genes; then malignant pleural mesothelioma and melanoma skin metastasis, which had 46 differentially expressed genes in common. Many similar (parallel) and dissimilar (anti-parallel) gene expression signatures were identified in this study. Parallel signature Differentially Expressed Genes (DEGs) are genes up-regulated together in either response or no-response group of different types of cancer, while the anti-parallel DEGs show up-regulation in the response group of one type of cancer while same genes are up-regulated in the no-response group of another different type of cancers. Target of Myb1 like 1 Membrane Trafficking Protein gene (TOM1L1) was found to be common among malignant pleural mesothelioma, melanoma skin metastasis, and renal cell carcinoma and was up-regulated in the cancer patients’ samples that did not respond to the anti-PD1/PDL1 immunotherapy. This

Urokinase Receptor gene (PLAUR), was common among squamous lung cancer, melanoma and melanoma skin metastasis, and was up-regulated in the patients' samples that responded to the anti-PD1/PDL1 immunotherapy.Conclusion In summary, this study identified DEGs biomarkers that can predict the response or no-response to the anti-PD1 immunotherapy across seven types of cancer. It also affirms the attention to important genes like TOM1L1 and PLAUR for their potential function in the cancer dissemination and metastasis.

Background
In March 2015, the Food and Drug Administration (FDA) approved the first application of immunotherapy (1). Following this approval, an extensive number of cancer types have been treated using the same immunotherapeutic approach. With remarkable success in treating cancer patients, this targeted immunotherapeutic approach harvested the 2018 Nobel Prize in physiology or medicine (2). This targeted immunotherapeutic approach is being described as a new phase in the war to end cancer. The approach shifted the era from directly targeting the cancerous cells and tissues using surgery, radiation, and chemotherapy, to an indirect approach that unleashes the body's immune system from being constrained by PDL1 checkpoint molecules expressed on the cancer cells (3). This cancer immunotherapeutic approach reached higher success rates and greatly improved the overall survival and the quality of life of a certain group of cancer patients when compared to conventional cancer chemotherapy (4). Additional advantages of this immunotherapy come from the reduced toxicity and lack of immune suppression accompanied by chemo and radiation therapy (5). Yet, despite the simplicity, the elegant approach, and the high success rate of this immunotherapy in treating cancer patients, the approach exhibited a response rate that varied widely among different types of cancers (6).
With the advancement of microarray technologies and the decreasing cost of genomic and transcriptomic profiling methods, a vast amount of data has been created [7]. This gave rise and made easier precision oncology studies to help understand and appoint new therapeutic targets [7], in addition to predicting biomarker signatures that help decisionmaking for clinicians to tailor cancer treatments to the potentially most benefited group of patients [7]. In this study, we aim to investigate the common transcriptomic profiling signatures among seven types of cancer as they relate with the response (or lack of) to the anti-PD1 immunotherapy. By taking advantage of the publicly available transcriptomic profiling studies, this high-throughput, precise and intense approach would provide a vital model tool in the form of transcriptomic signatures to discriminate responses (or lack of) to anti-PD1 immunotherapy.

Study characteristics.
The National Center for Biotechnology Information (NCBI) stores a vast amount of raw microarray data from different studies in their embedded repository database called Gene Expression Omnibus (GEO) (7). The raw data of this meta-analysis was downloaded from the GEO database. This project includes microarray profiling data of cancer patients' samples from four different studies with GEO Series accession number (GSE #) (GSE99070(8), GSE93157(9), GSE79691(10), and GSE67501(11)) as shown in Table 1.
These four studies include transcriptomic data of cancer tissue samples for seven types of cancers (malignant pleural mesothelioma, head and neck, lung non-squamous, squamous lung, melanoma, melanoma skin metastasis, and renal cell carcinoma). The microarray chip technology used in three of these studies were Illumina High-Density Silica Bead-Based Microarray Technology, while one study used GPL19965 nCounter PanCancer Immune Profiling Panel, (Table 1) shows the specific chips used from each study set with the Geo Sample Accession Number (GSM #) identifier provided. These GSM chips were from patients' samples that had either complete/partial response (CR/PR) or noresponse/progressed disease (PD) to anti-PD1 immunotherapy. All other microarray chips that had stable disease status or any other status were excluded from this study. The included GSM samples that responded to the anti-PD1 immunotherapy in all seven types of cancer equaled 34 samples, while the GSM samples from patients that did not respond to the anti-PD1 immunotherapy equaled 44 samples.

Common DEGs identified across different types of cancer.
Common DEGs between any two or more types of cancer were identified and shown in Table 2. All DEGs were significantly dis-regulated (p-value 0.05). The full list of significantly dis-regulated genes for each type of cancer is provided in supplementary excel sheet (Data1e). Malignant pleural mesothelioma shared common genes with five types of cancer included in this study ( Figure 1A), yet, the number of common DEGs between the malignant pleural mesothelioma and other types of cancer varied. For example, malignant pleural mesothelioma shared one common gene with lung nonsquamous cancer, and also shared the highest common DEGs with renal cell carcinoma 79 genes ( Figure 1A). Six DEGs were shared in common among malignant pleural mesothelioma and two other types of cancers (melanoma skin metastasis and renal cell carcinoma) (Table 2 and Figure 1A).
The highest number of cancers to share common genes was four [squamous lung cancer, melanoma, melanoma skin metastasis, and renal cell carcinoma]. These four types of cancer shared one gene in common ( Table 2). The highest number of common DEGs shared was between melanoma skin metastasis and renal cell carcinoma (105 genes) ( Table 2).

Pattern of expression and potential signatures of the DEGs in this study.
To further detail and identify potential gene expression signatures associated with response or no-response to the anti-PD1 immunotherapy, we looked at the pattern of expression of common DEGs among different types of cancer included in this study. The Venn diagram in Figure 1A shows common DEGs shared between malignant pleural mesothelioma and five types of cancer that were previously detailed in part in Table 2.
The graph in Figure 1B shows the pattern of expression of the MX1 gene which was found to be common between malignant pleural mesothelioma and lung non-Squamous cancer.
The MX1 gene shows dissimilar (anti-parallel) pattern of expression, the negative logarithmic fold change value of any gene [MX1 as an example here (logFC= -1. 3

)]
indicates that the gene expression is up-regulated in the samples of patients that had complete response (CR) or partial response (PR) to the anti-PD1 immunotherapy compared to the patients' samples that had no-response (indicated as progressed disease "PD").  Figure 3D).
In addition to the one common gene that lung non-squamous cancer share with malignant pleural mesothelioma ( Figure 1B), this study also found that lung non squamous cancer share common genes with four other types of cancer [(squamous lung cancer, melanoma, melanoma skin metastasis, and renal cell carcinoma) two genes, two genes, one gene, and five genes respectively] ( Figure 4A). The patterns of expression signatures (similar and dissimilar) of these common genes between lung non-squamous and four other types of cancer are illustrated in Figures 4B, 4C, 4D, 4E and the inset table ( Figure 4F).
The head and neck cancer in this study share one common gene (CX3CL1) with renal cell carcinoma ( Figure 5A, 5C, and 5D), this gene shows dissimilar regulation signature in these two types of cancer. Therefore, CX3CL1 up-regulation could be an indicator of response in head and neck cancer, while the opposite can be interpreted from CX3CL1 which was up-regulated in the no-response group of the renal cell carcinoma ( Figure 5C).
Head and neck cancer, melanoma skin metastasis and squamous lung cancer share one gene in common (S100A7), however, S100A7 gene was up-regulated in squamous lung cancer patients which responded to the anti-PD1 immunotherapy ( Figure 5B), on the other hand, S100A7 gene in melanoma skin metastasis and head and neck was up-regulated in the no-response patients ( Figure 5B).
Four types of cancers are the most to share common genes in this study as shown in Figure 6B. Melanoma, squamous lung cancer, melanoma skin metastasis, and renal cell carcinoma they share PLAUR gene in common. PLAUR gene up-regulation is associated with the response to the anti-PD1 in all types of cancer except for the melanoma skin metastasis (associated with no-response as shown in Figure 6D).
This study also highlights a noteworthy gene expression signatures that is shown between melanoma and renal cell carcinoma (six genes as shown in Figure 6A). All six genes have similar up-regulation pattern of gene expression and associated with the response to anti-PD1 immunotherapy ( Figure 6E, 6F). On the other hand, melanoma and melanoma skin metastasis shared seven common DEGs. All seven genes had dissimilar patterns of gene expression ( Figure 6C, 6F). The seven genes in melanoma were up-regulated in the response patients, while same seven genes in melanoma skin metastasis were upregulated in the no-response patients.
Melanoma skin metastasis and renal cell carcinoma share the most common DEGs in this study (105 genes) ( Figure 6A). Among the 105 genes, 28 up-regulated DEGs were associated with response to anti-PD1 in both melanoma skin metastasis and renal cell carcinoma ( Figure 7A, and appendix 2A). Another gene expression signatures of 26 up-regulated DEGs were associated with no-response in both melanoma skin metastasis and renal cell carcinoma ( Figure 7B, and Appendix 2B). In addition to the two similar DEGs signatures, two dissimilar gene expression signatures were also identified between melanoma skin metastasis and renal cell carcinoma. 24 genes were significantly upregulated and associated to the response to anti-PD1 in melanoma skin metastasis sample, while same 24 genes were up-regulated in the renal cell carcinoma but were associated to the no-response ( Figure 7C, and Appendix 2C). The opposite gene expression signature pattern which comprised 27 genes was identified between melanoma skin metastasis and renal cell carcinoma. ( Figure 7D, and Appendix 2D).

Discussion
The anti-PD1 immunotherapy has improved survival and quality of life of cancer patients and is giving hope to many more as a new first line of treatment against cancer. With all the success anti-PD1 inhibitors have achieved so far, the benefit of this immunotherapy could be surged by increasing our understanding of the genetic makeup of responders across tumor types, as well as by identifying tumor genomics biomarker models that could be linked to the best response to the treatment.
The gene expression microarray approach is a high-scale powerful tool that has been used in many cancer studies. It precisely identifies a wide range of gene expression data that has allowed us to understand the fundamental mechanisms of drugs responses(12), including the new anti-PD1 immunotherapy. In this study, we uncovered transcriptomic signatures similarities and differences across different types of tumor and linked them to the prognosis of the disease in response to anti-PD1 immunotherapy.
The late-stage (meta-analysis) data integration approach utilized in this study (12), where the preprocessing and normalization of each experiment is implemented before, then merging the results of the analyses at the end step (12). This approach would minimize potential background discrepancies as a result of variation of the platforms and processing method of each study experiment (preparation and quality of the mRNA used, cDNA quality and the methodological variation in the hybridization) (13).
This study also illustrates important tumor gene expression patterns that were either parallel (similar, i.e. genes are up-regulated or down-regulated together either in the response or no-response group across more than one type of cancer), or anti-parallel (dissimilar, i.e. genes are up-regulated in the response group of one type of cancer while up-regulated in the no-response group of another type of cancer). Either parallel or antiparallel signatures would be invaluable in segregating target patients into a potential response or no-response groups. Generalizing results from this study, although gene enrichment and functionality of these signatures has not been addressed in this study, previous studies suggest that genes of parallel expression signature likely be of a closer functionality than those that are expressed in anti-parallel directions (14). Cancer studies suggest that cancer invasiveness and metastatic ability are driven through changes of many gene expression signatures within the tumor cells (15). Many of these studies highlighted signatures that are associated either with the primary or metastatic cancer to improve prognosis (16). Our study is in line with these previous studies and potentially uncovers gene expression signatures between for example primary melanoma and melanoma skin metastasis ( Figure 6C, the signature shows a perfect example of antiparallel gene expression pattern which in our study suggests being a perfect biomarker to tailor the anti-PD1 immunotherapy). Our study also supports the notion of changes in gene expression signatures that appear to facilitate primary cells to disseminate and metastasize (15).
Another good example in our study to support the statement in the previous paragraph is the expression signature of PLAUR gene. This gene is common among four types of cancer (squamous lung cancer, melanoma, melanoma skin metastasis, and renal cell carcinoma, shown in Figure 6D, 6F). Yet, it shows the anti-parallel pattern only in the metastatic type of melanoma (melanoma skin metastasis). The PLAUR gene has been highlighted in many studies for its importance in regulating extracellular proteolysis signaling pathways (17). It has also been suggested that PLAUR is associated with resistance to targeted immunotherapies (18). Our study result shows that PLAUR is associated with response to anti-PD1 immunotherapy in 3 primary cancer types, yet is associated with no-response to anti-PD1 when it comes to a metastatic type of cancer (melanoma skin metastasis in our study).

Conclusions
In summary, the meta-analysis used in this study identified DEGs signatures associated with response or no response to the anti-PD1 immunotherapy across different types of cancer. The DEGs signatures highlighted in this study could work as a valuable biomarker in tailoring anti-PD1 immunotherapy to the right patient cancer. The study also determined genes of potential functionality in cancer metastasis and resistance to targeted immunotherapy.

Selection of gene expression datasets.
Nivolumab is the generic name of the human programmed death receptor-1 (PD-1) blocking antibody (3), "Nivolumab" used as a keyword to search for transcriptomic data in the main repository of microarray data Gene Expression Omnibus database (GEO database) (19). The search retrieved 8 series and a total of 75 samples. The selection of the studies was based on the type of samples used in each study found, we aimed to include studies that has samples taken from a whole cancer tissue, all other studies for example "single cells isolated from cancer tissues, only immune cells isolated from cancer tissue samples, etc." ... were excluded from this study. The total number of studies included in our meta-analysis was four as GEO study accession numbers (GSE#) shown in Table1. Within each study selected, only samples that were associated with complete (CR) or partial response (PR) were selected to be compared against samples that were associated with no-response referred in the studies as Progressed Disease (PD), all other sample that had naive treatment, stable disease etc. status were excluded. GEO sample accession number (GSM#) of all samples included in this meta-analysis is provided in (Table1).

Meta-analysis of microarray datasets.
A late-stage data integration (meta-analysis) approach was utilized in this study as described previously (12,20). Briefly, in the analysis step, each microarray study data were analyzed individually, for combining results step, the derived statistically significant DEGs were combined and compared to find commonly shared genes using a Perl program (Combine_MA.pl) provided in (Supplementary Excel, Combine_MA.pl.xls).
The NCBI embedded tool known as GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r/) was utilized to analyze the included data sets of each study. The GEO2R is a bridging module between the R statistical package LIMMA (Linear Model for Microarray Analysis) and the NCBI repository database. LIMMA (https://www.bioconductor.org/packages/devel/bioc/manuals/limma/man/limma.pdf) was used to screen and identify DEGs between the [response to anti-PD1] versus [no response to the anti-PD1] samples. Results from each study stored as a .xls file, all results then combined into one excel sheet provided in the Supplementary excel file (sheet Data1e).
The Data1e file then saved as a text file and loaded to the Combine_MA.pl program to be processed for identify common DEGs across the seven types of cancer included in this study. As the microarray chip analysis considered a first step in identifying genes of potential functional or diagnostic importance (21), and to increase the power of identifying these subsets of genes, multiple testing corrections was not employed in this experiment and all common DEGs with p-value<0.05 were considered and analyzed as a potential gene expression signature (21,22

Ethics approval and consent to participate
The patient's microarray data used in this study is publically available in the GEO dataset "NCBI" database, therefore, no committee approvals were needed to implement this study.

Consent for publication
Not applicable.

Availability of data
As above, all data in this study is publically available.