Towards Understanding the Key Signature Pathways Associated from Differentially Expressed Gene Analysis in an Indian Prostate Cancer Cohort

1. Introduction

In the recent past, there is an increased occurrence of Prostate Cancer (PCa) varying among different races and countries, particularly between western (Caucasian) and eastern (Asian) populations [1] [2]. Over the last few years, there has been an increase in PCa cases across all the different parts and regions in India [3] [4]. Next Generation Sequencing (NGS) has advanced understanding of the diseased phenotypes, particularly cancers with tumor heterogeneity and capture of large amounts of genetic architecture leading to isoforms, splice junctions, post transcriptional modifications with a great accuracy and sensitivity [5] [6]. Although studies have dealt with the genetics, genomics, and the environmental influence in causal of PCa, no association of genotype and phenotype employing the NGS approaches has been discussed in the Indian population. Recent studies from our lab using the whole-exome sequencing (WES) approach of PCa in India have yielded very promising results to identify causal genes specific to PCa with 30 causal genes to the Indian cohort and over 8 genes specific to Indian phenotype with variable degree of genetic disposition [7]. Whole transcriptome shotgun sequencing (WTSS) or RNA-sequencing (RNA-Seq) on the other hand is used to study regulation, developmental biology, clinical and health sector [8]. The expected outcome of RNA-seq is identification of differentially expressed genes (DEGs) in different conditions (for eg. wild type vs mutant or control vs tumor) which could provide a detailed mechanism of an underlying disease [9]. The DEGs are selected based on the cut-off and log2 fold change that depends on p-values obtained by statistical modeling [10]. RNA-Seq generates millions of short reads that are aligned to a reference genome using different alignment softwares and based on this, the characteristic of a particular dataset is calculated [11]. Although RNAs are labile and unstable in alkaline conditions, they can be easily detected and quantified at very low abundance for varying gene expression partners, which makes them suitable to use as biomarkers [12]. Compared with DNA and protein biomarkers, RNA biomarkers have more sensitivity and specificity and are very cost-effective [13]. Similarly, RNA biomarkers have the advantage of providing dynamic insights into cellular states and regulatory processes than DNA biomarkers [14]. With the advent of high-throughput sequencing technologies, different types of non-coding RNAs (e.g., small nuclear RNA, micro RNA, small nucleolar RNA and long non-coding RNA etc.) and protein–coding RNAs (i.e., mRNAs) have been detected [15]. Interestingly, there are lots of novel non-coding RNAs discovered recently out of which mRNAs, piwiRNAs, siRNAs, ceRNAs and miRNAs have been well documented as diagnostic and prognostic markers in different types of cancers (ovary, lung, breast, colorectal) [16]. In addition to small RNAs, long noncoding RNAs (lncRNAs) previously considered as ‘transcriptional noises’ have been known to have a diverse and significant role in various diseases, primarily cancers [17] [18]. Although well-known lncRNAs in the form of MALAT, HOTAIR, XIST have been reported in genitourinary cancers, the increasing number of cases of PCa worldwide and as well Indian population warrants the need for identifying lncRNA based biomarkers [19] [20]. Further with the advent of RNA-seq, identification of lncRNAs will hopefully facilitate the translational research to the bench side [21]. Therefore, with considerable data available on the Indian population, transcriptome analysis on Indian patients supplementing strategies in diseased phenotypes of PCa will be of great interest [22]. In the current study, we analyzed the RNA-seq data from a cohort of 6 patients and identified DEGs as well as lncRNAs which could lay emphasis as deterministic markers. We discuss the impending results obtained from this study which we attempt to further compare and validate from our downstream analysis.

2. Materials and Methods

2.1. Patients, clinical samples and criteria

As a pilot study, from a cohort of 60, we have screened 6 patients who underwent prostatectomy (4 cases and 2 controls) from Rukmani Birla Hospital (RBH) (Table 1). The inclusion criteria include age > 55 years, non-diabetic or any other metabolic diseases, PSA > 4ng/ml, malignant while the exclusion criteria included the persons with smoking or familial history of benign prostate cases. The control inclusion on the other hand included age above 55, Gleason score <6, PSA <4ng/ml, normal/Benign Prostatic Hyperplasia (BPH) and normal BPH as exclusion criterion. The study carried out through our CA Prostate consortium of India (CAPCI; htttps://bioclues.org/capci last accessed on December 30, 2022) received approval from Institutional Ethics Committee (IEC) of Rukmani Birla Hospitals, Jaipur, India, and informed consent was judiciously taken.

2.2. Tissue preparation and RNA sequencing

RNA was isolated using the RNeasy FFPE kit (Qiagen, Catalog No-73504) from BPH and malignant Formalin-Fixed Paraffin Embedded (FFPE) blocks and sent for sequencing (outsourced). From an approximate 0.5 mg of cross section, the RNA was prepared wherein the NEBNext® Ultra™ II Directional RNA Library Prep Kit was used for preparing the libraries following the manufacturer’s protocol. While 100ng of FFPE RNA was used as input, it was then subjected to end repair and Illumina specific adaptors were ligated. The adaptor ligated product was then barcoded and subjected to 15 cycles of PCR. The samples after PCR were cleaned up using AMPure XP beads with the final libraries checked for quality using Qubit Fluorometer and Agilent Tapestation. The obtained libraries were pooled and diluted to final optimal loading concentration before cluster amplification on the Illumina flow cell. Once the cluster generation is completed, the cluster flow cell is loaded on Illumina HiSeq X instrument to generate 60M, 150bp paired end reads.

2.3. Bioinformatics and downstream RNA-Sequencing analysis

Cufflinks-Cuffdiff pipeline was employed to yield significant changes at the level of transcript expression [23], as we used our benchmarked pipeline from our lab to run through the workflow [24]. The sequences were aligned to human genome reference (build hg38) using HISAT2 to produce the alignment results output in SAM (sequence alignment map/file). HISAT2 aligns a set of unpaired reads (in fastq or .fq format) to the reference genome using the Ferragina and Manzini (FM)-index [25]. Cufflinks uses this map (SAM) and assembles the reads into transcripts, estimates their abundances, and finally examines DEGs from the samples. This was followed by Cuffdiff to check the DEGs which compared the aligned reads from RNA-seq samples from two or more conditions and identified transcripts that are differentially expressed using a rigorous FPKM normalization/statistical analysis [26]. These tools have gained worldwide attention and have been used in a number of transcriptomics studies. In the current study, Cuffdiff was used to perform differential analysis between the control samples and the other 4 malignant samples respectively [27]. A p-value cutoff of 0.05 and less was used to identify the significantly expressed transcripts (Figure 1). We performed a real-time PCR validation for some selected DEGs. The primers were ordered accordingly and RT-PCR was performed. We saw a significant difference in the expression of a few genes in malignant samples compared to control ones. But, some of the genes did not show any difference which could be due to small sample size as well as poor quality of FFPE blocks.

2.4. Interaction networks, statistical analysis, gene ontology and cbioportal analyses

We generated an interaction network considering a flexible and intuitive approach for generating gene function hypotheses, evaluating gene lists, and choosing genes for functional studies [28]. To check this, GeneMANIA was employed which overcomes the limitations of earlier methods for predicting gene function on yeast and animal benchmarks using GeneMANIA (https://genemania.org/ last accessed on December 26, 2022) [36]. In addition, we also used the search tool for the retrieval of interacting genes (STRING) database (https://string-db.org/), and integrated DEGs to STRING to evaluate the interaction. Experimentally valid interactions with a score of minimum 0.4 (40% or more) was chosen to be an ideal one with the resulting file saved as a tab-separated values (TSV) file [29]. The raw data files from the STRING database were then imported into Cytoscape 3.5.1 and the cytoHubba (http://hub.iis.sinica.edu.tw/cytohubba/ last accessed on December 22, 2022) plugin was employed with clustering coefficient, betweenness and closeness centralities to calculate the significant modules in the PPI network [30]. Different Cytoscape plugins can score and rank the nodes using different algorithms [31]. CytoHubba is one such plugin which uses a simple interface to analyze the different networks. CytoHubba implements eleven nodes such as degree, betweenness, closeness, clustering coefficient, stress etc. to rank any network [32]. PANTHER (http://www.pantherdb.org/ last accessed on December 22, 2022) is a gene ontology based functional annotation tool which takes a variety of inputs such as Gene IDs, UniProtKB IDs, Ensembl IDs etc. and results in either functional analysis or statistical enrichment analysis [33]. Furthermore, we deemed to check the expression and mutational profile of some of the genes from our study in cbioportal (https://www.cbioportal.org last accessed on December 22, 2022) for Cancer Genomics provides visualization, analysis and download of large-scale cancer genomics data sets [34]. We used data from TCGA, Pan cancer Atlas of Prostate Adenocarcinoma where 489 samples/patients were screened [35]. The statistics for the aforementioned analyses were based on network analyses with every approach having a function F attached to it that gives each node v a numerical value (v). If a node’s score—F(u)—is higher than another node’s score—F(v) we say that node u has a higher ranking than node v [32] [38].

3. Results

3.1. Distinct DEGs were obtained

Although there has been a difficulty in extracting RNA from FFPE blocks which is well known [37], from the RNA-seq and downstream analyses, all the samples yielded good quality reads from FastQC with ca. 40M transcript reads with no exposure in tiles. With ensuing cufflinks-cuffdiff pipeline, we obtained approximately 70 DEGs among which 65 were up-regulated and 5 were down-regulated with an inherent p-value heuristics <=0.05 and <=-2 Log2FC >=Log2FC. (Table 2).

Among the top niche specific DEGs, collagen type VI α1 chain (COL6A1), a gene which is located on chromosome 21, encoding the α1 (VI) chain of type VI collagen (which is a primary extracellular matrix protein) was found which maintains the integrity of various tissues. Catenin delta-1 (CTNND1) functions as an oncogene and is known to be the driver of metastatic cancer progression [40]. A very important gene which we identified is six-transmembrane epithelial antigen of Prostate-2 (STEAP2), known to be over-expressed in aggressive PCa, which corroborates our study as it was identified in a high grade tumor [44]. Furthermore, we also observed that docking protein 5 (DOK5), a member of a subgroup of the DOK family, is known to be expressed using c-Ret in several neuronal tissues [49]. Amongst the lncRNAs, we identified SOX2OT mapped to the chromosome locus 3q26.3 and is highly expressed in embryonic stem cells [55]. The lncRNA FTX (five prime to XIST) possesses an X-inactive specific transcript and is involved in X-chromosome inactivation [57].

3.2. Protein interactions yielded innate pathways responsible for PCa

From the input genes shown with cross-hatched circles of uniform size, GeneMANIA added relevant genes which are shown with solid circles and their size is proportional to the number of interactions they have (Figure 2). In addition, we have identified several interacting partners that are coexpressed such as DZIP1, COL6A2, TAGLN, ZBTB33, LAMP1, LAMP2, IKBKB, DNAJB11 which acts as oncogenes in different cancers such as gastric cancer, renal cell carcinoma, colorectal cancer, to name a few.

When we used the Cytoscape-cytoHubba plugin, top networks were ranked by clustering coefficient with the top 20 hub genes in our different malignant samples . Some of the prominent genes we identified through the networks are DOK5, APP, CTNND1, STX6, STX10, STX16, BACE1 and BACE2 (Figure 3 A-D)) which are in agreement with regulation of various cancers based on their coexpression patterns in GeneMANIA.

3.3. Validation of RNA-seq result using TCGA dataset by cbioportal

We sought to ask whether any DEGs were relatively expressed in publicly available datasets from various studies. To check this, we used TCGA datasets and checked for alteration of frequencies and expression in the cohorts. We used the PanCancer Atlas [39] dataset for Prostate Adenocarcinoma, where 489 samples were screened for their functional roles and molecular aberrations (Figure 4 and Figure 5). Mutations and CNV analysis for some of the DEGs were done and the results are summarized below:

Furthermore, an attempt was made to analyze genomic alterations such as gene amplifications, deep deletions (that is equivalent to homozygous deletions), shallow deletions (heterozygous loss), and truncating mutations, inframe mutations, or missense mutations. Among the DEGs, amplification was the most prominent one in case of DOK5 and STEAP2 whereas deep deletions and mutations were observed in COL6A1, STX6 and CTNND1. We argue that similar analysis could check performance of all the DEGs which will highlight the alteration frequencies across the cohort (Figure 5).

3.4. Gene ontology yielded distinct pathways regulating biogenesis

All the DEGs obtained from our cuffdiff pipeline were subjected to GO analysis by Pantherdb which shows the role of DEGs in biological adhesions, biological regulations, biogenesis, cellular, metabolic and developmental processes, localization, locomotion, multicellular organismal process (Figure 6 and Figure 7).

For the molecular function category, the terms are binding factors, catalytic, adapter, transducer, transcription regulator activities which are associated with differential expressed genes. This indicates the role of DEGs in different important processes like transcription, cell migration, differentiation etc.

Comparative analysis of RNA-seq data with other publicly available datasets

To screen the potential DEGs across the datasets, we compared and analyzed the DEGs from our current study to that of the DEGs from cbioportal and methylation/array specific datasets that are publicly available from NCBI gene expression omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/ last accessed December 22, 2022). As no concrete list of RNA-Seq datasets were associated in sequence read archive in lieu of PCa phenotype, we compared our DEGs to GEO datasets, viz. GSE6919 and GSE45016 in addition to previously benchmarked RNA-Seq dataset of PCa in the Chinese population (Ren et al, 2012). From this, we obtained only 135 hits in SRA for RNA-seq of PCa and most datasets are either from PCa cell-lines or RIP-seq and so we deemed them not useful keeping in view of diffident phenotypes, experiments and perhaps correlation studies they may be heralded with. Nevertheless, we identified a few common genes, viz. STEAP2, DOK5, Il6ST, LMOD1, CTNND1 etc., which could be the key candidate genes in our study. Likewise, when our data was compared with GEO datasets, we observed IL6ST, BACE2, SOX2-OT, STEAP2, APP, SNX9, STX16, CTNND1 among the other genes that are expressed. This, we believe, strengthens our finding that the DEGs which we have mentioned in the current study could be a valid signature for PCa diagnosis. On the other hand, when we compared Prostate Adenocarcinoma (TGCA, PanCancer Atlas) dataset to that of our list of DEGs, we found IL6ST, ZBTB20. These key DEGs are common between Indian population and western population which can be validated further as it is beyond the scope of this current analysis.

Phenolyzer highlights important DEGs

A cross-sectional comparison of phenotypes and DEGs using phenolyzer would provide us indicators for extent of expression across diseases. Therefore, we asked how many DEGs among all including the samples, cbioportal and comparative analyses are distinctly associated with disease/phenotype terms. Our Phenolyzer results with disease names, viz. Prostate Cancer, PCa neoplasia which are clinical phenotype terms have largely been connected with each other ( seen in pink edges in Figure 8). On the other hand, in an attempt to reach consensus with identification of DEGs, we also employed the DESeq/EdgeR normalization method which resulted in the identification of approximately 1,230 genes as DEG: 490 were up-regulated and 1,215 were down-regulated. We identified distinct genes including KLK4, FN1, PBOV1, TPM2, and FLNA which are known to be involved in PCa pathways along with IGF1, TPD52, and SRSF1 that are involved in different cancer pathways. While KLK4 is a very important gene which is known to be involved in the progression of prostate cancer by promoting proliferation, migration and epithelial to mesenchymal transition, we have checked its expression in our sample by qRT-PCR and we did significant change in benign vs malignant samples. We also found a few novel lncRNAs such as LINC00940, FLJ16779 that have not been reported earlier besides SNHG19, NPBWR1 lncRNAs that are known to be involved in cancer and other diseases as well. Further attempt was made to understand the regulatory mechanisms underpinning PCa signature pathways by mapping lncRNAs with protein encoding genes.

4. Discussion

In the current study, we have identified many DEGs which are known to be involved in different cancer pathways including PCa. While some of the DEGs are specifically known to be associated with PCa, we also discovered a few novel lncRNAs which need further investigation. Expression of COL6A1 is significantly elevated in different tumors such as lung, prostate, cervical and pancreatic cancer compared to normal tissues. Interestingly, our previous WES studies in PCa have identified COL6A1 as one of the causal genes [7], whereas, we also identified this through our RNA-seq analysis. COL6A1 was shown to be physically interacting with DNAJB11, APP along with some other genes. DNAJB11 is involved in aberrant signaling pathways associated with different cancers. Similarly, APP is known to be associated with androgen-responsive genes and regulates proliferation and migration of PCa cells. Therefore, we argue that COL6A1 might act as a prognostic marker for PCa in the Indian population. Transcriptional factor Kaiso/ZBTB33 was identified as a CTNND1-specific binding partner and this complex is a modulator of the canonical Wnt/β-catenin signaling pathway [41]. There is a large amount of research focusing on the role of CTNND1 in cancer development and progression, but in PCa, it is still not well elucidated [42]. In our study, we have identified CTNND1 and ZBTB33, through interaction studies but what is more interesting is that we have earlier identified CTNND1 as a co-localization partner with one of the lncRNAs which is known to be highly expressed in PCa. Furthermore, CTNND1 is one of the main interacting partners of ACE2/TMPRSS2, the main receptors which are responsible for SARS-CoV-2 entry into the cell. We had hypothesized that CTNND1 interacts with EGFR and this interaction could uphold SARS-CoV-2 infection independent of its endocytosis and associate with cell viability [43]. The interacting partners of another important gene STEAP2 are KLK3, KLK2 and AR all of which are hallmarks of PCa. KLK3 is a protein coding gene and its protein product Prostate specific antigen (PSA) is a well-established biomarker of PCa [45] [46]. Similarly, human kallikrein 2 (KLK2), interacts with AR and drives PCa progression [47] [48]. Earlier studies have shown that DOK5 are expressed in T-cells and their expression is regulated upon T-cell activation. DOK5 is shown to be involved in invasion and metastasis of cancer specifically in Gastric cancer but it has not been well studied in PCa [50]. Since, we got a strong interaction of DOK5 in our clustering coefficient studies using cytoHubba plugin, we argue that it would be worth analyzing this gene further. Interestingly, through our cytoscape-cytoHubba analysis, we identified many soluble N-ethylmaleimide-sensitive factor (NSF) attachment protein receptor (SNARE) proteins that are key mediators of membrane fusion [51]. All the SNARE proteins share a common sequence of 60–70-residue called the SNARE motifs, that helps in mediating the interaction between vesicle SNAREs (v-SNAREs) and target membrane SNAREs (t-SNAREs) [52]. One of the t-SNARE proteins, syntaxin 6 (STX6) is particularly important in vesicle fusion. STX6 is up-regulated in a variety of cancers including breast, colon, liver, pancreatic, prostate, bladder, skin, testicular, tongue, cervical, lung and gastric cancers and it has been identified as a common transcriptional target of p53 family members (p53, p63 and p73). Along with STX6, we have identified STX10, VTI1A and STX16 which can be further studied. Overall, we also found other important DEGs, viz. BACE1 and BACE2 which belong to a class of proteases called β-secretases that are extensively studied in Alzheimer’s disease [53]. Not many studies have been done with respect to their role in cancer, but there are some recent studies which have shown their involvement in pancreatic and skin cancers [54].

A major chunk of lncRNAs are novel and regulated in distinct pathways

LncRNAs do not code for proteins but they are involved in almost all biological processes such as gene expression, epigenetic regulation, cell cycle regulation etc in different cancers, PCa being one of them. Previous studies have highlighted the oncogenic role of lncRNAs in metastasis, proliferation and development of PCa whereas they act as tumor suppressants but still there are several lncRNAs whose functions are still unknown. With the advancement in NGS, bioinformatics analysis has enabled identification of many lncRNAs which show dysregulated expression in PCa [62]. Their diverse role has made them a target for all stages of PCa development which includes screening, diagnosis, prognosis and treatment further establishing their role as biomarkers in PCa. For example, MALAT-1 (Metastasis associated lung adenocarcinoma transcript 1), an lncRNA which is used to predict metastasis and survival in non-small cell lung cancer [63], has also been correlated with PCa development and progression [64]. It has also been reported that MALAT-1 expression closely correlates with PSA levels, Gleason scores, and tumor sizes [64]. Similarly, PCAT-18 (prostate cancer-associated non-coding RNA transcript 18 and SChLAP1 (second chromosome locus associated with prostate-1) has also been used as diagnostic and prognostic biomarkers in PCa [65] [66]. Some of the lncRNAs we identified from our current study includes LINC01440, SOX2OT, ENSG00000232855, ENST00000647843.1 and FTX. FTX has been reported to be involved in the tumorigenesis of multiple cancer types [58]. Long intergenic non-protein coding RNA 1440 (LINC01440) is novel lncRNA which needs to be explored for its role in cancers. A recent study has implicated its role in a spinal disorder known as Ossification of ligamentum flavum (OLF) where it was found to be up-regulated in the diseased patients compared to the healthy population [59]. Different studies have shown that SOX2OT acts as an oncogene and is elevated in different tumor types [56]. But its role and significance is still not explored in PCa, which makes it a very promising candidate for further analysis. Another lncRNA, ENSG00000287903 (NONHSAT106693) was earlier screened in our Vitamin K deficiency cohort from our recent study [61] wherein qRT-PCR showed a significant up-regulation in malignant samples compared to control samples (supplementary information) which provides us a raison d’etre to further validate the remaining lncRNAs as well. Given the role of lncRNAs in PCa, it would be interesting to see if any of these lncRNAs can serve as biomarkers, albeit several downstream experiments and well-designed clinical trials are to be employed. Taken together, there are a few limitations for our study. (i)The qRT-PCR validation needs more tumor and adjacent normal tissue samples as our sample size is very small. (ii) Some additional experiments, such as immunohistochemistry and Western blot, when performed could confirm the protein levels in PCa. (iii) Given the scarcity of PCa data in India, the survival analysis across the entire transcriptome could not be drawn as it is largely towards identification of biomarkers indicating the prognostic power but given the sample size conundrum we have, we are limited to check this power. Nevertheless, In our previous study, we have performed WES on PCa samples which we have cited but again the sample size is limited. (iv) Due to lack of fresh-frozen prostate/radical prostatectomy tissues, only FFPE blocks were used and isolating RNA from them is an arduous and challenging task [37], because of which we could not perform qRT-PCR for all the selected DEGs.

5. Conclusions

Prostate cancer cases are increasing in India even as NGS studies are just beginning to be explored. In our current RNA-seq study and subsequent bioinformatics analyses, we sought to characterize many DEGs including lncRNAs that are specific to PCa of Indian sub-population. While we identified some of the important genes, viz. DOK5, COL6A1, CTNND1, STEAP2, APP., their role in PCa is still not clear. We envisage that characterizing their functional aspects would help us understand PCa progression. Since search for non-invasive and more sensitive biomarkers is on the anvil across all solid tumors, we firmly hope that these lncRNAs amongst the DEGs would serve as a precedent in development of NGS panels for PCa detection in Indian phenotype.