Cross Comparison of Genes Influencing Deafness in Mice
Group U1: Zachary Beddingfield, Shraddha Bhutiyani, Janae Crawford, Nicholas Holly, Ross Paolucci
Introduction and Background
Mice are a common species used in analyzing human age-progressive hearing loss. DBA-2J mice are one model for analyzing age-progressive hearing loss, while CBA is considered a very good model of human presbycusis hearing loss (citation, Vasilyeva et al). We considered two studies that performed gene expression analysis on DBA-2J and CBA mice. The first experiment, GSE62173, analyzed DBA-2J mice using a treatment of L-methionine and valproic acid which has been experimentally shown to attenuate hearing loss (citation). L-methionine is a known epigenetic regulator that promotes methylation, a modification that inhibits gene expression, while valproic acid has the antagonistic effect of inducing demethylation (Ornoy et al). The second experiment, GSE49543, analyzed older CBA mice for whom presbycusis was identified compared to younger CBA mice without presbycusis symptoms and who served as controls.
We seek to investigate the genetic mechanism behind the L-methionine and valproic acid treatment that has been demonstrated to attenuate hearing loss in DBA-2J mice. We predict that, because L-methionine and valproic acid treatment have antagonistic effects on methylation that induce changes in gene expression, and because L-methionine and valproic acid are known to reduce the effects of hearing loss, the genes that are differentially expressed in treated DBA-2J mice and control CBA mice will be comparable, as will their resulting pathway analysis and PCA plots.
Datasets
| Number of Replicates | Age | Disease State |
| 4 | Young | Control |
| 10 | Middle-aged | Control |
| 4 | Old | Mild presbycusis |
| 4 | Old | Severe presbycusis |
| 5 | Young | Control |
| 7 | Middle-aged | Control |
| 5 | Old | Mild presbycusis |
| 2 | Old | Severe presbycusis |
Table 1. GSE62173
| Number of Replicates | Age | Protocol |
| 5 | 4 weeks | Untreated |
| 5 | 12 weeks | Control vehicle |
| 6 | 12 weeks | L-methionine and valproic acid |
Table 2. GSE49543
Analytic Pipeline and Workflow
Documentation of Statistical Tests
We first used ExAtlas software to identify differentially expressed genes in each dataset. We ran the function “Significant Genes” supplied by ExAtlas for expression profile data, controlling for false discovery rate using the following parameters: Z-value of 2, False Discovery Rate (FDR) of 0.05, and fold change of 1.5.
Hierarchical clustering and principal component analysis were performed within ExAtlas using the same parameters described above. Additionally, a correlation value of 0.7, a filtering value of “Filter rows and columns”, and a heatmap clustering value of “Hierarchical clustering” were used.
Summary Table of Differentially Expressed Genes
Differential expression analysis using the methods described above resulted in the following data:
| Series | Differentially Expressed Gene Count |
| GSE62173 | 4088 |
| GSE49543 | 29 |
Of these, 11 differentially expressed genes were found to be shared by both datasets, meaning that about 39% of those genes differentially expressed in the control CBA mice were also differentially expressed in treated DBA-2J mice. These were:
| Shared Genes | Function |
| Bpifb1 | Responsible for binding bacterial lipopolysaccharides (LPS) and modulating cellular response to LP (Hou et al 2004). |
| Crct1 | Poorly understood protein-coding gene. |
| Crisp1 | Responsible for helping spermatozoa achieve functional maturation while they move along the reproductive organ tract (Haendler et al 1993). |
| Dcpp1 | Shown to be a mouse ortholog of a human gene that is upregulated in pancreatic cancer and promotes tumor progression and metastasis (Song et al 2017). |
| Defb4 | Possesses antimicrobial activity targeting Gram-negative and Gram-positive bacteria alike (Rohrl et al 2008). |
| Defb6 | Possesses antimicrobial activity that is especially harmful to E. coli (Yamaguchi et al 2001). |
| Gabra6 | An inhibitory neurotransmitter in the vertebrate brain that promotes neuronal inhibition by opening a chloride channel (Kato 1990). |
| Krtdap | Hypothesized based on sequence similarity to act as a regulator of keratinocyte differentiation with roles in embryonic skin morphogenesis (Tsuchida et al 2004). |
| Lce1d | Poorly understood protein-coding gene. |
| Myh4 | Plays a role in muscle contraction (Acakpo-Satchivi et al 1997). |
| Sprr3 | The cross-linked envelope protein of keratinocytes (Steinert et al 1998). |
Hierarchical Clustering with Heat Map
Pathway Analysis Results
We initially performed pathway analysis on each series using David. Because the results of our pathway analysis were too large, visual representation was not possible using the selected software. To account for this, we used Reactome to generate visualizations of the pathways represented by the differentially expressed genes for each series (Figure 4, 5). Then, we used David to generate a list of the most overrepresented pathways, applying Bonferroni adjustment before filtering for the lowest Bonferroni scores. David identified 177 unique pathways overrepresented in Series 62173 (Table 5) and 6 pathways overrepresented in Series 49543 (Table 6), but for brevity, only the first 10 pathways for Series 62173 with the lowest Bonferroni scores are included.
Next, we compared the pathways that overlap between both series. One specific pathway was identified as being overrepresented in both series, “extracellular space”.
| Term | Count | % | Bonferroni |
| KW-0325~Glycoprotein | 693 | 17.5 | 3.82E-71 |
| KW-1015~Disulfide bond | 607 | 15.3 | 4.11E-69 |
| GO:0005576~extracellular region | 399 | 10.1 | 7.46E-51 |
| KW-0964~Secreted | 375 | 9.5 | 1.23E-49 |
| CARBOHYD:N-linked (GlcNAc…) asparagine | 658 | 16.6 | 7.91E-48 |
| mmu04080:Neuroactive ligand-receptor interaction | 138 | 3.5 | 8.44E-48 |
| GO:0005615~extracellular space | 370 | 9.4 | 5.72E-33 |
| TOPO_DOM:Extracellular | 422 | 10.7 | 3.45E-26 |
| KW-0732~Signal | 762 | 19.3 | 9.33E-26 |
| IPR017970:Homeobox, conserved site | 71 | 1.8 | 3.04E-21 |
Table 5. Pathway Analysis in David of Series 62173. Only the first ten identified pathways with the lowest Bonferroni scores are included.
| Term | Count | % | Bonferroni |
| KW-0514~Muscle protein | 5 | 17.2 | 1.38E-05 |
| GO:0005615~extracellular space | 12 | 41.4 | 6.36E-04 |
| KW-0417~Keratinization | 3 | 10.3 | 0.00753 |
| GO:0006936~muscle contraction | 4 | 13.8 | 0.011699 |
| KW-0964~Secreted | 8 | 27.6 | 0.015576 |
| KW-0518~Myosin | 3 | 10.3 | 0.027073 |
Table 6. Pathway Analysis in David of Series 49543.
PCA Results
Principal component analysis was performed on our data as described in the methods section. While this analysis provided insight for Series 49543, where 65.7% of the variance was described between PC1 and PC2, the visualization was ineffective for Series 62173, where PC1 and PC2 only described a combined 16.7% of variance (24.3% with PC3). Clustering can be observed in the PCA of Series 49543, with two distinct groups visually apparent.
Related Results to Other Studies and Published Literature
While our cross-analysis is novel, researchers have investigated each series individually. Miya et al. performed differential expression analysis on Series 62173 and identified 244 differentially expressed genes, a substantially smaller number than our analysis. They used approaches that limited the outputs of differential expression analysis that were beyond the scope of our project which resulted in a list of more pronounced differentially expressed genes. Miya et al. does not make specific conclusions about their findings beyond noting that the list of differentially expressed genes is important concerning the observed attenuating effects on progressive hearing loss of DBA-2J mice.
D’Souza et al. analyzed Series 49543, specifically looking at the GABA inhibitory neurotransmitter system. They found that genes within the GABA system were often upregulated or downregulated in CBA mice who experienced hearing loss. This finding correlates with our study, which identifies Gabra6, a subset of the GABA system, as differentially expressed in both Series 62173 and Series 49543.
Interestingly, Tadros et al. performed differential expression analysis of CBA mice with and without hearing loss and identified a different set of four genes that were neither acknowledged by D’Souza et al. nor identified by our analysis. Their gene functions do not appear to be correlated with the major pathways identified in Series 49543 in our research. We hypothesize that this is an artifact of differences either in protocol or statistical pipeline procedures. However, it could alternatively be a sign that CBA mice lose hearing for a variety of unrelated reasons. We find that the first possibility appears more likely.
Overall Results and Noteworthy Conclusions
We observed substantial overlap in shared differential expression between Series 62173 and Series 49543 mice, with 39% of those genes differentially expressed in Series 49543 mice also being differentially expressed in Series 62173 mice (Table 3). This finding supports our initial hypothesis, which was that the differential gene expression of Series 62173 DBA-2J mice who received L-methionine and valproic acid treatment would be correlated to the differential gene expression of Series 49543 CBA mice from the control group who did not yet age to the point of experiencing hearing loss. However, further analysis complicates these conclusions. Heatmap analysis did not reveal clear trends for either series (Figures 2, 3). While some shared pathways are visually apparent between the differentially expressed genes of both series (Figures 4, 5), we must also consider that Series 62173 has overrepresented pathways across all sections and would likely overlap well with randomly selected pathways from the list provided by Reactome.
Regardless, we expect that the overlap in pathways is not an artifact of chance but rather a result of the enormous level of variation in gene expression present in Series 62173 mice. Our analysis revealed a great deal of variation in gene expression in Series 62173 DBA-2J mice. We believe this is insightful to the effects of L-methionine and valproic acid treatment supplied to these mice, which are known to have antagonistic effects on methylation and subsequently gene expression (Ornoy et al 2020). These observations raise the follow-up hypothesis that the L-methionine and valproic acid treatment greatly alters gene expression across many aspects of DBA-2J mice and coincidentally provides benefits against the age-progressive hearing loss that they experience.
Works Cited
Acakpo-Satchivi, L. J. R., Edelmann, W., Sartorius, C., Lu, B. D., Wahr, P. A., Watkins, S. C., Metzger, J. M., Leinwand, L., & Kucherlapati, R. (1997). Growth and muscle defects in mice lacking adult myosin heavy chain genes. The Journal of Cell Biology, 139(5), 1219–1229. https://doi.org/10.1083/JCB.139.5.1219
D’Souza, M., Zhu, X., & Frisina, R. D. (2008). Novel approach to select genes from RMA normalized microarray data using functional hearing tests in aging mice. Journal of Neuroscience Methods, 171(2), 279. https://doi.org/10.1016/J.JNEUMETH.2008.02.022
Haendler, B., Kratzschmar, J., Theuring, F., & Schleuning, W. D. (1993). Transcripts for cysteine-rich secretory protein-1 (CRISP-1; DE/AEG) and the novel related CRISP-3 are expressed under androgen control in the mouse salivary gland. Endocrinology, 133(1), 192–198. https://doi.org/10.1210/ENDO.133.1.8319566
Hou, J., Yashiro, K., Okazaki, Y., Saijoh, Y., Hayashizaki, Y., & Hamada, H. (2004). Identification of a novel left-right asymmetrically expressed gene in the mouse belonging to the BPI/PLUNC superfamily. Developmental Dynamics : An Official Publication of the American Association of Anatomists, 229(2), 373–379. https://doi.org/10.1002/DVDY.10450
Kato, K. (1990). Novel GABAA receptor alpha subunit is expressed only in cerebellar granule cells. Journal of Molecular Biology, 214(3), 619–624. https://doi.org/10.1016/0022-2836(90)90276-R
Miya, F., Mutai, H., Fujii, M., Boroevich, K. A., Matsunaga, T., & Tsunoda, T. (2015). Gene expression profiling of DBA/2J mice cochleae treated with l-methionine and valproic acid. Genomics Data, 5, 323. https://doi.org/10.1016/J.GDATA.2015.06.022
Noben-Trauth, K., Zheng, Q. Y., & Johnson, K. R. (2003). Association of cadherin 23 with polygenic inheritance and genetic modification of sensorineural hearing loss. Nature Genetics 2003 35:1, 35(1), 21–23. https://doi.org/10.1038/ng1226
Ornoy, A., Becker, M., Weinstein-Fudim, L., & Ergaz, Z. (2020). S-Adenosine Methionine (SAMe) and Valproic Acid (VPA) as Epigenetic Modulators: Special Emphasis on their Interactions Affecting Nervous Tissue during Pregnancy. International Journal of Molecular Sciences, 21(10). https://doi.org/10.3390/IJMS21103721
Song, H., Song, J., Kim, Y. J., Jeong, H. H., Min, H. J., & Koh, S. S. (2017). DCPP1 is the mouse ortholog of human PAUF that possesses functional analogy in pancreatic cancer. Biochemical and Biophysical Research Communications, 493(4), 1498–1503. https://doi.org/10.1016/J.BBRC.2017.10.015
Steinert, P. M., Candi, E., Kartasova, T., & Marekov, L. (1998). Small proline-rich proteins are cross-bridging proteins in the cornified cell envelopes of stratified squamous epithelia. Journal of Structural Biology, 122(1–2), 76–85. https://doi.org/10.1006/JSBI.1998.3957
Tadros, S. F., D’Souza, M., Zhu, X., & Frisina, R. D. (2014). Gene Expression Changes for Antioxidants Pathways in the Mouse Cochlea: Relations to Age-related Hearing Deficits. PLoS ONE, 9(2). https://doi.org/10.1371/JOURNAL.PONE.0090279
Tsuchida, S., Bonkobara, M., McMillan, J. R., Akiyama, M., Yudate, T., Aragane, Y., Tezuka, T., Shimizu, H., Cruz, P. D., & Ariizumi, K. (2004). Characterization of Kdap, a protein secreted by keratinocytes. The Journal of Investigative Dermatology, 122(5), 1225–1234. https://doi.org/10.1111/J.0022-202X.2004.22511.X
Vasilyeva, O. N., Frisina, S. T., Zhu, X., Walton, J. P., & Frisina, R. D. (2009). Interactions of hearing loss and diabetes mellitus in the middle age CBA/CaJ mouse model of presbycusis. Hearing Research, 249(1–2), 44–53. https://doi.org/10.1016/J.HEARES.2009.01.007
Study of TSLP-Induced Neutrophil Expression in MRSA Infection
Agniruudrra R Sinha, Ari Jimenez, Jiyeong Choi, Yilin Lu, Jyothi Guruprasad
Background
Staphylococcus aureus is a gram-positive bacteria that can cause a series of infectious diseases in multiple human body organs. This pathogen can be transmitted within a community or hospital, and it is particularly dangerous if it develops antibiotic resistance. Despite multiple drugs having been developed, effective treatment still remains a challenge due to the presence of various multi-drug resistant (MDR) strains. One of the biggest challenges the current healthcare industry faces is multidrug-resistant bacteria. Generally, MDR bacteria are associated with nosocomial infections, however, some MDR bacteria have become quite prevalent causes of community-acquired infections. This spread of MDR bacteria into the community is of immense importance and is associated with increased morbidity, mortality, healthcare costs, and antibiotic use. Methicillin-resistant Staphylococcus aureus (MRSA) is one of the most common antibiotic-resistant bacterial pathogens in the world, responsible for about 171,000 infections each year in Europe [2]. In 2019, MRSA caused more than 100,000 deaths [1]. Because of the increasing prevalence of antibiotic bacteria such as this one, it is becoming more important to determine ways to effectively treat and protect individuals from infection.
Polymorphonuclear neutrophils (PMN) are the most abundant cells of the immune system, serving as the first line of immune defense against infection [3]. They entrap and destroy microorganisms through various methods such as phagocytosis, intracellular degradation, the release of granules, formation of neutrophil extracellular traps, and are also known to be mediators of inflammation[3].
Research Goal/Aim
Some of the first steps in designing treatment against a pathogen are to understand how it functions and interacts with its host. In an attempt to gain insight into this interaction, our study compares RNA-seq data between two data sets that studied neutrophil expression in the presence of Methicillin-resistant Staphylococcus aureus (MRSA) infection in Homo sapiens.
The 1st data set investigates the transcriptional suppression of the immune response from the virulence toxin saeRS produced from MRSA. This research found that during infection of S. aureus, the transcription level of the whole blood system was down-regulated. While the 2nd data set observed how the interaction of neutrophils with cytokine thymic stromal lymphopoietin (TSLP), which is highly expressed in keratinocytes of skin, can be used to increase the potency of neutrophil action (killing) on MRSA [4,5]. The aim of our study is to show the effect of cytokine thymic stromal lymphopoietin (TSLP) on neutrophil-associated degradation of methicillin-resistant Staphylococcus aureus through transcriptome analysis. Through this study, we hope to find out if there are any specific genes that get upregulated or deregulated by TSLP in neutrophils that enhance the degradation of MRSA.
Methods
Figure 1. Overview of the research
The datasets were downloaded from the NCBI’s Gene Expression Omnibus (GEO) database.
Dataset 1:
Series GSE193219: 6 samples (3 Controls, 3 Samples)
| GSM5776862 | D5_PMNs_Buffer (phagocyte infection) | GSM5776870 | D4_PMNs_WT (phagocyte infection) |
| GSM5776858 | D4_PMNs_Buffer (phagocyte infection) | GSM5776874 | D5_PMNs_WT (phagocyte infection) |
| GSM5776866 | D6_PMNs_Buffer (phagocyte infection) | GSM5776878 | D6_PMNs_WT (phagocyte infection) |
Dataset 1 contains fresh PMN samples. The control was treated with PMNs buffer, while the treated samples contained PMSs treated with S. aureus and the RNA sample was collected.
Dataset 2:
Series GSE73313: 4 samples (2 Controls, 2 Cases)
| GSM1890591 | DN1_CTL_4H | GSM1890592 | DN1_TSLP_4H |
| GSM1890599 | DN2_CTL_4H | GSM1890600 | DN2_TSLP_4H |
Dataset 2 samples contained RNA samples extracted from PMNs. The control group was treated with a salt buffer while the treated group was treated with S. aureus with TSLP.
The datasets were imported onto Galaxy for further analysis. On Galaxy, the data were normalized using Kallisto, followed by an analysis of differentially expressed genes using the DESeq2 function. The tabular data generated from this step was used for further analysis. Principal components analysis (PCA) was obtained as a result of DESeq2.
DESeq2 uses Wald test for hypothesis testing by taking parameters that have been estimated by maximum likelihood. DESeq2 implements the Wald test by taking the LFC and dividing it by its standard error, thus resulting in a z-statistic. The z-statistic is compared to a standard normal distribution, and a p-value is computed. Benjamini-Hochberg method was used in DESeq2 to obtain the adjusted p-values by multiplying each ranked p-value by m/rank [6].
With the tables generated on Galaxy, the heatmaps were made showing the dendrograms as well. The pathway analysis was done with Reactome, using the ssGSEA method. Other methods, PADOG, and Camera could be used to analyze dataset1 but they were not able to analyze dataset2 due to a lack of the number of samples.
Results
1. Summary of Differentially Expressed Genes
Figure 1. Bar graph of most significantly expressed genes in Dataset1
Figure 2. Bar graph of most significantly expressed genes in Dataset2
Upon preliminary investigation of the RNAseq data, it was seen that FTH1 and PAICS were most significantly expressed across dataset1 while SRR and SEPIN10 were more significantly expressed in Dataset2. Furthermore, there was no common gene significantly expressed in both datasets.
2. Principal component analysis cannot distinguish TSLP-treated samples from control samples.
Principal component analysis showed a clear variable that classifies the control group and the experimental group in dataset 1. While in the second dataset, the most significant variable cannot clearly distinguish the two groups. This result suggested that the treatment of TSLP could alter the gene expression level on the infected cells close to the normal condition. Combined with the previous research from the second dataset, we concluded that the gene expression level under uninfected conditions has more potency of neutrophil action (killing) on MRSA. And the TSLP treatment increased this potency by altering the gene expression pattern back to normal condition.
Figure 3. PCA plot Dataset1
| Treated | Control |
| Data25: D6_PMNs_WT | Data28: D6_PMNs_Buffer |
| Data26: D5_PMNs_WT | Data29: D5_PMNs_Buffer |
| Data27: D4_PMNs_WT | Data30: D4_PMNs_Buffer |
Data35: Reference – GRCh37_latest_rna.fna.gz
Figure 4. PCA plot Dataset2
| Control | Treated |
| Data1: DN1_CTL_4H | Data2: DN1_TSLP_4H |
| Data3: DN2_CTL_4H | Data4: DN2_TSLP_4H |
Data5: Reference – GRCh37_latest_rna.fna.gz
3. Hierarchical clustering and heat maps showed different gene expression patterns.
Based on the heatmap generated and clustered by hierarchical clustering, it was clear that the expression level did not significantly change in both of the datasets after the infection. Most of the genes were down-regulated for all of the samples compared to the reference genome. And only a few genes shown in the heat maps displayed an up-regulation pattern. For the first dataset, the samples were taken from three patients, and each patient contributed one control sample and one treatment sample. However, based on the first heatmap, we could not observe a consistent differential expression pattern among the three sample groups. This could suggest that the gene expression changes are unique to individuals. For the second dataset, the heat map showed a similar expression pattern from the samples within individual patients, which suggested that the TSLP treatment could help neutrophils maintain their gene expression level.
Figure 5. Heatmap of Dataset 1
Figure 6. Heatmap of Dataset 2
4. Pathway Analysis
The pathway analysis was done with Reactome, using the ssGSEA method. As ssGSEA is calculating the enrichment score (ES), the magnitude of changes should be considered in the analysis. The pathway analysis of dataset1 provided a visual summary with the list of pathways sorted by p-values (Figure5). Focusing on the most affected pathways, the lists of pathways are compared.
Figure 7. Visual summary and a list of pathways of Dataset1 as a result of pathway analysis
Figure 8, 9, and 10 show the same pathway in different datasets analyzed with different methods. Neutrophil degranulation is upregulated in the first dataset while cytokine signaling in the immune system is downregulated. GTP hydrolysis and joining of the 60S ribosomal subunit, SRP-dependent cotranslational protein targeting to membrane, Nonsense-Mediated Decay (NMD), and some other pathways are downregulated in both datasets. Otherwise, Defective SLC26A4 causes Pendred syndrome (PDS), sperm motility and taxes, signaling by RNF43 mutants, and other pathways are upregulated in both datasets.
Figure 8. Pathway analysis of Dataset1, using the Camera method
Figure 9. Pathway analysis of Dataset1, using the ssGSEA method
Figure 10. Pathway analysis of Dataset2, using the ssGSEA method
Discussion
Generally, our data did not observe significant differential expressed genes in both data sets. While there are still some interesting findings that could suggest the influence of TSLP treatment and the mechanism behind this. The first dataset can be easily clustered into two groups based on the most significant variable. However, there are overlaps between the two groups in the second dataset. This suggested that the treatment of TSLP can alter the gene expression level in neutrophils, keeping the gene expression level similar to uninfected samples. We also found that the changes in gene expression levels varied among different individuals. While most of the gene expression patterns maintain the same within the same patient, some of the genes displayed different expression levels without TSLP treatment. This is parallel with the previous research mentioned in the two datasets. Interestingly, samples with TSLP treatment displayed almost the same gene expression pattern within the same individual, which is consistent with the result of PCA that TSLP could maintain the gene expression under infection. For the pathway analysis, we expected the same pattern in neutrophil or cytokine-related pathways in both datasets. Even though the results showed some biological pathways in common, they were not involving PMNs or cytokine TSLP. We assume it might be due to the minor magnitude of difference in regulation or the limited number of samples we have. It may also be possible that TSLP affects the neutrophil environment and not the neutrophils themselves, hence not many pathways relevant to PMNs were found in pathway analysis. This is in accordance with the literature wherein it was mentioned that gene expression levels remained mostly unchanged with the addition of TSLP[4]. A study with enough samples for different methods of pathway analysis might bring meaningful results.
Overall, our data suggested that the maintenance of gene expression levels in the neutrophils is critical for keeping the killing activity of the cell. Therefore, in-depth research on how to maintain the gene expression level and the subsequent result should continue.
Reference
1. An estimated 1.2 million people died in 2019 from. University of Oxford. (2022, January). Retrieved November 4, 2022, from https://www.ox.ac.uk/news/2022-01-20-estimated-12-million-people-died-2019-antibiotic-resistant-bacterial-infections#:~:text=aeruginosa)%20led%20directly%20to%20929%2C000,between%2050%2C000%20and%20100%2C000%20deaths.
2. Larsen, J., Raisen, C. L., Ba, X., Sadgrove, N. J., Padilla-González, G. F., Simmonds, M. S., Loncaric, I., Kerschner, H., Apfalter, P., Hartl, R., Deplano, A., Vandendriessche, S., Černá Bolfíková, B., Hulva, P., Arendrup, M. C., Hare, R. K., Barnadas, C., Stegger, M., Sieber, R. N., … Larsen, A. R. (2022). Emergence of methicillin resistance predates the clinical use of antibiotics. Nature, 602(7895), 135–141. https://doi.org/10.1038/s41586-021-04265-w
3. Souto JC, Vila L, Brú A. Polymorphonuclear neutrophils and cancer: intense and sustained neutrophilia as a treatment against solid tumors. Med Res Rev. 2011 May;31(3):311-63. doi: 10.1002/med.20185. PMID: 19967776.
4. West, E. E., Spolski, R., Kazemian, M., Yu, Z. X., Kemper, C., & Leonard, W. J. (2016). A TSLP-complement axis mediates neutrophil killing of methicillin-resistant staphylococcus aureus. Science Immunology, 1(5). https://doi.org/10.1126/sciimmunol.aaf8471
5. Zwack, E. E., Chen, Z., Devlin, J. C., Li, Z., Zheng, X., Weinstock, A., Lacey, K. A., Fisher, E. A., Fenyö, D., Ruggles, K. V., Loke, P., & Torres, V. J. (2022). staphylococcus aureus induces a muted host response in human blood that blunts the recruitment of neutrophils. Proceedings of the National Academy of Sciences, 119(31). https://doi.org/10.1073/pnas.2123017119
6. Hawinkel, S., Mattiello, F., Bijnens, L., & Thas, O. (2019). A broken promise: microbiome differential abundance methods do not control the false discovery rate. Briefings in bioinformatics, 20(1), 210-221.
Validating the Efficacy of Anti-PD1/ PDL1 Therapies using Gene Expression Data & Analysis
Neha Jain, Akshita Singh, Anirudh Suri
Introduction
Checkpoints are immunological processes that inhibit T-Cell actions to control their cytotoxicity. Cancerous cells and pathways stimulate such checkpoints, thus, reducing the effects of T cells on neoplastic cells. Anti-PD1/PDL1 and anti-CTLA4 therapies are checkpoint inhibitors that have developed owing to this knowledge. Immune-checkpoint blockade (ICB) drugs have only recently begun to take prominence in the cancer therapy landscape and has the potential to revolutionize the way tumor malignancies are treated. The first ICB to be introduced in the market was ipilimumab in the year 2011.
Being a new and relatively unconventional therapy (compared to chemotherapy) there is still plenty inhibition surrounding the adoption of ICB as the main course of treatment (Korman et al., 2021)
Problem Statement
Owing to the rapid rise in bioinformatics and analysis tools we can penetrate into the vast biological datasets available that aid in broadening our fundamental understanding of complex biological phenomena and its implications with respect to healthcare and therapeutics.
We aim to compare the similarity scores of the patients’ RNA seq data before and during treatment to see if there is a significant difference across patients. The results of this analysis might give us some insight into the efficacy of the treatment.
We will also look at the genes with the most radical changes between pre and during-treatment and examine what these genes do in cancer treatment in general. By doing so, we can correlate this back with the immunology and see if there is a large number of these genes which have a positive outcome; we will use this data to validate the treatment effectiveness. We will be downsizing analysis to five patients’ data for computational purposes.
Datasets
For this analysis, we are sourcing all our data from one data source (GSE91061). This data set contains 109 RNAseq runs of a 65-patient cohort- 58 on treatment and 51 pre-treatment. We look to compare and analyze 5 patients’ on and pre-treatment RNAseq stemming from their Gene Expression Data. All the patients were treated with checkpoint inhibitory drugs (Anti PD1/PDL1 and/or CTLA4 therapies) as a cure for advanced malignancies, precisely, melanoma and small-cell lung cancer. All the FASTQ files were sourced using the Illumina Genome Analyzer at the Cleveland Clinic, OH, USA. The datasets were curated over a three-year tenure.
| Dataset/ Accession ID | SRA ID | Sample Detail |
|---|---|---|
| GSM2420331 | SRR5088885 | Pt9_Pre |
| GSM2420332 | SRR5088886 | Pt9_On |
| GSM2420359 | SRR5088913 | Pt31_Pre |
| GSM2420358 | SRR5088912 | Pt31_On |
| GSM2420375 | SRR5088929 | Pt34_Pre |
| GSM2420374 | SRR5088928 | Pt34_On |
| GSM2420301 | SRR5088855 | Pt62_Pre |
| GSM2420302 | SRR5088856 | Pt62_On |
| GSM2420313 | SRR5088867 | Pt92_Pre |
| GSM2420316 | SRR5088870 | Pt92_On |
Bioinformatics Analysis Pipeline
The datasets were first loaded as FASTQ files onto the Galaxy servers via their accession IDs from the GEO database. The abundance of RNA-Seq transcripts was then quantified using Kallisto Quant (a tool hosted on the Galaxy servers). The transcript files were modified to facilitate the mapping of the ENST target ID to ENSG. This data was then piped into the Deseq2 tool (also hosted on Galaxy) to determine the differentially expressed features/genes.
Tximport was used (along with python) to produce the differentially expressed gene matrix which was then normalized to run Principle Component Analysis (PCA) plots along with Hierarchical Clustering
To automatically remove outliers and genes whose means are below a certain threshold Deseq2 uses cool’s distance. Benjamini and Hochberg is also used to control for false discovery which is also used to calculate differentially expressed genes by adjusting the p-value. After obtaining a list of the differentially expressed genes, a PCA plot and heatmap were plotted to better visualize gene expression differences. The Database for Annotation, Visualization and Integrated Discovery (DAVID) was used to furnish pathways analysis and find genes that are predominantly expressed.
Data and Statistical Analysis
Our initial statistical analysis after getting raw counts from Kallisto was done using deseq2 from galaxy in order to get our differentially expressed genes. Deseq2 uses cook’s distance to automatically remove outliers and also removes genes if their base means are below a certain threshold. It uses Benjamini and Hochberg to control for false discovery which is also used to calculate gene differential expression through adjusting the p-value. It starts by calculating p-values to use to represent statistical significance in the gene differentiation of the different factor levels, on treatment and pre treatment, in the data sets. It also calculates the fold change, the standard deviation of error, and a Wald Stat test in order to compare the distances between the samples for multiple hypothesis testing. Deseq2 performs a Wald test instead of t-test for multiple hypothesis testing, we used an alpha value of .05. The p-values of certain genes were significant(had an alpha value of greater than .05), and we used these genes’ log fold change to determine the upregulated and down regulated genes and pathways. Before clustering we used tximport to generate a gene differential expression matrix. We compared the mean and medians of all the samples by running boxplot from matplotlib and we saw that the means and medians of the samples were pretty close, but there were quite a few outliers. Moreover, when we only looked at the genes with significant p-values, not all the means and medians were still sufficiently close[figures 2 and 3]. So we then normalized this data by scaling the features using z-scoring from scipy.stats and dropping any rows with low variance to the rest of the genes. This can all be seen in our jupyter notebook which are linked as google colabs in the section called Jupyter Notebooks below. We used sci-kit learn to perform both PCA and hierarchical clustering which will be described further in the next few sections.
Hierarchical Clustering and PCA
PCA is a technique to reduce the dimensionality of the data, and take data with multiple dimensions and reduce it down to a more interpretable set of principle components. The data set starts out very large and is reduced to these principal components. We ran a PCA with our principle components being our genes with significant p-values[figure 10] and a PCA with our principle components being our sample data[figure 5]. A new distance matrix is calculated on every iteration for each principal component. This distance matrix is calculated on our normalized data using Euclidean distance. For our PCA on the samples the first three principal components captured 67.19% of variance and our PCA on the genes with significant P-values found that the first three principal components captured 80.33% of variance(calculations in Jupyter Notebook linked below).
Hierarchical clustering is very sensitive to noise, thus we scaled our data before running hierarchical clustering as described in the data and statistical analysis section. We used complete linkage agglomerative clustering with Euclidean distance as the metric to create the distance matrix of the distance between the clusters to perform our hierarchical clustering. We then used the elbow method to choose the optimal number of k clusters as this method gives us the least amount of clusters we need in order to have each observation, gene or sample, closest to its mean using the sum of square distances of all the observations. Since we ran Hierarchical clustering for the genes with significant p-values and for all our samples separately we did this method twice as seen in figures 7 and 8 below. We found for our samples the optimal number of clusters was 8 clusters, and for our genes with significant p-values, our optimal number of clusters was 18. We first ran PCA on the samples[figure 5] and the principle components and genes with significant p-values as the principle components[figure 10] separately. We then highlighted the principal components which belong to the same k clusters in the same colors in figures 6 and 11. In figure 6, we can see three of the principle components that are samples are the same color (0. blue), we can see from the previous PCA plot[figure 5] that these are samples 92 pre, 92 on, and 9 on. This indicates that patient 92 had the least variance in their genomic data from their pre-trial and on-trial samples. Every other patient had their pre-trial and on trial data grouped in a unique separate cluster, showing us the means of each of these clusters were significantly different. The efficacy of the treatment was least on patient 92 from the paper these results were from (and they had the earliest mortality), patient 92 had the least amount of variance between their data pre and on drugs and least change in their genomic of their tumors. We also noticed patients 34 and 31 had the most amount of variance in our data between their on and pre trial tumor genomics, this was also the case in the paper we read. In the paper, as well as in our data, patient 62 had slightly more variation between the pre-trial and on-trial data then patient 9. Our heatmap resulting from our hierarchical clustering[figure 9] illustrates the variance in pre vs on trial tumor genomic data for all of our patients, with lighter colors indicating a larger variation in that specific cluster from either the patients overall averaged data. We also saw 91% of the genes with significant p-values were clustered in a cluster with a mean near the overall mean of the data in our PCA plot (calculations in jupyter notebook) in figure 10. Thus there were around 9% of these genes with significant variance among all the samples, these genes can be seen in the heatmap above, in which it is also shown for which samples they varied the most[figure 4].
matplotlib and sci-kit learn Clusters color coded, generated with matplotlib and sci-kit learn
Figure 10. PCA Plot of Genes with Significant p-value, Figure 11. PCA Plot of Genes with Significant p-value where hierarchical . generated with matplotlib and sci-kit learn clusters are color coded, generated with matplotlib and sci-kit learn
Jupyter Notebooks
Analysis on Genes with Significant P-Values
https://colab.research.google.com/drive/1LF5rxyu9mtnwrKj960KNgIYN73G7XlWm#scrollTo=ieM4tFHYaB4Q
Initial Statistical Analysis
https://colab.research.google.com/drive/11qf06KyouQ5Vh2B2FTH3PfjCs3QxFPWY#scrollTo=ieM4tFHYaB4Q
Analysis on Data Samples
https://colab.research.google.com/drive/1nd5pGaD-1iuYJ-yw_qLTS4uwxfcdW4fs#scrollTo=JlVlEGz87BVD
Pathway analysis
The pathway analysis for upregulated genes for patient pre versus on data pointed to the cAMP signaling pathway. (Table) The pathway analysis results for the upregulated genes agree with the published literature for this data, wherein they discussed the upregulation of immune checkpoint genes upon conducting transcriptome analysis (Riaz et al., 2017),(Sasi et al., 2021).
The cAMP pathway controls pro- and anti-inflammatory responses from the immune system (Raker et al., 2016). Immune cells produce anti-inflammatory chemicals when this pathway is active at high levels (thereby suppressing the innate and adaptive immune systems) (Sasi et al., 2021). For the downregulated genes, the analysis generated the complement and coagulation cascades pathway.
Figure 12: cAMP signaling pathway (taken from KEGG pathway)
Figure 13: Complement and coagulation cascades pathway (taken from KEGG pathway)
Table 2: Pathway analysis for upregulated and downregulated genes (taken from DAVID)
Interestingly, olfactory transduction and aldosterone-regulated sodium reabsorption were also a part of the pathways associated with our upregulated genes. The literature that published the dataset we used did not describe these pathways. On conducting a literature search, we saw that, specifically, the OR2C3 gene was suggested to have a functional part in melanoma progression (Ranzani et al., 2017). The olfactory pathway has also been enriched for other diseases that involved Nivolumab as a part of their treatment (Hodkinson et al., 2021).
Another literature survey was conducted to investigate the correlation between the aldosterone-regulated sodium reabsorption pathway and our study. The results of this survey pointed to a paper that described the side effects of using checkpoint inhibitor therapy (Trainer et al., 2016). Researchers reported nivolumab-induced adrenalitis that subsequently led to adrenal failure in a patient suffering from malignant melanoma (Trainer et al., 2016). Hence the aldosterone-regulated sodium reabsorption pathway could indicate potential immune-related adverse events.
Discussion and Literature Comparison
The literature published for the dataset used, analyzed genes that, regardless of response, displayed a large change in gene expression on therapy to see a pharmacological reaction for Nivolumab. For the purpose of our study, we examined all the genes irrespective of the degree of change in gene expression (Riaz et al., 2017).
We identified the cAMP pathway to be associated with upregulated genes. The cAMP cascade leads to the transcription of immune checkpoint PD-L1 and protein expression via the PKA enzyme (Sasi et al., 2021). Hence, comparing the pathways derived from the upregulated genes from our study with the existing literature, the general results converged at pathways associated with immune regulation.
For the downregulated genes, the pathways identified by the original paper were translation, cell-cycle regulation, melanin pathways, and mitotic division (Riaz et al., 2017). Our research pointed to the complement and coagulation cascades pathway for downregulated genes. Immune checkpoint inhibitors (like PD-L1) have an inhibitory relationship with the coagulation cascade, as they can aggravate inflammation (Joseph et al., 2020). PD-L1 causes coagulopathy by interrupting the inflammation-coagulation pathways (Joseph et al., 2020). Thereby these findings agree with the downregulation of complement and coagulation cascade-associated genes.
Relevant DEGs genes highlighted in the original literature parallel our derived list of DEGs. Common checkpoint genes CD274 (PD-L1), TNFRSF9 (4-1BB), PDCD1 (PD-1) LAG3, and CD80 (CTLA-4-L) were described to have higher expression irrespective of response to therapy in the literature source (Riaz et al., 2017). Immune-related genes sourced from pre-patient data were associated with cytokine signaling, lymphocyte activation, immune cytolytic activity, and chemotaxis. The genes found in our dataset of DEGs include HAVCR2 (TIM-3), TNFRSF4 (OX40), and TIGIT (Riaz et al., 2017).
In our initial hypothesis, we aimed to see the efficacy of the given treatment and identify the genes that are responsible for doing so. On comparing data derived from the original literature (Table S6) with our hierarchical clustering, we can see that there is a correlation between variance and mortality (Riaz et al., 2017). Our results point to the notion that higher genomic variance is indicative of increased chances of patient survival.
Our study also highlighted two previously non-reported pathways associated with this dataset that potentially connect to two main findings: 1. The OR2C3 has a functional role in melanoma (for the olfactory transduction pathway) (Ranzani et al., 2017), and 2. The drug used for the study, Nivolumab, may have detrimental side effects in treating patients with melanoma (Trainer et al., 2016).
Citations
Korman, A. J., Garrett-Thomson, S. C., & Lonberg, N. (2021). The foundations of immune checkpoint blockade and the ipilimumab approval decennial. Nature Reviews Drug Discovery 2021 21:7, 21(7), 509–528. https://doi.org/10.1038/s41573-021-00345-8
Seymour, L., Bogaerts, J., Perrone, A., Ford, R., Schwartz, L. H., Mandrekar, S., Lin, N. U., Litière, S., Dancey, J., Chen, A., Hodi, F. S., Therasse, P., Hoekstra, O. S., Shankar, L. K., Wolchok, J. D., Ballinger, M., Caramella, C., & de Vries, E. G. E. (2017). iRECIST: guidelines for response criteria for use in trials testing immunotherapeutics. The Lancet. Oncology, 18(3), e143–e152. https://doi.org/10.1016/S1470-2045(17)30074-8
Riaz, N., Havel, J. J., Makarov, V., Desrichard, A., Urba, W. J., Sims, J. S., Hodi, F. S., Martín-Algarra, S., Mandal, R., Sharfman, W. H., Bhatia, S., Hwu, W. J., Gajewski, T. F., Slingluff, C. L., Chowell, D., Kendall, S. M., Chang, H., Shah, R., Kuo, F., … Chan, T. A. (2017). Tumor and Microenvironment Evolution during Immunotherapy with Nivolumab. Cell, 171(4), 934-949.e15. https://doi.org/10.1016/J.CELL.2017.09.028
Raker, V. K., Becker, C., & Steinbrink, K. (2016). The cAMP pathway as therapeutic target in autoimmune and inflammatory diseases. Frontiers in Immunology, 7(MAR), 123. https://doi.org/10.3389/FIMMU.2016.00123/BIBTEX
Sasi, B., Ethiraj, P., Myers, J., Lin, A. P., Jiang, S., Qiu, Z., Holder, K. N., & Aguiar, R. C. T. (2021). Regulation of PD-L1 expression is a novel facet of cyclic-AMP-mediated immunosuppression. Leukemia, 35(7), 1990–2001. https://doi.org/10.1038/S41375-020-01105-0
Joseph, J. J., Rajan, A., Gulley, J. L., Ito, S., & Kessler, C. M. (2020). Acquired Coagulopathy With Immune Checkpoint Inhibitors: An Underrecognized Association Between Inflammation and Coagulation. JTO Clinical and Research Reports, 1(3). https://doi.org/10.1016/J.JTOCRR.2020.100049
Ranzani, M., Iyer, V., Ibarra-Soria, X., Del Castillo Velasco-Herrera, M., Garnett, M., Logan, D., & Adams, D. J. (2017). Revisiting olfactory receptors as putative drivers of cancer. Wellcome Open Research, 2. https://doi.org/10.12688/WELLCOMEOPENRES.10646.1Hodkinson, B. P., Schaffer, M., Brody, J. D., Jurczak, W., Carpio, C., Ben-Yehuda, D., Avivi, I., Forslund, A., Özcan, M., Alvarez, J., Ceulemans, R., Fourneau, N., Younes, A., & Balasubramanian, S. (2021). Biomarkers of response to ibrutinib plus nivolumab in relapsed diffuse large B-cell lymphoma, follicular lymphoma, or Richter’s transformation. Translational Oncology, 14(1), 100977. https://doi.org/10.1016/J.TRANON.2020.100977l
Raieli, S. (2021, December 13). Clustering techniques with gene expression data. Retrieved November 7, 2022, from https://medium.com/leukemiaairesearch/clustering-techniques-with-gene-expression-data-4b35a04f87d5
Mikulskibartosz, B. (2019, June 03). PCA - how to choose the number of components? Retrieved November 7, 2022, from https://www.mikulskibartosz.name/pca-how-to-choose-the-number-of-components/
Pedregosa, F., Varoquaux, Ga”el, Gramfort, A., Michel, V., Thirion, B., Grisel, O., … others. (2011). Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 12(Oct), 2825–2830.
McKinney, W., & others. (2010). Data structures for statistical computing in python. In Proceedings of the 9th Python in Science Conference (Vol. 445, pp. 51–56).
Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., … SciPy 1.0 Contributors. (2020). SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods, 17, 261–272. https://doi.org/10.1038/s41592-019-0686-2
Waskom, M., Botvinnik, Olga, O’Kane, Drew, Hobson, Paul, Lukauskas, Saulius, Gemperline, David C, … Qalieh, Adel. (2017). mwaskom/seaborn: v0.8.1 (September 2017). Zenodo. https://doi.org/10.5281/zenodo.883859
Kluyver, T., Ragan-Kelley, B., Fernando P’erez, Granger, B., Bussonnier, M., Frederic, J., … Willing, C. (2016). Jupyter Notebooks – a publishing format for reproducible computational workflows. In F. Loizides & B. Schmidt (Eds.), Positioning and Power in Academic Publishing: Players, Agents and Agendas (pp. 87–90).
Galaxy. (n.d.). Retrieved November 7, 2022, from https://usegalaxy.org/
David functional annotation bioinformatics Microarray analysis. (n.d.). Retrieved November 7, 2022, from https://david.ncifcrf.gov/
Comparison of Gene Expression Changes in Response to Haloperidol Treatment in Brain and Liver Mice Samples
Maya Jaffe, Savannah Richardson, Yuzheng (Leo) Yang, Travis Land, Alex Kovensky
INTRODUCTION
For our project, we decided to compare the effect of Haloperidol in the brain and liver of mice while using placebo mice as a control. The study from which our datasets were taken focused primarily on understanding how the striatum in the brain was impacted by chronic Haloperidol dosage by utilizing RNA-seq while also generating gene expression microarrays of the whole brain and liver to further understand the overall effects of Haloperidol [1]. We focused on comparing the control livers’ and brains’ microarrays to the treated tissues to analyze the gene expression that occurs under this chronic dosage. Understanding the effects of this treatment outside of the target region will provide further information into the benefits and adverse effects of Haloperidol.
BACKGROUND
Genome-wide association studies (GWAS) for schizophrenia have identified around 500 encoding genes between 100 loci; however, except for the D2 receptor, it is unclear whether these genes are responsible for antipsychotic effects or if they are antipsychotic drug targets [1]. Haloperidol is a conventional antipsychotic commonly used in treating and managing positive symptoms of schizophrenia, including hallucinations and delusions [2]. Whilst the brain has a high density of dopamine receptors, which serve as the target for haloperidol, the liver was used to look at effects of the drug outside of the brain and to see if there are any potential gene targets present. Although uncommon, haloperidol can have adverse effects on the liver and has been linked to clinically apparent acute liver injury [3]. While haloperidol can have benefits for the psychiatric health and brain of treated individuals, the associated adverse effects on the liver make understanding the tradeoffs of this treatment important.
RESEARCH GOAL
We analyzed the gene expression changes in mice brain and liver samples after being treated with haloperidol. To understand the effects the drug has on relative tissues, our goal was to identify differentially expressed genes and analyze any similarities and differences between the two sets. We also wanted to understand what pathways were involved and how those compared between the various tissue samples.
METHODS
DATA SETS
The study from which our data is sourced looked at the effect of haloperidol, an antipsychotic used to treat schizophrenia, on mouse brain and liver tissue [1]. Mice were implanted with slow releasing haloperidol tablets, which produce about the same levels of drug in plasma and brain similar to humans. After 30 days, the mice were killed by cervical dislocation to reduce effects on gene expression and tissue was subsequently frozen. 38 mice in total were used for the study. 20 mice had brain tissue harvested with half being treated with the drug and half being placebo. Liver tissue was harvested from 18 mice, with half treated with haloperidol and half being a placebo. The platform used is GPL11533 ([MoGene-1_1-st] Affymetrix Mouse Gene 1.1 ST Array [transcript (gene) version]) [4].
For our first dataset, we used 1 brain placebo and 1 treated brain sample for a total of 8 samples. For our second dataset, we used 1 liver placebo and 1 treated liver sample for a total of 8 samples (Table 1).
BIOINFORMATICS PIPELINE
DATA AND STATISTICAL ANALYSIS
The statistical method that we used to analyze our data is a correlation test with Pearson algorithm and a fold-change threshold value of 1. This correlation method allowed us to compare the similarities and differences between our brain and liver datasets. We decided to use a fold change of 1, which is typically used when analyzing RNA-seq data; any higher fold-change threshold was unable to properly generate anything in ExAtlas for our datasets [5].
For our heatmap and matrices, we decided to use a fold change threshold of 2 and a PC clustering value of 0.7. Once again, these values are typical for analyzing RNA-seq data while also generating the best results in ExAtlas. In order to control for FDRs, we used the standard value of 0.05 [6].
In performing hierarchical clustering, we used the function pheatmap in R and clustered via Euclidean distances. The clustering method used was a “complete linkage” method. Complete linkage clustering is an agglomerative method, where each element starts in its own cluster and the maximum distance between clusters is found. The similarities of two clusters are determined by the similarities of the clusters’ most dissimilar members. The clustering trees were then cut based on the largest clusters to see similarities within gene expression between samples [7-9].
Following differential gene expression, principal component analysis (PCA) was performed. PCA reduces the dimensionality of the data set by standardizing the data, calculating a covariance matrix, then finding the eigenvalues and eigenvectors of the principal components. The samples are plotted based on the principal components that have the highest combined variation [10].
RESULTS
In the PCA plot for the brain dataset, the experimental samples (drug) and control samples (placebo) are well-clustered, indicating a difference in gene expression between the experimental and control groups. Alternatively, clustering was poor between the experimental and control samples in the liver dataset, suggesting that drug treatment does not have much notable effect on gene expression in liver cells.
The pathway analysis is performed by Inputting our differentially expressed gene matrices into the Reactome. The reactome pathway database and the IntAct molecular interactions database were used in the pathway analysis. The PADOG method is utilized for analysis, which is a weighted gene set analysis method that down-weighs genes that are present in many pathways. The pathway analysis allowed for identification of over-represented pathways within both the brain and liver datasets (Table 2, Figure 2), as well as pathway overlap between datasets [11-18].
While some pathways were identified as overlapping between the brain and liver datasets, all of these pathways have a p-value and FDR > 0.05 in at least one dataset. These overlapping pathways include:
- Autophagy (R-MMU-9612973)
- Selective autophagy (R-MMU-9663891)
- Macroautophagy (R-MMU-1632852)
- mRNA splicing- major pathway (R-MMU-72163)
- Processing of Capped Intron-Containing Pre-mRNA (R-MMU-72203)
- mRNA splicing (R-MMU-72172)
- Cytokine Signaling in the Immune system (R-MMU-1280215)
- Immune System (R-MMU-168256)
- Metabolism (R-MMU-1430728)
- Metabolism of RNA (R-MMU-8953854)
- Metabolism of proteins (R-MMU-392499)
- Citric acid (TCA) cycle and respiratory electron transport (R-MMU-1428517)
- Prolactin receptor signaling pathway (R-MMU-1170546)
- Growth hormone receptor signaling pathway (R-MMU-982772)
DISCUSSION & LITERATURE COMPARISON
Based on the heatmaps and PCA plots generated, haloperidol has a clear effect on brain cell gene expression, but no notable effect in liver cell gene expression. Schizophrenia shows a different gene expression profile, specifically related to brain aging, in comparison to normal patients [19]. Based on haloperidol’s usage as an antipsychotic for treatment of positive symptoms of schizophrenia, the drug having an effect on brain cell gene expression is not particularly surprising [2].
Relatedly, our pathway analysis results through reactome indicate that there is not much overlap between over-represented pathways in brain and liver cells, as shown by the high P value and FDR for said overlapped pathways. This is also not surprising since the liver and brain cells are highly differentiated cells with almost completely different functions.
Overall our results indicate that haloperidol mainly has effects on brain cells and its effect on liver cells is not detectable. This finding is further supported by later studies such as that from Dean and Scarr, which concluded that antipsychotic drugs like haloperidol and chlorpromazine use changes in gene expression as a method for reducing psychotic symptoms after identifying significant fold change in over 150 genes located in the brain [20]. Besides blocking some of the dopamine receptors in brain cells, long-term antipsychotic treatment is also thought to be associated with brain volume reduction [21]. In Dr. Beng-Choon Ho’s study, they investigated the data from two hundred and eleven patients with schizophrenia who underwent repeated neuroimaging – beginning soon after illness onset – in an Iowa longitudinal study. Ho et al. found that the amount of antipsychotic treatment is negatively related to the volume of gray matter in the patient’s brain. Although the mechanism of how the antipsychotic causes reduction of brain volume is unclear, this question could be addressed in future research. Furthermore, later studies could relate the gene expression profiles in the brain to phenotypic results using weighted gene co-expression network analysis (WGCNA), as recent studies show haloperidol can cause tardive dyskinesia [22]. This would allow for a deeper dive into how specifically the overrepresented pathways are expressed phenotypically, leading to possible options for treatment or symptom reduction.
REFERENCES
- Kim, Y., Giusti-Rodriguez, P., Crowley, J. et al. Comparative genomic evidence for the involvement of schizophrenia risk genes in antipsychotic effects. Mol Psychiatry 23, 708–712 (2018). https://doi.org/10.1038/mp.2017.111
- Rahman S, Marwaha R. Haloperidol. [Updated 2022 Jul 4]. In: StatPearls [Internet]. Treasure Island (FL): StatPearls Publishing; 2022 Jan-. Available from: https://www.ncbi.nlm.nih.gov/books/NBK560892/
- LiverTox: Clinical and Research Information on Drug-Induced Liver Injury [Internet]. Bethesda (MD): National Institute of Diabetes and Digestive and Kidney Diseases; 2012-. Haloperidol. [Updated 2018 Mar 25]. Available from: https://www.ncbi.nlm.nih.gov/books/NBK548393/
- Alternative splicing regulates mouse embryonic stem cell pluripotency and differentiation. Salomonis N, Schlieve CR, Pereira L, Wahlquist C, Colas A, Zambon AC, Vranizan K, Spindler MJ, Pico AR, Cline MS, Clark TA, Williams A, Blume JE, Samal E, Mercola M, Merrill BJ, Conklin BR. Proc Natl Acad Sci U S A. 2010 Jun 8;107(23):10514-9. Epub 2010 May 24.
- Matthew E. Ritchie, Belinda Phipson, Di Wu, Yifang Hu, Charity W. Law, Wei Shi, Gordon K. Smyth, limma powers differential expression analyses for RNA-sequencing and microarray studies, Nucleic Acids Research, Volume 43, Issue 7, 20 April 2015, Page e47, https://doi.org/10.1093/nar/gkv007
- J Bioinform Comput Biol . 2015 Dec;13(6):1550019. doi: 10.1142/S0219720015500195. Epub 2015 Jun 9. (ExAtlas)
- Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, Feng T, Zhou L, Tang W, Zhan L, Fu x, Liu S, Bo X, Yu G (2021). “clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.” The Innovation, 2(3), 100141. doi: 10.1016/j.xinn.2021.100141.
- Yu G, Wang L, Han Y, He Q (2012). “clusterProfiler: an R package for comparing biological themes among gene clusters.” OMICS: A Journal of Integrative Biology, 16(5), 284-287. doi: 10.1089/omi.2011.0118.
- AltAnalyze and DomainGraph: analyzing and visualizing exon expression data. Emig D, Salomonis N, Baumbach J, Lengauer T, Conklin BR, Albrecht M. Nucleic Acids Res. 2010 Jul 1;38 Suppl:W755-62. Epub 2010 May 31. co-first author
- Single-cell analysis of mixed-lineage states leading to a binary cell fate choice. Olsson A, Venkatasubramanian M, Chaudhri VK, Aronow BJ, Salomonis N, Singh H, Grimes HL. Nature. 2016 Aug 31;537(7622):698-702. co-corresponding*
- Marc Gillespie, Bijay Jassal, Ralf Stephan, Marija Milacic, Karen Rothfels, Andrea Senff-Ribeiro, Johannes Griss, Cristoffer Sevilla, Lisa Matthews, Chuqiao Gong, Chuan Deng, Thawfeek Varusai, Eliot Ragueneau, Yusra Haider, Bruce May, Veronica Shamovsky, Joel Weiser, Timothy Brunson, Nasim Sanati, Liam Beckman, Xiang Shao, Antonio Fabregat, Konstantinos Sidiropoulos, Julieth Murillo, Guilherme Viteri, Justin Cook, Solomon Shorser, Gary Bader, Emek Demir, Chris Sander, Robin Haw, Guanming Wu, Lincoln Stein, Henning Hermjakob, Peter D’Eustachio, The reactome pathway knowledgebase 2022, Nucleic Acids Research, 2021;, gkab1028, https://doi.org/10.1093/nar/gkab1028
- Griss J, Viteri G, Sidiropoulos K, Nguyen V, Fabregat A, Hermjakob H. ReactomeGSA – Efficient Multi-Omics Comparative Pathway Analysis. Mol Cell Proteomics. 2020 Sep 9. doi: 10.1074/mcp. PubMed PMID: 32907876.
- Jassal B, Matthews L, Viteri G, Gong C, Lorente P, Fabregat A, Sidiropoulos K, Cook J, Gillespie M, Haw R, Loney F, May B, Milacic M, Rothfels K, Sevilla C, Shamovsky V, Shorser S, Varusai T, Weiser J, Wu G, Stein L, Hermjakob H, D’Eustachio P. The reactome pathway knowledgebase. Nucleic Acids Res. 2020 Jan 8;48(D1):D498-D503. doi: 10.1093/nar/gkz1031. PubMed PMID: 31691815.
- Fabregat A, Korninger F, Viteri G, Sidiropoulos K, Marin-Garcia P, Ping P, Wu G, Stein L, D’Eustachio P, Hermjakob H. Reactome graph database: Efficient access to complex pathway data. PLoS Comput Biol. 2018 Jan 29;14(1):e1005968. doi: 10.1371/journal.pcbi.1005968. eCollection 2018 Jan. PubMed PMID: 29377902.
- Fabregat A, Sidiropoulos K, Viteri G, Marin-Garcia P, Ping P, Stein L, D’Eustachio P, Hermjakob H. Reactome diagram viewer: data structures and strategies to boost performance. Bioinformatics. 2018 Apr 1;34(7):1208-1214. doi: 10.1093/bioinformatics/btx752. PubMed PMID: 29186351.
- Sidiropoulos K, Viteri G, Sevilla C, Jupe S, Webber M, Orlic-Milacic M, Jassal B, May B, Shamovsky V, Duenas C, Rothfels K, Matthews L, Song H, Stein L, Haw R, D’Eustachio P, Ping P, Hermjakob H, Fabregat A. Reactome enhanced pathway visualization. Bioinformatics. 2017 Nov 1;33(21):3461-3467. doi: 10.1093/bioinformatics/btx441. PubMed PMID: 29077811.
- Fabregat A, Sidiropoulos K, Viteri G, Forner O, Marin-Garcia P, Arnau V, D’Eustachio P, Stein L, Hermjakob H. Reactome pathway analysis: a high-performance in-memory approach. BMC Bioinformatics. 2017 Mar 2;18(1):142. doi: 10.1186/s12859-017-1559-2. PubMed PMID: 28249561.
- Wu G, Haw R. Functional Interaction Network Construction and Analysis for Disease Discovery. Methods Mol Biol. 2017;1558:235-253. doi: 10.1007/978-1-4939-6783-4_11. PubMed PMID: 28150241.
- Sabunciyan, S. Gene Expression Profiles Associated with Brain Aging are Altered in Schizophrenia. Sci Rep 9, 5896 (2019). https://doi.org/10.1038/s41598-019-42308-5
- Dean, B., Scarr, E. Common changes in rat cortical gene expression after chronic treatment with chlorpromazine and haloperidol may be related to their antipsychotic efficacy. Neuroscience Applied, (2022). https://doi.org/10.1016/j.nsa.2022.101015
- Ho, B. C., Andreasen, N. C., Ziebell, S., Pierson, R., & Magnotta, V. (2011). Long-term antipsychotic treatment and brain volumes: a longitudinal study of first-episode schizophrenia. Archives of general psychiatry, 68(2), 128–137. https://doi.org/10.1001/archgenpsychiatry.2010.199
- Paola Giusti-Rodríguez, James G Xenakis, James J Crowley, Randal J Nonneman, Daniela M DeCristo, Allison Ryan, Corey R Quackenbush, Darla R Miller, Ginger D Shaw, Vasyl Zhabotynsky, Patrick F Sullivan, Fernando Pardo Manuel de Villena, Fei Zou, Antipsychotic Behavioral Phenotypes in the Mouse Collaborative Cross Recombinant Inbred Inter-Crosses (RIX), G3 Genes|Genomes|Genetics, Volume 10, Issue 9, 1 September 2020, Pages 3165–3177, https://doi.org/10.1534/g3.120.400975
Identifying drivers of metastasis in small-cell and non-small cell lung cancer
Joel Markus Vaz, Harini Mudradi, Asmita Lagwankar, Shreya Rajasekar, Siddhi Ranbhor
INTRODUCTION AND BACKGROUND
Lung cancer is the second most diagnosed cancer in both men and women in the United States. Lung cancer is the most common cancer worldwide, accounting for 2.1 million new cases and 1.8 million deaths yearly. Lung cancer (LC) care cost accounts for $13.4 billion in the U.S. (Miller, 2005). It is a molecularly heterogeneous disease that can mainly be divided into two major subtypes: Small Cell Lung Cancer (SCLC) and Non-Small Cell Lung Cancer (NSCLC). NSCLC accounts for most lung cancers but has a higher survival rate than SCLC. The survival rate for NSCLC is about 20% higher than that for SCLC. Early detection of lung carcinoma is challenging because of the lack of specific symptoms and rapid tumor growth, making it necessary to identify the drivers of metastasis.
SCLC is a highly malignant tumor derived from cells exhibiting neuroendocrine characteristics. SCLC often starts in the bronchi or the airways that lead from the trachea into the lungs and then branches off into progressively smaller structures. After affecting the bronchi, SCLC quickly grows and spreads to other body parts, including the lymph nodes. This type is typically caused by tobacco smoking. NSCLC is divided into three subtypes: adenocarcinoma, squamous cell carcinoma, and large cell lung carcinoma, based on the area of lungs affected. The gene mutations and pathogenesis of cancer depend on the type. 8% of cases are due to inherited genetic factors. (D’Angelo & Pietanza, 2010) Lung epithelium dysplasia occurs after repeated exposure to carcinogens, specifically cigarette smoke. If exposure persists, genetic mutations are caused, and protein synthesis occurs. EGFR, c-MET, NKX2-1, LKB1, PIK3CA and BRAF are some of the genes that play a role in pathogenesis. MYC, BCL2 and p53 are responsible for SCLC’s most common genetic mutations. Similarly, EGFR, KRAS, P10, TP53, and CDK4 are the most commonly altered genes in the oncogenic pathways in NSCLC (dela Cruz et al., 2011).
The aim of our project is to compare primary and metastasized small cell lung cancer against non-small cell lung cancer to understand the underlying differences and visualize the phenotypic alterations between the two subtypes of lung cancer with the help of PCA, hierarchical plots and pathway analysis. We have used transcriptome data to analyze the oncogenic pathway at a broader level than just the genetic alterations.
RESEARCH GOALS
We undertook an analysis that relates common and unique mechanisms of invasiveness in SCLC and squamous cell NSCLC at a transcriptome level. With this analysis, we set out to see whether-
- Both SCLC and NSCLC exhibit both typical and distinctive metastasis-related indicators.
- Compared to NSCLC, SCLC has a greater risk of aggressiveness.
- Compared to NSCLC, SCLC has an extremely poor prognosis and outlook.
METHODS
Datasets
The patient’s raw RNA-seq data were obtained from the Cancer Cell Line Encyclopedia (CCLE) dataset (Ghandi et al., 2019). 28 cell lines were chosen, seven each from primary and metastasis SCLC and NSCLC. These samples were paired-end reads done on a non-strand-specific assay protocol. The data was collected from NCBI BioProject PRJNA523380 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA523380/)
Data pre-processing
The analysis workflow is defined in Fig. 1. Quality control checks were initially done on the samples using FastQC to assess their PHRED scores and adaptor sequences. The sample reads were aligned to the hg38 reference genome using Bowtie2 (Langmead & Salzberg, 2012). Maximum fragment length was considered as 500, with 0 mismatches as a limit for seed alignment. The resultant BAM files were sorted and indexed using SAMtools (Danecek et al., 2021). HTSeq was used to generate raw count matrices from the aligned read files with a MAPQ quality cutoff of 10 (Anders et al., 2014). The log2-normalized Transcripts Per Kilobase Million mapped reads (TPM) were generated from the count matrix using an in-house R script.
Differential gene expression
Differential expression was performed using the R package, ‘DESeq2’ (Love et al., 2014).DESeq2 uses negative binomial distribution to model the counts into a generalized linear model (GLM). Wald test was applied to test for the hypothesis. For differential expression, we used a log2(Fold Change) (log2FC) cutoff of 1 and a p-value cutoff of 0.05.
Hierarchical clustering
Hierarchical clustering of the genes and heatmaps was performed using the ‘pheatmap’ package in R. For hierarchical clustering, only the genes that were differentially expressed in at least one comparison were chosen. The TPM values were z-score normalized across genes before clustering.
PCA clustering
Higher dimensional clustering was performed using the R package ‘Seurat’ (Hao et al., 2021). TPM normalization of the samples was used instead of the standard Seurat normalization, the variable features were identified using the ‘vst’ method, and the data was scaled before running for PCA using only the top 2000 variable genes.
Pathway analysis
Three pathway analysis tools and databases were used to determine pathways significantly across comparisons. The graphical web tool, ShinyGO, was utilized to identify enriched pathways in the KEGG databases (Ge et al., 2020), while the PANTHER knowledge base was implemented for GO pathways (Thomas et al., 2022). For the enrichment of cancer Hallmark pathways, single-sample Gene Set Enrichment Analysis (ssGSEA) was carried out for each sample using the Enrichr Python wrapper, GSEApy (Subramanian et al., 2007).
Fig. 1: Proposed Workflow
RESULTS
Our study aimed to identify common and individual mechanisms of metastasis between SCLC and NSCLC, more specifically, squamous cell carcinoma. We set out to use 7 samples each for NSCLC and SCLC metastasis and primary tumor data. Initial multi-dimensional clustering analysis revealed an outlier in NSCLC metastasis, forcing us to remove it for our downstream analysis. We performed differential expression first to identify the set of significantly up or downregulated genes across three comparisons: Metastasis vs primary NSCLC and SCLC, and all samples across NSCLC vs SCLC (Fig. 2-4). Our findings established that 218 and 177 genes were upregulated in SCLC and NSCLC, respectively. Out of these, 12 were commonly upregulated and enriched in metastatic samples across the two cancer types (Fig. 5). Seven of such genes were observed to be downregulated (Fig. 6). Regarding common genes between the primary/metastasis comparisons and NSCLC vs SCLC, we observed that there were more common genes when cross compared than when compared between the same cancer type (38 + 15 = 53 cross-comparison over 12 between cancer types in upregulated, and 42 over 7 in downregulated). This would suggest more similarities between the cancer subtypes than primary and metastasis samples within the cancer type itself, meaning tumor cells tend to diverge more from their characteristics as they progress from primary to metastatic cells than how they were collected as one cancer type.
Out of the top 5 differentially expressed genes across the three comparisons (Fig. 7), it is interesting to note that gene BRDTP1 was significantly upregulated in SCLC while downregulated in NSCLC, implying its shift in role from one cancer type to another, if at all it is known to be associated in lung cancer.
Fig. 2: Volcano plot for Metastasis vs. Primary SCLC samples. Red points display upregulated genes, and blue displays down-regulated ones.
Fig. 3: Volcano plot for Metastasis vs. Primary NSCLC samples. Red points display upregulated genes, and blue displays down-regulated ones.
Fig. 4: Volcano plot for Squamous cell carcinoma vs. SCLC samples. Red points display upregulated genes, and blue displays downregulated ones.
Fig. 5: Venn diagram for downregulated genes across three comparisons
Fig. 6: Venn diagram for upregulated genes across three comparisons
Fig. 7: Top 5 upregulated and downregulated genes for the three comparisons.
We next performed multi-dimensional clustering using Seurat for the 27 samples. NSCLC and SCLC were segregated well on the PC1 axis, barring a few points. Among the metastasis vs primary comparisons, we observe a greater scattering of points, which only validates our previous hypothesis that primary and metastasis samples don’t cluster as well as individual cancer subtypes.
Fig. 8: PCA of all samples using top 2 PCs and top 2000 variable features.
Given that our results in PCA and differential gene expression directed towards significant differences between primary and metastasis mechanisms between SCLC and NSCLC and transcriptome profiles not leading to clusters, we attempted to profile these samples using hierarchical clustering. We used TPM normalized data, which was further z-score normalized across the genes as an input for hierarchical clustering. Unsurprisingly, NSCLC and SCLC samples clustered together when differential genes were used, while primary and metastasis samples of both the cancer types were scattered throughout. While this might further match our previous findings, we should also take into account that there was a ten-fold increase in the differentially expressed genes while cross-comparing samples (~2300 genes) as compared to differential expression within the cancer type (~200 in each), which might have skewed our analysis in favor of better clustering for the cross-comparison. It is probable that intra-cancer sample comparison, without using inter-cancer DE genes, could aid in a more precise clustering.
Fig. 9: Hierarchical clustering of differentially expressed genes for metastasis and primary samples in SCLC and NSCLC.
Pathway analysis was performed using the Panther database and ShinyGo. We used Panther to generate tables on pathway analysis, and ShinyGo was used for the visualization of those pathway analysis. The most frequent pathways for each comparison is tabulated in Table 1.
Table 1: Table of selected pathways generated from PANTHER for biological, cellular and molecular processes.
| Biological processes | Cellular processes | Molecular function | |
| Pri SCLC Upregulated | response to stimulus (GO:0050896) | plasma membrane (GO:0005886) | molecular_function (GO:0003674) |
| Pri SCLC Downregulated | negative regulation of chronic inflammatory response (GO:0002677) | extrinsic component of postsynaptic specialization membrane (GO:0098892) | cytokine activity (GO:0005125) |
| Met_pri_NSCLC_DEs Upregulated | epithelial cell differentiation (GO:0030855) | – | – |
| Met_pri_NSCLC_DEs Downregulated | respiratory burst involved in defence response (GO:0002679) | midbody (GO:0030496) | protein serine/threonine kinase activity (GO:0004674) |
| NSCLCvsSCLC Upregulated | cell communication (GO:0007154) | cellular anatomical entity (GO:0110165) | protein-containing complex binding (GO:0044877) |
| NSCLCvsSCLC Downregulated | transmembrane transport (GO:0055085) | membrane (GO:0016020) | nucleic acid binding (GO:0003676) |
Fig. 10 refers to cAMP Signaling Pathway involving genes LH, FSH, TSH, SST, NMDAR, VAV2. To diminish SIRT6 expression by boosting its ubiquitin-proteasome-dependent degradation and therefore promote radiation-induced apoptosis in lung cancer cells, the cAMP signaling pathway can mediate the PKA-dependent suppression of the Raf-MEK-ERK pathway. (T.-P. Lu et al., 2010)
Fig. 10: cAMP Signaling Pathway.
Cytopenias Associated with Cancer which involves HLA-DR, CD44, CD127, IL-1, and CD36 genes are shown in Fig. 11. By displacing and killing stem and progenitor cells, harming the bone marrow microenvironment, hindering the generation of hematopoietic growth factors, or triggering the production of cytokines that inhibit hematopoiesis, metastatic disease in the bone marrow can disturb hematopoiesis. (Zuckerman, 1998).
Fig. 11: Hematopoietic cell lineage.
Fig. 12 shows ECM Receptor Interaction involving genes Tenascin HLA-DR, CD44, CD127, IL-1, CD36, CD44, CD36. Abnormal ECM impacts how cancer progresses by directly encouraging cellular transformation and metastasis. Importantly, however, stromal cell activity is also deregulated by ECM abnormalities. This promotes tumor-associated angiogenesis and inflammation, which creates a tumorigenic surrounding (P. Lu et al., 2012).
Fig. 12: ECM Receptor Interaction.
Fig. 13 shows the neuroactive ligand-receptor interaction involving genes ADR, KISS1, SST, FSH, LHB, GRI, GHR. The neuroactive ligand-receptor interaction pathway is one of the top 10 most enriched pathways due to the downregulation of several GPCRs and the upregulation of a handful of others including Kissr1, Mchr1, Ptger3, Lpar2, and Gabbr1. We can also infer the gene alterations from the pathway analysis figure above. ADR, KISS1, SST, PRLHR, FSH, LHB, TSH, GRIN, GRI, and GHR are getting altered. In tissues, this pathway is also enriched, mainly due to the downregulation of several of its DEGs (Nguyen et al., 2021).
Fig. 13: Neuroactive Ligand Receptor Interaction Pathway
Fig. 14 shows the enrichment analysis and protein-protein network of the common differentially expressed genes using the String tool. Functional enrichments and various cellular processes can be analyzed and derived from the network. Each node signifies a gene and returns its functions. String uses information from different sources like genomic context predictions, high throughput experiments, co-expression and automated text mining.
Fig. 14: Network depicting the protein-protein interactions of differentially expressed genes
Finally, Hallmark pathways were analyzed using ssGSEA (Fig. 15). Single-sample GSEA facilitates enrichment scores for each sample, resulting in an absolute comparison of scores across each sample. The normalized enrichment scores were utilized for the analysis. It was observed that most NSCLC samples were enriched for all Hallmark pathways compared to SCLC. For metastasis-related pathways like epithelial-mesenchymal transition (EMT) and TGF-β signaling, metastasis samples were upregulated, as expected.
Fig. 15: Heatmap of the normalized ssGSEA scores for MSigDB Hallmark gene sets.
DISCUSSION
In our study, we focused on decoding mechanisms of metastasis at a transcriptome level. The phenomenon of metastasis undergoes much quicker than the process of cell division. This makes it difficult to explain why tumor cells metastasize while having the same genetic makeup as normal cells. These tumor cells need not undergo genetic mutations to progress, and our analysis deduced phenotypic mechanisms of SCLC and NSCLC between and across these cancer types.
From the clustering analysis, we observed clustering between NSCLC and SCLC samples but not individual metastasis and primary samples. It was later established that there were more common upregulated/downregulated genes across the cancer subtypes than there were within. We inferred that NSCLC and SCLC might have similar profiles; however, profiling metastasis from primary, even in individual cancer subtypes, is difficult as clustering was poor and common pathways were fewer. This could be due to greater heterogeneity in primary and metastasis tumor samples. This heterogeneity is what leads to aggressiveness in cancer cells (Marino et al., 2019). Moreover, we initially asked if SCLC was more aggressive than NSCLC. Since the clustering was slightly worse in SCLC, there was a greater heterogeneity in the transcriptome profiling, which suggests greater aggressiveness in return. This could map to poorer prognosis and clinical outcome, as aggressive cells are more invasive in their host environment.
REFERENCES
Anders, S., Pyl, P. T., & Huber, W. (2014). HTSeq – A Python framework to work with high-throughput sequencing data. https://doi.org/10.1101/002824
Danecek, P., Bonfield, J. K., Liddle, J., Marshall, J., Ohan, V., Pollard, M. O., Whitwham, A., Keane, T., McCarthy, S. A., Davies, R. M., & Li, H. (2021). Twelve years of SAMtools and BCFtools. GigaScience, 10(2). https://doi.org/10.1093/gigascience/giab008
D’Angelo, S. P., & Pietanza, M. C. (2010). The molecular pathogenesis of small cell lung cancer. Cancer Biology & Therapy, 10(1), 1–10. https://doi.org/10.4161/cbt.10.1.12045
dela Cruz, C. S., Tanoue, L. T., & Matthay, R. A. (2011). Lung Cancer: Epidemiology, Etiology, and Prevention. Clinics in Chest Medicine, 32(4), 605–644. https://doi.org/10.1016/j.ccm.2011.09.001
Ge, S. X., Jung, D., & Yao, R. (2020). ShinyGO: a graphical gene-set enrichment tool for animals and plants. Bioinformatics, 36(8), 2628–2629. https://doi.org/10.1093/bioinformatics/btz931
Ghandi, M., Huang, F. W., Jané-Valbuena, J., Kryukov, G. v., Lo, C. C., McDonald, E. R., Barretina, J., Gelfand, E. T., Bielski, C. M., Li, H., Hu, K., Andreev-Drakhlin, A. Y., Kim, J., Hess, J. M., Haas, B. J., Aguet, F., Weir, B. A., Rothberg, M. v., Paolella, B. R., … Sellers, W. R. (2019). Next-generation characterization of the Cancer Cell Line Encyclopedia. Nature, 569(7757), 503–508. https://doi.org/10.1038/s41586-019-1186-3
Hao, Y., Hao, S., Andersen-Nissen, E., Mauck, W. M., Zheng, S., Butler, A., Lee, M. J., Wilk, A. J., Darby, C., Zager, M., Hoffman, P., Stoeckius, M., Papalexi, E., Mimitou, E. P., Jain, J., Srivastava, A., Stuart, T., Fleming, L. M., Yeung, B., … Satija, R. (2021). Integrated analysis of multimodal single-cell data. Cell, 184(13), 3573-3587.e29. https://doi.org/10.1016/j.cell.2021.04.048
Langmead, B., & Salzberg, S. L. (2012). Fast gapped-read alignment with Bowtie 2. Nature Methods, 9(4), 357–359. https://doi.org/10.1038/nmeth.1923
Love, M. I., Huber, W., & Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15(12), 550. https://doi.org/10.1186/s13059-014-0550-8
Lu, P., Weaver, V. M., & Werb, Z. (2012). The extracellular matrix: A dynamic niche in cancer progression. Journal of Cell Biology, 196(4), 395–406. https://doi.org/10.1083/jcb.201102147
Lu, T.-P., Tsai, M.-H., Lee, J.-M., Hsu, C.-P., Chen, P.-C., Lin, C.-W., Shih, J.-Y., Yang, P.-C., Hsiao, C. K., Lai, L.-C., & Chuang, E. Y. (2010). Identification of a Novel Biomarker, SEMA5A, for Non–Small Cell Lung Carcinoma in Nonsmoking Women. Cancer Epidemiology, Biomarkers & Prevention, 19(10), 2590–2597. https://doi.org/10.1158/1055-9965.EPI-10-0332
Marino, F. Z., Bianco, R., Accardo, M., Ronchi, A., Cozzolino, I., Morgillo, F., Rossi, G., & Franco, R. (2019). Molecular heterogeneity in lung cancer: from mechanisms of origin to clinical implications. International Journal of Medical Sciences, 16(7), 981–989. https://doi.org/10.7150/ijms.34739
Miller, Y. E. (2005). Pathogenesis of Lung Cancer. American Journal of Respiratory Cell and Molecular Biology, 33(3), 216–223. https://doi.org/10.1165/rcmb.2005-0158OE
Nguyen, X.-X., Renaud, L., & Feghali-Bostwick, C. (2021). Identification of Impacted Pathways and Transcriptomic Markers as Potential Mediators of Pulmonary Fibrosis in Transgenic Mice Expressing Human IGFBP5. International Journal of Molecular Sciences, 22(22), 12609. https://doi.org/10.3390/ijms222212609
Subramanian, A., Kuehn, H., Gould, J., Tamayo, P., & Mesirov, J. P. (2007). GSEA-P: a desktop application for Gene Set Enrichment Analysis. Bioinformatics, 23(23), 3251–3253. https://doi.org/10.1093/bioinformatics/btm369
Thomas, P. D., Ebert, D., Muruganujan, A., Mushayahama, T., Albou, L., & Mi, H. (2022). PANTHER: Making genome‐scale phylogenetics accessible to all. Protein Science, 31(1), 8–22. https://doi.org/10.1002/pro.4218
Zuckerman, K. S. (1998). Hematopoietic Abnormalities in Patients with Cancer. Cancer Control, 5(2_suppl), 6–11. https://doi.org/10.1177/107327489800502s02
Gene Expression Data Analysis of Mutations in Arabidopsis thaliana
Submitted by Group 6:
- Daniel Groves
- Haoming Shi
- Jiahong Zhang
- Varsha Srinivasan
- Yassin Watson
Introduction
Arabidopsis thaliana, the thale cress, mouse-ear cress, or arabidopsis, is a small flowering plant native to Eurasia and Africa. A. thaliana is a popular organism in plant biology and is considered a weed; it is found along the shoulders of roads and in disturbed land [1]. In this analysis, we compared the expression pathway of two A. thaliana mutations of interest: mdf-1 and prp8a-14.
For background on each mutation, Meristem Defective Factor (MDF) has been found to have a role in the organization and maintenance of root meristems [2]. MDF is also required for the correct splicing of numerous transcripts, including transcription factors and cell cycle-associated genes. A loss of MDF function is associated with an impaired transition to mitosis, growth defects, and spontaneous cell death in meristems [3]. Secondly, Arabidopsis pre-RNA processing8a (prp8a-14) modulates splice-site selection. Previous RNA-seq studies have shown that prp8a-14 resulted in hundreds of differentially spliced transcripts and thousands of transcripts with significantly altered levels. Among differentially spliced transcripts, prp8a-14 significantly altered 5′- and 3′-splice-site utilization to favor sites resulting in shorter introns [4].
Each of these mutations alters components involved in RNA splicing. Therefore, we analyzed the expression similarities and differences between experimental groups that induced each mutation separately. The first dataset we selected analyzed the transcriptional and alternative splicing changes in mdf-1 compared to wild type via RNA sequencing. The second dataset contained the result of next-generation sequencing that compared RNA-Seq data from prp8a-14 and Col-0 wild type. Additionally, we examined wild-type groups in each series to identify common expression patterns in control versus mutated groups.
Methods and Materials
The two raw GEO datasets were downloaded from the Sequence Read Archive (SRA), and the reference data and corresponding annotations were downloaded from NCBI. Kallisto [5] was used on the Galaxy server to quantify the abundance of paired-end RNA-Seq transcript data. The program DESeq2 was then used to analyze the differential expression profiles of both datasets.
DESeq2 models differential expression using a generalized Negative Binomial (NB) model by inputting the corresponding transcript per million (TPM) values from the Kallisto abundance files for each gene. The model fits each normalized TPM count and estimates a log2-fold change for each gene. A Wald chi-squared test was then performed to extract a p-value for each gene from the NB regression fit.
After generating statistics from DESeq2, genes with a log2FC > 3 and adjusted p-value < 0.001 were filtered and marked as upregulated genes, while genes with a log2FC < 1/3 and adjusted p-value < 0.001 were identified as downregulated genes using R for both datasets.
Principal Component Analysis (PCA), volcano plots, hierarchical clustering, and pathway analysis were performed to investigate differentially expressed genes. The relative expression level of the two mutated samples, namely, mdf-1 and prp8a-14, represented by the unit TPM, were compared with the respective control sets using PCA. The TPM values were normalized and transformed using the formula log2(TPM+1) to facilitate comparison between samples [6]. PCA was performed to reduce dimensionality. This was done by constructing linear combinations of gene expressions, also known as principal components (PCs), to explain gene expression variation [7].
In addition to TPM values, fold change and adjusted p-values were used in hierarchical clustering. logFC transformation was performed using log2(logFC) to normalize mean counts. Genes were considered differentially expressed if logFC >= 2 and adjusted p-value < 0.05. Hierarchical clustering heatmaps were generated for genes that had logFC >= 2. Volcano plots were generated to visualize and identify statistically significant gene expression changes from normal and mutated conditions. Hierarchical clustering heatmaps were generated to facilitate the visualization and identification of statistically significant gene expression changes among the thousands of genes from normal and mutated conditions.
The gene lists were given as input into ShinyGo to search for pathways and determine the function of differential genes within A. thaliana. In addition, KEGG was utilized to show the details of pathway regulation.
Results
The differential genes were extracted considering the logFC values. and adjusted p-values. It was found that the mdf-1 mutation has 2409 differentially expressed genes and the prp8a-14 mutation has 39 differentially expressed genes.
| Upregulated Gene Count | Downregulated Gene Count | |
| mdf-1 | 1605 | 804 |
| prp8a-14 | 35 | 4 |
PCA
The transformed TPM values of standard and mdf-1/prp8a-14 mutated genes were taken, and PCA was performed using Python. PC retention plots were generated by comparing untreated control genes and mdf-1/prp8a-14 mutated genes (Figures 2,3). The proportion of variance for each PC was calculated by dividing the corresponding eigenvalue by the sum total of the eigenvalues. The PC retention plot for mdf-1 mutations (Figure 2) showed that 76.25% of the total variance was accounted for in PC1 and 22.08% was accounted for in PC2, which means PC1 and PC2 account for a total of 98.33% of the total variance. The PC retention plot for prp8a-14 mutations (Figure 3) showed that 87.80% of the total variance was accounted for in PC1 and 11.62% was accounted for in PC2, which means PC1 and PC2 together account for a total of 99.42% of the total variance. Since PC1 and PC2 alone accounted for the vast majority of the variance, it was concluded that a 2D PCA plot, using only PC1 and PC2, would be a good approximation of the mdf-1/prp8a-14 datasets.
Both the mdf-1 mutation (Figure 4) and the prp8a-14 mutation (Figure 5) showed clustering compared to their replicates. The 2D PCA plot for the mdf-1 mutation again showed that 76.25% of the total variance was accounted for in PC1, while 22.08% was in PC2. The prp8a-14 mutation 2D PCA plot showed that PC1 accounted for 87.80% of the total variance, while PC2 accounted for only 11.62%. Additionally, the PCA plot for the mdf-1 mutation data showed close clustering of the controls and mdf-1 mutations. Similarly, the replicates for the controls and the prp8a-14 mutations cluster together.
Volcano Plots
To better investigate mutations that showed an absolute logFC >= 2, volcano plots were constructed using Python to show upregulated and downregulated genes after mutations (Figure 6,7; adjusted p-value < 0.05, logFC >= 2). The volcano plot highlighted all significant genes for mdf-1 (Figure 6) and prp8a-14 (Figure 7) mutations. In the plots below, the grey dots represent the genes that do not pass through the adjusted p-value and logFC thresholds. The red dots denote upregulated genes, and the green dots denote downregulated genes.
Hierarchical Clustering
To show the spread and intensity of mutation detected, we performed hierarchical clustering by constructing heatmaps for mdf-1/prp8a-14 showing the upregulated and downregulated genes after mutation (Figures 8, 9). Each hierarchical cluster was cross-compared to determine whether there were any shared upregulated or downregulated genes across species with standard normalized values ranging from 0 to 1. Figure 8 reveals differential gene expression between mdf-1 mutations and controls in the heatmap. A handful of genes are upregulated (indicated in red), while a more significant number is downregulated (indicated in blue). Likewise, Figure 9 visually represents the differential gene expression between prp8a-14 mutations and controls. A similar trend is observed for the prp8a-14 mutations
Pathway Analysis
To further analyze the functions of differential genes in each dataset, differential gene lists are extracted and given as input into ShinyGO (http://bioinformatics.sdstate.edu/go/). Based on the search result from the differential gene set from mdf-1, the pathway for sesquiterpenoid and triterpenoid biosynthesis and the pathway for cultin, submarine, and wax biosynthesis were observed to be the implicated pathways with high fold enrichment. Regarding the prp8a-14 dataset, the top pathways implicated are monoterpenoid biosynthesis and sesquiterpenoid and triterpenoid biosynthesis. Moreover, pathways were searched in KEGG for detailed information on gene regulation.
Discussion
In this study, we compared the transcriptional and alternative splicing changes in mdf-1 and prp8a-14 mutations to the wild type of A. thaliana to determine the shared genes between species that are differentially expressed. After cross-comparing the two datasets, it was found that AT5G37990, AT4G11320, and AT5G38020 are upregulated in mdf-1/prp8a-14 mutations, while only 1 shared gene, AT1G61800, was downregulated between two mutations (Figures 8, 9).
Based on the initial comparison from ShinyGO, it was found that the pathway for sesquiterpenoid and triterpenoid biosynthesis is dominant in both mdf-1 and prp8a-14 datasets. To have a detailed look at the pathway, we turned to KEGG by casting our differential genes into the pathway for sesquiterpenoid and triterpenoid biosynthesis. Further, genes like AT4G13280 and AT5G36150 are significantly dominant in the expression for mdf-1 and prpa8a-14 datasets, which may indicate some relations between these two mutations.
Sesquiterpenoids (C15 terpenoids) are a group of terpenoids consisting of three isoprene units. They are derived from farnesyl diphosphate (FPP) and can be cyclized to produce various skeletal structures. Sesquiterpenoid biosynthesis begins with the loss of diphosphate from FPP under the action of sesquiterpene synthesis enzymes, generating an allylic cation that is highly susceptible to intramolecular attacks. Cyclization of the farnesyl cation may occur onto either of the remaining double bonds, resulting in the formation of 6-, 10-, or 11-membered rings. Additionally, the enhanced expression in this regulation pathway could catalyze sesquiterpenoids to (Z)-γ-Bisabolene and triterpenoid to tirucalla-7, 24-dien-3β-ol. (Z)-γ-Bisabolene is expressed in root and tirucalla-7, 24-dien-3β-ol is related to plant metabolite. Based on the initial literature review for the regulation and function of the two mutations [4, 8], it can be gathered that both mutation types resulted in changes in root cell growth and elongation. Therefore, these two molecules may play a role in prp8a and mdf mutation types and need a designed experiment in the future to prove their functions in the mutation.
References
[1] Hoffmann, M. H. (2002). Biogeography of arabidopsis thaliana (L.) heynh. (Brassicaceae). Journal of Biogeography, 29(1), 125–134. https://doi.org/10.1046/j.1365-2699.2002.00647.x
[2] Helen. I. Thompson, H. T. (2020). Meristem-defective, a novel splicing factor essential for root meristem development in Arabidopsis thaliana. https://doi.org/10.26226/morressier.5ebd45acffea6f735881b0f9
[3] Barton, M. K. (1993). Genetic analysis of meristem structure and function in Arabidopsis thaliana. Cellular Communication in Plants, 69–73. https://doi.org/10.1007/978-1-4757-9607-0_11
[4] Llinas, R. J., Xiong, J. Q., Clark, N. M., Burkhart, S. E., & Bartel, B. (2022). An arabidopsis pre-RNA PROCESSING8A (prp8a) missense allele restores splicing of a subset of mis-spliced mRNAs. Plant Physiology, 189(4), 2175–2192. https://doi.org/10.1093/plphys/kiac221
[5] Bray, N. L., Pimentel, H., Melsted, P., & Pachter, L. (2016). Near-optimal probabilistic RNA-seq quantification. Nature Biotechnology, 34(5), 525â527. https://doi.org/10.1038/nbt.3519
[6] Zhao, Y., Li, MC., Konaté, M.M. et al. TPM, FPKM, or Normalized Counts? A Comparative Study of Quantification Measures for the Analysis of RNA-seq Data from the NCI Patient-Derived Models Repository. J Transl Med 19, 269 (2021). https://doi.org/10.1186/s12967-021-02936-w
[7] Shuangge Ma, Ying Dai, Principal component analysis based methods in bioinformatics studies, Briefings in Bioinformatics, 12 (6), November 2011, 714–722, https://doi.org/10.1093/bib/bbq090
[8] Casson, S. A., Topping, J. F., & Lindsey, K. (2009). MERISTEM‐DEFECTIVE, an RS domain protein, is required for the correct meristem patterning and function in Arabidopsis. The Plant Journal, 57(5), 857-869.
Evaluating the Impact of Two Different Types of Drugs as Adjunct Therapy for the Treatment of Ovarian Cancer Using RNA-Seq Analysis
Group 7: Zijin Guo, Yue Han, Ying Li, Saanika Tambe
INTRODUCTION AND BACKGROUND
Epithelial Ovarian Cancer (EOC) is one of the most malignant cancers and is the fourth most common lethal cancer in women across the globe [1]. The diagnosis is usually at a late stage as presently there are no effective screening strategies for the early detection of EOC [2]. The treatment methods are heavily dependent on the stage at which the disease is diagnosed [3].
As most patients with advanced ovarian cancer cannot be treated by surgery alone, there is an increasing focus on chemotherapy [4] and targeted therapy [3] for treatment. Common chemotherapeutic drugs are platinum-containing drugs like cisplatin or carboplatin in combination with the taxane, paclitaxel [4]. Targeted therapy includes the use of targeted biologic agents such as bevacizumab [3].
Use of Bevacizumab
Bevacizumab is an anti-VEGF (Vascular Endothelial Growth Factor) agent and a critical angiogenic promoter [5]. When used as a single-agent therapy, bevacizumab has shown to limit the progression of cancer in Phase 2 studies [4] and is under trial for frontline therapy in combination with platinum-taxane chemotherapy [5].
Use of Luteolin
Luteolin is a flavonoid that has been demonstrated to have anti-tumor activity in various types of cancer [6]. In case of ovarian cancer, chemotherapy resistance often obstructs effective treatment. Luteolin’s anti-tumor mechanisms include cell-cycle arrest, metastatic inhibition and angiogenesis.
OBJECTIVE
Our main objective is to study the effects of Bevacizumab and Luteolin on ovarian cancer cell lines individually using RNA-Seq analysis. Both these are adjunct therapies that are used with chemotherapeutic drugs such as cisplatin or carboplatin. We performed gene expression analysis for the two datasets of bevacizumab and luteolin. The cell lines used were from human ovarian cancer tumor cells.
In practice, a combination of various drugs and treatment methods are used to treat ovarian cancer. Luteolin and bevacizumab are considered for use in combination with the traditional chemotherapy method to make the treatment more effective.
METHODS
To study the two drugs that are used to treat ovarian cancer in different pathways, we choose two datasets, both from NCBI Gene Expression Omnibus. Both of the datasets are aimed to study the gene expression of the ovarian cancer cell lines after they are treated with the drugs they are studying. One dataset is from the study of Luteolin treatment. The A2780 cell line is treated with luteolin and DMSO (control), and then gene expression profiling analysis is performed on the RNA-seq data. Another dataset is from the study about gamma irradiation with bevacizumab. The OVISE cells were treated with bevacizumab and incubated for 192 h and then total RNA was extracted from these cells 2 hours after treatment with 2 Gy of gamma irradiation for RNA sequencing. The RNA-seq data is then used to identify the differentially expressed genes. The samples of this comparative gene expression profiling analysis are outlined in [Tab. 1].
Data sets
TABLE I
A LIST OF RNA-Seq DATASETS USED IN THIS STUDY
| Dataset GEO Accession | Sample GEO Accession | Sample Type | Cell Line | Treatment |
| GSE212598 | GSM6538260 | SRA | A2786 | DMSO |
| GSE212598 | GSM6538261 | SRA | A2787 | DMSO |
| GSE212598 | GSM6538262 | SRA | A2788 | DMSO |
| GSE212598 | GSM6538263 | SRA | A2789 | luteolin |
| GSE212598 | GSM6538264 | SRA | A2790 | luteolin |
| GSE212598 | GSM6538265 | SRA | A2791 | luteolin |
| GSE203044 | GSM6152910 | SRA | OVISE | PBS+gamma irradiation |
| GSE203044 | GSM6152911 | SRA | OVISE | PBS+gamma irradiation |
| GSE203044 | GSM6152912 | SRA | OVISE | PBS+gamma irradiation |
| GSE203044 | GSM6152913 | SRA | OVISE | bevacizumab+gamma irradiation |
| GSE203044 | GSM6152914 | SRA | OVISE | bevacizumab+gamma irradiation |
| GSE203044 | GSM6152915 | SRA | OVISE | bevacizumab+gamma irradiation |
Workflow
Above is the bioinformatics pipeline we developed to identify differentially expressed genes and perform gene expression profiling analysis on the identified genes. After datasets were obtained, we used FastQC to perform a quality test for raw sequencing data and MultiQC to aggregate the results. The reference genome we chose is GRCh38.p13 (GenBank GCA_000001405.28), which served as reference transcriptome for all of the following analysis. Kallisto quant [7] was used to quantify abundances of transcripts from RNA-seq data. The count tables generated by Kallisto quant were then used to run DESeq2 [8] to generate differentially expressed features as well as normalized counts and useful statistical test results that can be used to perform the false discovery rate control following Benjamini-Hochberg procedure [9]. To see the clusters of samples based on their similarity, principal component analysis (PCA) was performed and a plot showing the first two principal components was made. For hierarchical clustering, heatmaps were generated by heatmap tools in Galaxy. Pathway analysis for different datasets was performed using PANTHER.
Data and Statistical Analysis
Principal component analysis and hierarchical clustering were performed using DESeq2 and heatmap tools in Galaxy. PCA plots, normalized counts, and heatmaps were generated. Following Benjamini-Hochberg procedure, a level of δ = 0.05 was chosen to reduce the rate of Type I error. Genes that have a fold change significantly different from 1 were regarded as being differentially expressed. All genes were ranked by p-value from small to large. The genes with ranks lower than or equal to j were determined as differentially expressed where pj <= (j/m) * δ and m is the number of total tests. Results were uploaded to PANTHER for pathway analysis to identify affected cellular pathways.
RESULT
Differentially expressed genes
Following the workflow and statistical approaches described in the previous section, a list of differentially expressed genes was obtained for each dataset. For the bevacizumab dataset, a total of 5 genes were found to be differentially expressed, of which 1 was upregulated and 4 were downregulated. For the luteolin dataset, a total of 11268 genes were found to be differentially expressed, of which 6119 were upregulated and 5149 were downregulated.
PCA for Luteolin and Bevacizumab and Their Respective Controls
Fig. 1. PCA for Luteolin and DMSO (Control).
Fig. 2. PCA for Bevacizumab and PBS (Control).
The gene counts table generated by Kallisto quant was uploaded to DESeq2 on the Galaxy server. The PCA plots were generated. A PCA plot of all 12 samples was drawn, but significant differences between the two control groups were observed. As a result, PCA was performed on the 2 datasets separately. For Luteolin and DMSO, 97% of the variance is retained by the first two principal components. A clear separation is present between the control and treated group for luteolin treatment as illustrated in [Fig. 1].
For treatment with bevacizumab, no clear separation between the control and treated group is observed and a lower 50% variance is retained by the first two principal components, which is obvious in [Fig. 2]. This is not surprising given that only 5 genes were found to be differentially expressed in this set of data.
Hierarchical Clustering for Luteolin and Bevacizumab and Their Respective Controls
Fig. 3. Heatmaps for hierarchical clustering for luteolin and bevacizumab and their respective controls. (a) Red color indicates more gene counts. (b) Blue color indicates less gene counts.
A hierarchical clustering is performed on differentially expressed genes in the two datasets – Luteolin & DMSO and Bevacizumab & PBS. Euclidean distance is used to quantify the relatedness of different data points and a complete-linkage clustering method is used. To construct the map, the gene counts are normalized and log2 transformed.
Again, as expected, we do not observe much difference between bevacizumab and its control due to the small number of differentially expressed genes (five in number). In [Fig. 3, left], the patterns of the color on the left and on the right are not very different.
For Luteolin vs DMSO, there is some amount of differential expression observed in the genes. However, not much concrete inference can be taken from this alone because the patterns in [Fig. 3, right] are not clear enough. The darker blue on luteolin side might indicate some down-regulated gene clusters.
Overall, some clusterings of up/downregulation of genes were observed. To study the functionality of those genes and make biological inferences, we went on to perform pathway analysis.
Pathway Analysis
We used the PANTHER Overrepresentation test to analyze the Luteolin vs. DMSO dataset and the Bevacizumab vs. PBS dataset, with the input file being the significant differentially expressed genes. The annotation data set was the Reactome pathway, with a test type of Fisher’s Exact, and we used the Bonferroni correction for multiple hypothesis testing.
There were no overlap pathways between the two datasets because the bevacizumab dataset did not show any significant overrepresentation pathways that had a p-value smaller than 0.05 and this was not surprising given the small number of differentially expressed genes.
On the other hand, the Luteolin dataset showed some significant overrepresentation of upregulated pathways in metabolism of RNA (R-HSA-8953854), cell cycle(R-HSA-1640170), and transcription(R-HSA-74160). The Luteolin dataset also showed some significant overrepresentation of downregulated pathways in mainly keratinization (R-HSA-6805567), GPCR ligand binding(R-HSA-500792), and ClassA/1 pathway(R-HSA-373076).
TABLE II
TOP 10 UP-REGULATED GENES DURING LUTEOLIN TREATMENT
TABLE III
TOP 10 DOWN-REGULATED GENES DURING LUTEOLIN TREATMENT
RESULT DISCUSSION
Since the two datasets chosen are both published recently (one in May 2022 and the other in September 2022), we were not able to find the publication related to the specific datasets. However, we went into relevant literature about luteolin and bevacizumab, including past publications by authors of the two datasets.
Bevacizumab is a Vascular Endothelial Growth Factors (VEGF)-A targeting monoclonal antibody that inhibits angiogenesis [6]. Solid tumor growth requires formation of new blood vessels, which is activated by VEGF binding to VEGF receptor tyrosine kinases (VEGFR1-3). And bevacizumab acts as an inhibitor of VEGF-A. This signaling pathway occurs in the endothelial cells, which explains the lack of differential expression in tumor cells in this experiment.
Although pathway analysis of bevacizumab did not return much useful information, several differentially expressed genes in the dataset may be interesting. The MDM2 gene encodes a protein that is an important negative regulator of the p53 tumor suppressor [10]; TOP2A gene expression has been found to be a key indicator of chemotherapy response [11].
Luteolin, the other drug in this study, is a flavonoid recognized as an anticancer agent. It has been shown to inhibit cancer development by inhibition of proliferation of tumor cells, protection from carcinogenic stimuli, and activation of cell cycle arrest, as well as inducing apoptosis through different signaling pathways [12]. We found several pathways related to cell cycle being upregulated, which agrees with the literature. Pathways related to response to stress and stimuli were also found to be upregulated, which explain why metabolism of RNA, transcription, and metabolism of proteins are also upregulated. These correspond to the literature that luteolin increases levels of intracellular reactive oxygen species by activation of lethal endoplasmic reticulum stress response and mitochondrial dysfunction, and by activation of ER stress-associated protein expressions [12]. Many pathways related to signaling (e.g. GPCR ligand binding), however, were found to be downregulated, which disagrees with literature.
In conclusion, we were not able to find many differentially expressed genes or pathways in the OVISE cell during treatment with bevacizumab, but many were found in the A2780 cell line during treatment with luteolin. The possible biological implications were discussed by performing a literature search on the role of bevacizumab and luteolin in cancer treatment. However, due to the mechanism of action bevacizumab, we were not able to perform a comparison between the two types of treatment from the gene expression data.
REFERENCES
[1] G. C. Jayson, E. C. Kohn, H. C. Kitchener, and J. A. Ledermann, “Ovarian cancer,” The Lancet, vol. 384, no. 9951, pp. 1376–1388, Oct. 2014, doi: 10.1016/S0140-6736(13)62146-7.
[2] U. A. Matulonis, A. K. Sood, L. Fallowfield, B. E. Howitt, J. Sehouli, and B. Y. Karlan, “Ovarian cancer,” Nature Reviews Disease Primers 2016 2:1, vol. 2, no. 1, pp. 1–22, Aug. 2016, doi: 10.1038/nrdp.2016.61.
[3] D. Jelovac and D. K. Armstrong, “Recent progress in the diagnosis and treatment of ovarian cancer,” CA Cancer J Clin, vol. 61, no. 3, pp. 183–203, May 2011, doi: 10.3322/CAAC.20113.
[4] L. R. Kelland, “Emerging drugs for ovarian cancer,” http://dx.doi.org/10.1517/14728214.10.2.413, vol. 10, no. 2, pp. 413–424, May 2005, doi: 10.1517/14728214.10.2.413.
[5] R. A. Burger, “Experience with bevacizumab in the management of epithelial ovarian cancer,” Journal of Clinical Oncology, vol. 25, no. 20, pp. 2902–2908, Jul. 2007, doi: 10.1200/JCO.2007.12.1509.[6] H. Wang, Y. Luo, T. Qiao, Z. Wu, and Z. Huang, “Luteolin sensitizes the antitumor effect of cisplatin in drug-resistant ovarian cancer via induction of apoptosis and inhibition of cell migration and invasion,” J Ovarian Res, vol. 11, no. 1, Nov. 2018, doi: 10.1186/S13048-018-0468-Y.
[6] Garcia J, Hurwitz HI, Sandler AB, Miles D, Coleman RL, Deurloo R, Chinot OL. Bevacizumab (Avastin®) in cancer treatment: A review of 15 years of clinical experience and future outlook. Cancer Treat Rev. 2020 Jun;86:102017. doi: 10.1016/j.ctrv.2020.102017. Epub 2020 Mar 26. PMID: 32335505.
[7] N. L. Bray, H. Pimentel, P. Melsted, and L. Pachter, “Near-optimal probabilistic RNA-seq quantification,” Nature Biotechnology, vol. 34, no. 5, pp. 525–527, 2016.
[8] M. I. Love, W. Huber, and S. Anders, “Moderated estimation of fold change and dispersion for RNA-seq data with deseq2,” Genome Biology, vol. 15, no. 12, 2014.
[9] Y. Benjamini and Y. Hochberg, “Controlling the false discovery rate: A practical and powerful approach to multiple testing,” Journal of the Royal Statistical Society: Series B (Methodological), vol. 57, no. 1, pp. 289–300, 1995.
[10] Oliner JD, Kinzler KW, Meltzer PS, George DL, Vogelstein B (July 1992). “Amplification of a gene encoding a p53-associated protein in human sarcomas”. Nature. 358 (6381): 80–3. Bibcode:1992Natur.358…80O. doi:10.1038/358080a0. hdl:2027.42/62637. PMID 1614537. S2CID 1056405.
[11] Burgess DJ, Doles J, Zender L, Xue W, Ma B, McCombie WR, Hannon GJ, Lowe SW, Hemann MT. Topoisomerase levels determine chemotherapy response in vitro and in vivo. Proc Natl Acad Sci U S A. 2008 Jul 1;105(26):9053-8. doi: 10.1073/pnas.0803513105. Epub 2008 Jun 23. PMID: 18574145; PMCID: PMC2435590.
[12] Lin Y, Shi R, Wang X, Shen HM. Luteolin, a flavonoid with potential for cancer prevention and therapy. Curr Cancer Drug Targets. 2008 Nov;8(7):634-46. doi: 10.2174/156800908786241050. PMID: 18991571; PMCID: PMC2615542.
Analysis of Differentially Expressed Genes in Human Fibroblasts vs Human Keratinocytes when exposed to Interferon-gamma
Group U4: Enysa Moore, Rebecca Oh, Trip Hash, Caleb Connors, Sekou Noble-Kuchera, Cristina Enériz Bueno
November 7, 2022
BIOL 4150-A
Introduction
Interferons are substances naturally produced by white blood cells and certain other cell types that help the body’s immune system fight infection and other diseases. However, they can also be produced artificially in laboratories to use them for treating diseases.
Figure 1: basic model of interferon mechanism in human immune response. [20]
In 1965, the gamma interferon was discovered. Due to its cell type-specific antiviral and antiangiogenic activities, early clinical trials of this cytokine began to evaluate its therapeutic potential. The initial studies focused on its tolerability and pharmacology and served to determine its antitumor and anti-infection activities. In later years, interferon-gamma has been used to treat a variety of clinical indications, some very crucial such as skin cancer or severe atopic dermatitis.
In this study, we analyzed the gene expression profiles of two different types of human skin cells, fibroblast and keratinocytes, after IFN-γ treatment aiming to answer the following experimental questions:
1.) How do fibroblasts respond to treatment with interferon-gamma?
2.) How do keratinocytes respond to treatment with interferon-gamma?
3.) Are there any significant differences in gene expression or influenced cellular pathways between the two treated and untreated cell types?
Background
Interferon-gamma (IFN-γ) is a type of interferon (in fact, the only type II interferon) that was first discovered by E.F. Wheelock as a product of human leukocytes stimulated with phytohemagglutinin. From that moment on, it has been historically known as “the immune interferon” due to its antiviral, immunoregulatory, and anti-tumor properties.
Researchers found that IFN-γ is both involved in the innate and adaptive immune response. It is produced mainly by natural killer cells and natural killer T cells as part of the innate immune response, and by cytotoxic T lymphocyte effector T cells once antigen-specific immunity develops as part of the adaptive immune response. IFN-γ is also known to be produced by non-cytotoxic innate lymphoid cells, a family of immune cells first discovered in the early 2010s.
Due to its positive effects and its importance in the functioning of the immune system, the administration of IFN-γ has been studied, among others, as a treatment for skin-related diseases. Although not officially approved, it has been shown to be effective in treating patients with moderate to severe atopic dermatitis. Specifically, recombinant IFN-γ therapy has shown promise in patients with lowered IFN-γ expression, such as those with a predisposition to herpes simplex virus, and pediatric patients.
Besides, a potential use in immunotherapy has been observed. In its natural form, IFN-γ has been shown to increase the anti-proliferative state of cancer cells while increasing immune recognition and removal of pathogenic cells by upregulating MHC I and MHC II expression. Moreover, it has been found that IFN-γ reduces metastasis in tumors by upregulating fibronectin, which negatively impacts tumor architecture.
Regardless of these findings, IFN-γ has not been approved yet for the treatment in any cancer immunotherapy. However, improved survival was observed when IFN-γ was administered to melanoma patients, making it a focus of attention for researchers in the skin and cancer field.
As a multifunctional, immunomodulatory cytokine, IFN-γ alters transcription in up to 30 genes, producing a variety of physiological and cellular responses. Because of this reason, a cross-comparison in the gene expression between human fibroblast and keratinocytes, both treated with interferon-γ, was made in this study. Our aim was to find shared differentially expressed genes and determine if similar cellular pathways were affected in the two treated and untreated cell types.
Methods
Figure 2: pipeline used in data processing and analysis.
Datasets
Each dataset was downloaded from the NCBI’s Gene Expression Omnibus (GEO) database and both use the Affymetrix Human Genome U95 Version 2 Array. Fibroblasts and keratinocytes were both exposed to an interferon-gamma treatment. Samples used in this cross-comparison analysis are delineated in Table 1 and Table 2.
Table 1 (Microarray dataset GSE3920): Affymetrix Human Genome U133A Array of fibroblasts treated with interferon-gamma for five hours.
Table 2 (Microarray dataset GSE440): Affymetrix Human Genome U95 Version 2 Array of keratinocytes treated with interferon-gamma for up to 48 hours.
Microarray Data Analysis
ExAtlas was used for the meta-analysis and visualization of gene expression for the datasets shown above using the Gene Set Enrichment (GSE) tool, and the downloaded replicates were renamed in order to create a matrix. The matrix was generated using ExAtlas allowing for differentially expressed genes, both under-expressed and over-expressed, to be targeted. The conditions for generating the differentially expressed genes were as follows:
- Baseline: HU_con_1
- FDR Threshold: .05
- Fold-Change Value: 2
- Expression Threshold: 0
- Specific Gene Filter: 0
Differentially expressed genes and their log values for gene expression were determined using GEO2R, extracting those with an α-value below 10-4 for pathway analysis. As already stated, for the dataset with different cell types and interferon treatments (GSE 3920), only fibroblasts treated with IFN-y and control fibroblasts were considered (the bottom 5 datasets).
Pathway analysis was performed using Cytoscape, with the differentially expressed genes found above. The pathway analysis was done by following the differential expression analysis guide on Cytoscape. (step-by-step guide on Cytoscape here [19]). For this, a STRING protein query was used, our p values were integrated with the data, and the relevant pathways were found and extracted. Redundant paths were filtered out and removed.
Results
1.) GEO2R Run for Differentially Expressed Genes:
Table 3: GEO2R Run for GSE440
Table 4: GEO2R Run for GSE 3920
The displayed GEO2R results show analyzed differentially expressed genes between the control samples and the treatment samples. It shows the up-regulated and down-regulated genes, along with the total amount of differentially expressed genes for each dataset.
2.) Differentially Expressed Genes Found in Both Data Sets:
| Gene | Description |
| CXCL9,10,11 | CXCL9,10,11 are closely related cytokines that specifically bind to their receptor CXCR3. This axis regulates immune cell migration, differentiation, and activation, leading to tumor suppression. However, there are some reports that show the involvement of this axis in tumor growth and metastasis. Thus, their pro- or anti-tumor effects remain controversial and a better understanding is necessary to develop effective cancer control [11] [12]. |
| IRF1 | Transcriptional regulator and tumor suppressor. Plays a role in response to bacteria, immune response, and DNA damage response. Defects in this gene have been associated with gastric cancer, lung cancer, and myelogenous leukemia [13]. |
| IFIT3 | IFIT3 is an interferon-induced protein that enables identical protein binding. It is involved in the negative regulation of cell population proliferation and virus response [14]. |
3.) Hierarchical Clustering (heat maps):
Figure 3: hierarchical clustering heat map for GSE 3920
Figure 4: hierarchical clustering heat map for GSE 440.
4.) Principle Component Analysis (PCA):
Figure 5: PCA plot of GSE 3920 (Endothelial Cells treated with IFN-g). Left: two controls were compared against three IFN-γ treated samples.
Table 5: PC values for GSE 3920. 99.681 percent of the variance was captured with three principal components.
Figure 6: PCA plot of GSE440 (Keratinocytes treated with IFN-g). Five controls were compared against four IFN-γ treated samples.
Table 6: PC values for GSE 400. 70.517 percent of the variance was captured with three principal components.
5.) Pathway Analysis:
Figure 7: pathway analysis for GSE 440.
Figure 8: pathway analysis for GSE 3920.
6.) Cytoscape Analysis Results – Top Pathways:
Table 7: GSE 440 top pathways, and the genes that influence them.
Table 7: GSE 440 top pathways reduced by Cytoscape
Table 9: GSE 3920 top pathways, and the genes that influence them.
Table 10: GSE 3920 top pathways reduced by Cytoscape
Both datasets showed differential expression in immune signaling response, tuberculosis immune response, and type 2 interferon signaling pathways, as expected from the literature. Differentially expressed pathways in the endothelial cell dataset include cytokine signaling, cell signaling responses, and other immune or interferon-related pathways. Pathways in the keratinocyte dataset include cytokine TNF signaling, signal transduction, chemokine receptor, oncostatin M signaling, lipopolysaccharide response, and tuberculosis immune response pathways. Overall, the pathways that were differentially expressed were very similar.
Discussion
In this study we wanted to see if there was a difference in the differentially expressed genes in Keratinocytes in comparison to Fibroblasts when exposed to Interferon-gamma. The comparison was made using datasets GSE3920 and GSE440, and while we found that both cell types expressed genes involved in the immune signaling response as well as the Interferon type 2 signaling pathway, there was an unexpecting finding that both had genes involved in the tuberculosis immune response pathway. Both the Keratinocytes and Fibroblasts shared differential expressions in genes IRF1, IFIT3, and PSMB8, all involved in the tuberculosis immune system response pathway.
Interferon gamma (INF-γ) is a pleiotropic cytokine that has been utilized to treat both tumors and autoimmune diseases due to its immune system regulation and host defense properties [9]. The indication that INF-γ causes genes in the Tuberculosis immune system response pathway to be expressed in both keratinocytes and fibroblast led us to hypothesize that tuberculosis patients may benefit from treatment with INF- γ. Further studies have expanded upon this hypothesis noting that the use of Interferon-gamma may combat the damage that multi-drug resistant tuberculosis (MDR-TB) causes by regulating the patient’s immune through activation of microphases, inducing the release of nitrous oxide, and serving as an inducer for second class histocompatibility complex (MHC Class II) [10].
Through this study we found that further research into treating MDR-TB with Interferon- γ should be conducted in hopes of finding a method of curing tuberculosis in the future as well as providing care to patients that are suffering. Our study confirmed that in Keratinocytes and Fibroblasts, two tissue types found in abundance within human lungs, there is gene expression in the tuberculosis immune response pathway when exposed to INF- γ, however we concluded that more research needs to be conducted in order to see long-term potential side effects from treatment as well as if the gene continues to be expressed over longer periods of time.
References
[1] Wheelock EF (July 1965). “Interferon-Like Virus-Inhibitor Induced in Human Leukocytes by Phytohemagglutinin”. Science. 149 (3681): 310–311. Bibcode:1965Sci…149..310W. doi:10.1126/science.149.3681.310. PMID 17838106. S2CID 1366348.
[2] Schoenborn JR, Wilson CB (2007). “Regulation of Interferon‐γ During Innate and Adaptive Immune Responses”. Regulation of interferon-gamma during innate and adaptive immune responses. Advances in Immunology. Vol. 96. pp. 41–101. doi:10.1016/S0065-2776(07)96002-2. ISBN 978-0-12-373709-0. PMID 17981204.
[3] Artis D, Spits H (January 2015). “The biology of innate lymphoid cells”. Nature. 517 (7534): 293–301. Bibcode:2015Natur.517..293A. doi:10.1038/nature14189. PMID 25592534. S2CID 4386692.
[4] Wells M, Seyer L, Schadt K, Lynch DR (December 2015). “IFN-γ for Friedreich ataxia: present evidence”. Neurodegenerative Disease Management. 5 (6): 497–504. doi:10.2217/nmt.15.52. PMID 26634868.
[5] Seyer L, Greeley N, Foerster D, Strawser C, Gelbard S, Dong Y, et al. (July 2015). “Open-label pilot study of interferon gamma-1b in Friedreich ataxia”. Acta Neurologica Scandinavica. 132 (1): 7–15. doi:10.1111/ane.12337. PMID 25335475. S2CID 207014054.
[6] Lynch DR, Hauser L, McCormick A, Wells M, Dong YN, McCormack S, et al. (March 2019). “Randomized, double-blind, placebo-controlled study of interferon-γ 1b in Friedreich Ataxia”. Annals of Clinical and Translational Neurology. 6 (3): 546–553. doi:10.1002/acn3.731. PMC 6414489. PMID 30911578.
[7] Brar K, Leung DY (2016). “Recent considerations in the use of recombinant interferon gamma for biological therapy of atopic dermatitis”. Expert Opinion on Biological Therapy. 16 (4): 507–514. doi:10.1517/14712598.2016.1135898. PMC 4985031. PMID 26694988.
[8] Kak G, Raza M, Tiwari BK (May 2018). “Interferon-gamma (IFN-γ): Exploring its implications in infectious diseases”. Biomolecular Concepts. 9 (1): 64–79. doi:10.1515/bmc-2018-0007. PMID 29856726. S2CID 46922378.
[9] Jorgovanovic D, Song M, Wang L, Zhang Y (2020-09-29). “Roles of IFN-γ in tumor progression and regression: a review”. Biomarker Research. 8 (1): 49. doi:10.1186/s40364-020-00228-x. PMC 7526126. PMID 33005420.
[10] Ivashkiv LB. IFNγ: signalling, epigenetics and roles in immunity, metabolism, disease and cancer immunotherapy. Nat Rev Immunol. 2018 Sep;18(9):545-558. doi: 10.1038/s41577-018-0029-z. PMID: 29921905; PMCID: PMC6340644.
[11] Tokunaga R, Zhang W, Naseem M, Puccini A, Berger MD, Soni S, McSkane M, Baba H, Lenz HJ. CXCL9, CXCL10, CXCL11/CXCR3 axis for immune activation – A target for novel cancer therapy. Cancer Treat Rev. 2018 Feb;63:40-47. doi: 10.1016/j.ctrv.2017.11.007. Epub 2017 Nov 26. PMID: 29207310; PMCID: PMC5801162.
[12] Blanco, G., Puiggros, A., Sherry, B., Nonell, L., Calvo, X., Puigdecanet, E., Chiu, P. Y., Kieso, Y., Ferrer, G., Arnal, M., Rodríguez-Rivera, M., Gimeno, E., Abella, E., Rai, K. R., Abrisqueta, P., Bosch, F., Ferrer, A., Chiorazzi, N., & Espinet, B. (2019). Deciphering the CXCL9-CXCL10-CXCL11/CXCR3 axis in CLL-like monoclonal B-cell lymphocytosis and chronic Lymphocytic leukemia: A new target for immune activation? Blood, 134(Supplement_1), 3029–3029. https://doi.org/10.1182/blood-2019-122061
[13] IRF1 Interferon Regulatory Factor 1 [Homo Sapiens (Human)] – Gene – NCBI. https://www.ncbi.nlm.nih.gov/gene/3659. Accessed 7 Nov. 2022.
[14] IFIT 3 interferon induced protein with tetratricopeptide repeats 3 [Homo Sapiens (Human)] – Gene – NCBI. https://www.ncbi.nlm.nih.gov/gene/3659. Accessed 7 Nov. 2022 https://www.ncbi.nlm.nih.gov/gene/3437.
[15] Miller C, Maher S, Young H. Clinical use of Interferon-γ. Annals of the New York Academy of Sciences. 2009 Dec 14. Vol. 1182. pp. 69-79. https://doi.org/10.1111/j.1749-6632.2009.05069.x
[16] Berns S, Isakova J, Pekhtereva P. Therapeutic potential of interferon-gamma in tuberculosis. ADMET. 2022 Feb 14. Vol 10. pp. 63-73. doi: 10.5599/admet.1078
[17] Indraccolo S, Pfeffer U, Minuzzo S, Esposito G et al. Identification of genes selectively regulated by IFNs in endothelial cells. J Immunol 2007 Jan 15;178(2):1122-35. PMID: 17202376
[18] Banno T, Adachi M, Mukkamala L, Blumenberg M. Unique keratinocyte-specific effects of interferon-gamma that protect skin from viruses, identified using transcriptional profiling. Antivir Ther 2003 Dec;8(6):541-54. PMID: 14760888
[19] DE Genes Network Analysis. https://cytoscape.org/cytoscape-tutorials/protocols/differentially-expressed-genes/#/. Accessed 7 Nov. 2022.
[20] Oiseth S, Jones L, Maza E: Interferons. Concise Medical Knowledge. https://www.lecturio.com/concepts/interferons/. Accessed 7 Nov. 2022.
Investigating the role of CBX2 in breast cancer through the analysis of differentially expressed genes
Sojeong Gwon, Fang Shi, Joseph Tsenum, Cheng Zhang, Siming Zhao
Introduction
Proteins that regulate the gene expression process can be potential therapeutic targets for cancer. CBX2 is one of the important components of Polycomb Repressive Complex 1 (PRC1), which is involved in histone modification and chromatin remodeling (Jangal, Lebeau & Witcher, 2019). CBX2 expression is elevated in triple-negative breast cancer (TNBC), and its knockdown in TNBC models reduced cell numbers by inhibiting proliferation (Bilton et al., 2022). Similar effects have been observed in estrogen receptor (ER)-positive breast cancer (Bilton et al., 2022). Our chosen dataset is composed of sequenced mRNA profiles of two different types of breast cancer cells (MDA-MB-231 for TNBC and MCF-7 for ER-positive) where CBX2 has been knocked down with siRNA. Our project aims to analyze the effect of CBX2 knockdown on gene expression in MDA-MB-231 and MCF-7 cells. MCF-7 and MDA-MB-231 are common breast cancer cell lines for ER-positive and TNBC, respectively. MCF-7 is estrogen- and EGF- (Epidermal Growth Factor) dependent, while MDA-MB-231 is hormone-independent. By analyzing the difference in the transcriptomic profiles of these two cell lines with CBX2 knockdown, we can further identify potential CBX2-regulated genes across two different breast cancer types.
Methods
1. Datasets
- CBX2 Knockdown in MCF-7 cells
- GSE198418
- Platform: Illumina NovaSeq 6000 (Homo sapiens)
- Control: siSCR (non-silencing scrambled control siRNA)
- CBX2 Knockdown: siCBX2#3 (CBX2 targeting siRNA)
| GSM number | SRR number | Treatment |
| GSM5946434 | SRR18299859 | siSCR Replicate1 |
| GSM5946435 | SRR18299860 | siCBX2#3 Replicate1 |
| GSM5946436 | SRR18299861 | siSCR Replicate2 |
| GSM5946437 | SRR18299862 | siCBX2#3 Replicate2 |
| GSM5946438 | SRR18299863 | siSCR Replicate3 |
| GSM5946439 | SRR18299864 | siCBX2#3 Replicate3 |
- CBX2 Knockdown in MDA-MB-231 cells
- GSE198417
- Platform: Illumina NovaSeq 6000 (Homo sapiens)
- Control: siSCR (non-silencing scrambled control siRNA)
- CBX2 Knockdown: siCBX2#3 (CBX2 targeting siRNA)
| GSM number | SRR number | Treatment |
| GSM5946425 | SRR18299888 | siSCR Replicate1 |
| GSM5946427 | SRR18299886 | siCBX2#3 Replicate1 |
| GSM5946428 | SRR18299880 | siSCR Replicate2 |
| GSM5946430 | SRR18299882 | siCBX2#3 Replicate2 |
| GSM5946431 | SRR18299883 | siSCR Replicate3 |
| GSM5946433 | SRR18299885 | siCBX2#3 Replicate3 |
2. Overall Workflow
Kallisto Quant and DESeq2 were performed on Galaxy (https://usegalaxy.org), a graphical user interface providing various data analysis tools (Afgan et al., 2018).
3. Quantifying Expression
Kallisto Quant is an alignment-free (pseudo alignment) expression estimation tool that quantifies the abundancies of RNA-seq transcripts (Bray, Pimentel, Melsted & Pachter, 2016). We first imported our datasets to Galaxy using Galaxy’s wrapper tool called “Faster Download and Extract Reads in FASTQ”, a tool specifically made for retrieving FASTQ files off SRA. Then, we performed Kallisto Quant in the paired-end format to generate TPMs (transcript per million) for each sample. All other settings were run at the default setting, and no sequence bias correction or bootstrapping was performed. The reference file used was the prebuilt index file for Homo Sapiens provided by the Pachter lab at Caltech, which developed Kallisto. The Kallisto conversion was run in batches, and each file’s tabular output was used for downstream processing. Kallisto version 0.46.2 was used in this analysis.
4. Differential Gene Expression (DEG) Analysis and Principal Component Analysis (PCA)
DESeq2 (v1.34.0) is a differential gene expression analysis tool based on the negative binomial distribution (Love, Huber & Anders, 2014). The quantified results of each sample generated by Kallisto Quant in the tabular format were imported as input files. Tximport is the default method to input the data. The .gtf file from the Ensembl homo sapiens transcriptome shared by the Pachter lab was used as the annotation file. DESeq2 was performed to identify genes that are differentially expressed with CBX2 Knockdown relative to the controls on MCF-7 and MDA-MB-231 cell lines, respectively. All other settings were run by default.
Benjamini-Hochberg procedure was used for false discovery rate (FDR) control for our differential expression gene analysis. Genes significantly differentially expressed were defined as genes with an adjusted p-value of below 0.05 and a log2 fold-change of greater than 0.58, representing a minimum of a 1.5-fold change. We also performed principal components analysis (PCA) using covariance on MCF-7 datasets based on our differential expressed genes identified.
5. Data import and summarization
We used Tximport to import TPMs, estimated counts, and transcript lengths generated from Kallisto Quant (Soneson, Love, & Robinson, 2015). Tximport also converted transcript ID into Gene ID using a table of transcript-to-gene data provided by the Kallisto developers, available on their GitHub (https://github.com/pachterlab/kallisto-transcriptome-indices/releases). All other Tximport settings were kept at defaults. The outputs were consolidated into a matrix table of both counts and abundances for all the genes covered by the reference genome used in Kallisto. Only the abundance data were used for downstream analysis. Non-DEGs, according to DESeq2, were filtered out of the dataset. The R script for Tximport is provided in Supplementary Methods.
6. Hierarchical Clustering
We performed agglomerative hierarchical clustering using the default complete linkage on our differential expressed genes from MCF-7 datasets. First, data were scaled using R’s scale() function, which converts the original data in a column of a data frame into z-scores, such that the column has a mean of 0 and a standard deviation of 1. This scaling function is a prerequisite for clustering using R’s hierarchical clustering function.
Next, a Euclidean distance matrix was generated using dist(), and the distance matrix was used as an input for hclust(), which clustered the data using the complete-linkage agglomerative method, where each element begins as an individual cluster and then is sequentially clustered until the entire dataset is consolidated. Clustering was performed along both the sample and gene axes. A heatmap of differential expression genes was also generated using the heatmap.2 function from R’s gplots package. R scripts used in this analysis are provided in Supplementary Methods.
7. Pathway analysis
An online pathway analysis tool, Reactome, was used to perform pathway analysis (Gillespie et al., 2022). The dataset from Tximport (a matrix table of abundances for each sample) was supplied to Reactome, and the pathway analysis with the down-weighting of overlapping genes (PADOG) method was used to perform the differential expression analysis between two groups by a default setting. The control groups for the MCF-7 cell type were compared to the CBX2 knockdown treatment groups.
Results
Differentially Expressed Gene (DEG) analysis using Kallisto and DESeq2
Using DESeq2 on the Galaxy platform, differentially expressed genes (DEGs) were identified for both cell lines. A summary of the total number of DEGs for each cell line is shown in figure 2, separated by up- and down-regulated genes.
For the MCF-7 cell line, there were a total of 5024 differentially expressed genes, with 2560 up-regulated genes in the treatment cell line and 2464 down-regulated genes. For the MD-MBA-231 cell line, there were no differentially expressed genes identified.
A list of the top 10 up- and down-regulated DEGs for the MCF7 cell lines treated with siCBX2-knockdown is listed below:
Principal Component Analysis (PCA)
Using significantly differentially expressed genes, principal components were determined using the Galaxy DESeq2 tool’s optional visualization function, which generates plots using R.
In figure 4, the six samples for the MCF-7 cell line were plotted using the first two principal components. The replicates for the MCF-7 cells treated with the silent control siRNA are marked pink, while the replicates for the MCF-7 cell lines treated with siCBX2 are marked blue.
PC1 accounted for 66% of the variance in samples, while PC2 accounted for 29% of the variance in samples. There is a clear separation between the control and treatment groups, with a small batch effect for replicate 3.
In figure 5, the six samples for the MDA-MB-231 cell line were plotted using the first two principal components. Again, controls are marked in pink, while samples treated with the siRNA knockdown are marked in blue.
The first principal component accounts for 50% of the sample variance, while the second principal component accounts for 41. However, no meaningful difference is observed between control and siCBX2-treated samples, and each replicate of the treatment cell line is most similar to its corresponding control replicate. At this point, the MD-MBA-231 was omitted from further analysis due to a failure to successfully identify DEGs using the computational methods detailed above.
Hierarchical clustering of MCF-7 cell line data by sample and gene
Two-way hierarchical clustering plots were generated using the heatmap.2 function from R’s gplots package.
The hierarchically clustered heatmap (figure 6) shows the distribution of up- and down-regulated genes across all six samples. The green denotes genes where expression in the associated sample is higher relative to the normalized data, which indicates up-regulation in the associated condition relative to its comparison. Red denotes down-regulation relative to its comparison. A more intense color indicates a greater degree of up- or down-regulation.
Pathway analysis for MCF-7 cell lines
Pathway analysis visualizations were generated using Reactome’s web application. In the figure below, pathways are grouped by biological function. In the figure below, the size of each node reflects the number of biological agents (including proteins, genes, and other molecules) involved with each pathway, while lines represent a connection between a pathway and its sub-pathways.
Figure 8 visualizes the same information using Reactome.org’s Reacfoam view. In this view, each “cell” represents one of Reactome’s designated pathways, and each embedded cell represents a sub-pathway.
Five of the most significant up- and down-regulated pathways are highlighted below for clarity.
Discussion
Discrepancies in gene expression analysis of the MD-MBA-231 cell line
While our workflow successfully generated a robust list of DEGs for the MCF-7 cell line using the 3 replicates per treatment condition, it identified no DEGs for the MD-MBA-231 cell line. This result was supported by a PCA analysis of the results, which showed that the majority of the variance in the dataset was between replicates rather than between experimental conditions. The expression profiles of each treated sample were most similar to its corresponding replicate of the control sample.
This result differs from the originally published paper, which was successful in identifying up-regulated genes and pathways in MD-MBA-231 cells treated with siCBX2, compared to the controls. The analysis pipeline used by the authors of the original paper differed from our pipeline in several ways:
- For alignment, Bilton et al. used HiSAT2, whereas our pipeline used Kallisto. Previous studies (Liu et al., 2022) have indicated that Kallisto may lack sensitivity for lower-expression genes that would be detected using HiSAT2. While reducing computational needs, Kallisto’s pseudo alignment method was demonstrated to sometimes produce a different number of DEGs compared to HiSAT2, which uses the burrows-wheeler transform to align reads to the reference. While our pipeline used a pre-built reference genome off of hg38 by Kallisto’s developers, HiSAT2 allows for a wider range of reference genome input formats and allows users to build their own reference genomes that may be more specific to the samples being aligned.
- For quantifying gene expression, Bilton et al. used fragments per kilobase of transcript per million fragments mapped (FPKM), whereas Kallisto measures abundance in transcripts per million (TPM). Again, prior studies have shown that TPM and FPKM, which use different approaches to normalize for transcript length, reads per sample, and sequencing biases, do not necessarily produce comparable results in downstream analysis like hierarchical even when applied to the same dataset (Dillies et al., 2013; Zhao et al., 2021).
In the mTORC1 genes highlighted in the paper, we did not observe the same statistically significant changes in expression between control and treatment groups. Two genes highlighted in the paper were TSC1 and PRKAA2. In our DESeq2 results, TSC1 had a log-fold change of only 0.067758, with an adjusted p-value of 0.99. PRKAA2 had a log-fold change of -0.05702, with an adjusted p-value of 0.99.
Due to the lack of any differentially expressed genes in the MD-MBA-231 dataset, we are not able to draw any biological conclusions about the effect of CBX2 knockdown on MD-MBA-231 cells or compare those results to the results of the MCF-7 cell line. However, our results underscore the importance of rigor and transparency in RNA-seq methodology. Using different tools for the various steps of gene expression analysis can yield results that lead to significantly different conclusions about the biology and mechanism of action of a particular treatment. This suggests that it may be important to attempt to replicate RNA-seq analysis using additional tools or provide a clearer justification for the use of certain tools over others.
Up-regulated pathways
Bilton et al. (2022) have reported that CBX2 knockdown in MDA‐MB‐231 cells upregulates the expression of the mTORC1 inhibitors TSC1 and PRKAA2. The mammalian target of rapamycin (mTOR) is known to regulate cell proliferation, autophagy, and apoptosis, and studies have shown that mTORC1 is often activated in tumors and is predominantly associated with cell growth and metabolism (Zou et al., 2020 ). TSC1 forms a protein complex that inhibits Rheb, which is a crucial activator of mTORC1 signaling (Manning & Huang, 2008). Bilton et al. (2022) observed similar effects of CBX2 knockdown in MCF-7 cells as well. In our current study, we corroborated the previous finding of upregulated TSC1 (but not PRKAA2) in MCF-7 cells that have CBX2 knocked down.
Another notable gene found in our study is TP53BP1, which is significantly upregulated in CBX2 knockdown cells compared to the control. TP53BP1 participates in the DNA repair pathway and maintains genomic stability (Mirza-Aghazadeh-Attari et al., 2019). Previous studies (Ward et al., 2005; Difilippantonio et al., 2008; Morales et al., 2006) suggest that TP53BP1 may act as a critical tumor suppressor because of its role in DNA repair.
Down-regulated pathways
Reactome analysis identified some significant down-regulated pathways that may be significantly involved in breast cancer progression. For example, the most down-regulated pathway is R-HSA-1483101, the synthesis of phosphatidylserine. The presence of phosphatidylserine on the surface of cells is an engulfment signal that encourages phagocytosis in apoptotic cells (Yu et al., 2020). Studies have suggested that PS exposure is significantly upregulated on the surface of tumor cells, and drugs that target PS have been developed and are being tested as part of drug cocktails to target various types of cancers (Chang et al., 2020).
Similarly, the second most down-regulated pathway is the macroautophagy pathway (R-HSA-1632852). Macroautophagy is the process by which portions of a cell’s cytoplasm are packaged and delivered to the lysosome for degradation. In normal cells, it recycles organelles and molecules into hydrolytic enzymes, providing new sources of energy while also removing intracellular “garbage”. In cancer cells, mutations along the autophagy pathway can have manifold effects but are particularly impactful in breast cancer stem cells (BCSC). Wolf et al. (2013) demonstrated using RNAi screens that genes involved in regulating macroautophagy are critical for maintaining pluripotency in breast cancer-derived stem cells and also improved the tumorigenic potential of breast cancer cells in tissue graft experiments. Mechanism studies have suggested that this may have been due to the involvement of macroautophagy in the IL6/STAT3 and TGFB/SMAD pathways (Niklaus et al, 2021).
A gene of significant interest in breast cancer research that displayed statistically significant differential gene expression in our dataset is GAPDH. GAPDH expression has been shown to increase breast cancer cell proliferation and tumor aggressiveness (Guo et al., 2013). Prior studies have shown GAPDH is increased in a wide range of cell cancer types, but most relevant to our study, it has been demonstrated to be significantly upregulated in prior experiments with MCF7 cell lines (Revillion et al. 2000). In our dataset, GAPDH has a statistically significant logFC of -2.441. Although the mechanism of GAPDH’s association with cancer proliferation is not yet known, researchers have hypothesized that it may be involved in cell cycle regulation because levels of GAPDH gene expression and protein quantity fluctuate in different stages of the cell cycle.
Another gene that has been associated with breast cancer proliferation that was down-regulated in our MCF-7 data set was ABCE1, a gene that has previously been shown to be overexpressed in breast cancer tissues. ABCE1 is an RNase L inhibitor, an interferon that targets and destroys intracellular DNA. It is involved in the 2-5A/RNase L pathway, which inhibits cell apoptosis via the modulation of cell metabolism. Huang et al. (2014) used siRNA to knock down ABCE1 in MCF-7 cells, which resulted in a decrease in cell viability and health, and a change in cell morphology, which they hypothesized was due to the increase in RNAse L expression as a consequence of the ABCE1 knockdown. In our data analysis, the statistically significant logFC of the ABCE1 gene in the treatment cell line was -0.6539.
The downregulation of these pathways and genes, all of which have been shown to be associated with greater progression of breast cancer when expressed at high levels, demonstrates the involvement of CBX2 in the expression of various genes that are directly involved in cancer cell proliferation, metastasis, and development. The results of our analysis affirm the claim of the original authors that CBX2 is a potential target for breast cancer therapeutics, due to its influence on various significant biological pathways in breast cancers across multiple types.
References
Afgan, E., Baker, D., Batut, B., Van Den Beek, M., Bouvier, D., Čech, M., … & Blankenberg, D. (2018). The Galaxy platform for accessible, reproducible and collaborative biomedical analyses: 2018 update. Nucleic acids research, 46(W1), W537-W544.
Bilton, L. J., Warren, C., Humphries, R. M., Kalsi, S., Waters, E., Francis, T., … & Wade, M. A. (2022). The Epigenetic Regulatory Protein CBX2 Promotes mTORC1 Signalling and Inhibits DREAM Complex Activity to Drive Breast Cancer Cell Growth. Cancers, 14(14), 3491
Bray, N. L., Pimentel, H., Melsted, P., & Pachter, L. (2016). Near-optimal probabilistic RNA-seq quantification. Nature biotechnology, 34(5), 525-527.
Chang, W., Fa, H., Xiao, D., & Wang, J. (2020). Targeting phosphatidylserine for Cancer therapy: prospects and challenges. Theranostics, 10(20), 9214–9229. https://doi.org/10.7150/thno.45125
Difilippantonio, S., Gapud, E., Wong, N., Huang, C. Y., Mahowald, G., Chen, H. T., Kruhlak, M. J., Callen, E., Livak, F., Nussenzweig, M. C., Sleckman, B. P., & Nussenzweig, A. (2008). 53BP1 facilitates long-range DNA end-joining during V(D)J recombination. Nature, 456(7221), 529–533. https://doi.org/10.1038/nature07476
Dillies, M. A., Rau, A., Aubert, J., Hennequet-Antier, C., Jeanmougin, M., Servant, N., Keime, C., Marot, G., Castel, D., Estelle, J., Guernec, G., Jagla, B., Jouneau, L., Laloë, D., Le Gall, C., Schaëffer, B., Le Crom, S., Guedj, M., Jaffrézic, F., & French StatOmique Consortium (2013). A comprehensive evaluation of normalization methods for Illumina high-throughput RNA sequencing data analysis. Briefings in bioinformatics, 14(6), 671–683. https://doi.org/10.1093/bib/bbs046
Gillespie, M., Jassal, B., Stephan, R., Milacic, M., Rothfels, K., Senff-Ribeiro, A., … & D’Eustachio, P. (2022). The reactome pathway knowledgebase 2022. Nucleic acids research, 50(D1), D687-D692.
Guo, C., Liu, S., & Sun, M. Z. (2013). Novel insight into the role of GAPDH playing in tumor. Clinical & translational oncology : official publication of the Federation of Spanish Oncology Societies and of the National Cancer Institute of Mexico, 15(3), 167–172. https://doi.org/10.1007/s12094-012-0924-x
Huang, B., Zhou, H., Lang, X., & Liu, Z. (2014). siRNA‑induced ABCE1 silencing inhibits proliferation and invasion of breast cancer cells. Molecular medicine reports, 10(4), 1685–1690. https://doi.org/10.3892/mmr.2014.2424
Huang, J., & Manning, B. D. (2008). The TSC1–TSC2 complex: a molecular switchboard controlling cell growth. Biochemical Journal, 412(2), 179-190.
Jangal, M., Lebeau, B., & Witcher, M. (2019). Beyond EZH2: is the polycomb protein CBX2 an emerging target for anti-cancer therapy? Expert opinion on therapeutic targets, 23(7), 565-578.
Liu, X., Zhao, J., Xue, L., Zhao, T., Ding, W., Han, Y., & Ye, H. (2022). A comparison of transcriptome analysis methods with reference genome. BMC genomics, 23(1), 1-15.
Love, M. I., Huber, W., & Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome biology, 15(12), 1-21.
Mirza-Aghazadeh-Attari, M., Mohammadzadeh, A., Yousefi, B., Mihanfar, A., Karimian, A., & Majidinia, M. (2019). 53BP1: A key player of DNA damage response with critical functions in cancer. DNA repair, 73, 110–119. https://doi.org/10.1016/j.dnarep.2018.11.008
Morales, J. C., Franco, S., Murphy, M. M., Bassing, C. H., Mills, K. D., Adams, M. M., … & Carpenter, P. B. (2006). 53BP1 and p53 synergize to suppress genomic instability and lymphomagenesis. Proceedings of the National Academy of Sciences, 103(9), 3310-3315.
Niklaus, N. J., Tokarchuk, I., Zbinden, M., Schläfli, A. M., Maycotte, P., & Tschan, M. P. (2021). The Multifaceted Functions of Autophagy in Breast Cancer Development and Treatment. Cells, 10(6), 1447. https://doi.org/10.3390/cells10061447
Revillion F, Pawlowski V, Hornez L et al (2000) Glyceraldehyde-3-phosphate dehydrogenase gene expression in human breast cancer. Eur J Cancer 36:1038–1042
Soneson, C., Love, M. I., & Robinson, M. D. (2015). Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Research, 4.
Ward, I. M., Difilippantonio, S., Minn, K., Mueller, M. D., Molina, J. R., Yu, X., Frisk, C. S., Ried, T., Nussenzweig, A., & Chen, J. (2005). 53BP1 cooperates with p53 and functions as a haploinsufficient tumor suppressor in mice. Molecular and cellular biology, 25(22), 10079–10086. https://doi.org/10.1128/MCB.25.22.10079-10086.2005
Wolf, J., Dewi, D. L., Fredebohm, J., Müller-Decker, K., Flechtenmacher, C., Hoheisel, J. D., & Boettcher, M. (2013). A mammosphere formation RNAi screen reveals that ATG4A promotes a breast cancer stem-like phenotype. Breast cancer research : BCR, 15(6), R109. https://doi.org/10.1186/bcr3576
Yu, M., Li, T., Li, B., Liu, Y., Wang, L., Zhang, J., … & Shi, J. (2020). Phosphatidylserine-exposing blood cells, microparticles and neutrophil extracellular traps increase procoagulant activity in patients with pancreatic cancer. Thrombosis Research, 188, 5-16.
Zhao, Y., Li, M. C., Konaté, M. M., Chen, L., Das, B., Karlovich, C., … & McShane, L. M. (2021). TPM, FPKM, or normalized counts? A comparative study of quantification measures for the analysis of RNA-seq data from the NCI patient-derived models repository. Journal of translational medicine, 19(1), 1-15.
Zou, Z., Tao, T., Li, H., & Zhu, X. (2020). mTOR signaling pathway and mTOR inhibitors in cancer: Progress and challenges. Cell & Bioscience, 10(1), 1-11.
Supplementary Methods
tximport.R
clustering.R
Analysis of Neuron and Macrophage Gene Expression in Parkinson’s Disease
Introduction and Background
Parkinson’s disease (PD) is a neurodegenerative disorder that affects aging individuals and impacts their movement, such as having uncontrollable tremors or difficulty with balance and coordination [2]. While the disorder largely impacts movement, PD can also cause depression, skin disease, and urinary problems [8]. The cause of PD is largely attributed to the death of neurons in the basal ganglia, which causes the neurons to no longer be able to produce dopamine [8]. However, recent studies have also shown that inflammation plays a central role in the disorder with resulting reactions from surrounding macrophages [1]. Therefore, it is crucial to identify differentially expressed genes and cellular pathways in affected macrophages and neurons as it may provide insight into the presence of genes that may aid in PD risk, onset, or progression. Additionally, identifying any overlapping genes and cellular pathways between affected macrophages and neurons may elucidate the link between how neurons and macrophages are related in Parkinson’s Disease.
In this study, we analyzed two separate datasets, one of affected neurons and one of affected macrophages, and determined differentially expressed genes and affected pathways of both as well as any overlapping genes and affected pathways.
Data Acquisition
Transcriptome datasets were downloaded from NCBI’s Gene Expression Omnibus (GEO) database.
| Treatment | Replicate | ID | Dataset |
|---|---|---|---|
| Control | Healthy #1 | GSM5602426 | GSE184956 |
| Healthy #2 | GSM5602428 | GSE184956 | |
| Treated | Parkinson’s #1 | GSM5602427 | GSE184956 |
| Parkinson’s #2 | GSM5602429 | GSE184956 |
| Treatment | Replicate | ID | Dataset |
|---|---|---|---|
| Control | Healthy #1 | GSM6250335 | GSE206308 |
| Healthy #2 | GSM6250334 | GSE206308 | |
| Treated | Parkinson’s #1 | GSM6250339 | GSE206308 |
| Parkinson’s #2 | GSM6250338 | GSE206308 |
Analytical Pipeline
The first part of the analytical pipeline was performed using Galaxy. The datasets were first downloaded using fasterqdump [14] before being converted into BAM files using bowtie2 [15]. The BAM files were then used in featurecounts [16], in order to normalize gene expression of the RNA-seq data. Lastly, DESeq2 was performed to determine differentially expressed genes of both data sets [17].
The second part of analysis of differentially expressed genes in order to create heatmaps was performed in R using the heatmap package [21].
The last part of the analytical pipeline was pathway analysis, and this was performed with Reactome [18].
Differentially Expressed Genes
We applied thresholds of FDR < 0.01 and log2-fold-change > 2 to find 175 differentially expressed genes between the control and Parkinson’s-associated macrophages and 28 genes between the control and Parkinson’s-associated neurons. Note: for generating the following heatmaps in Figure 1(a), we instead used a slightly higher log2-fold-change value to isolate the top 35 genes.
Comparing between the differentially expressed genes in the neuron and macrophage datasets, we found a single gene (TNC) to be up-regulated in neurons and macrophages from Parkinson’s disease. Tenascin C (TNC) encodes an extracellular matrix protein with a spatially and temporally restricted tissue distribution. This protein is homohexameric with disulfide-linked subunits, and contains multiple EGF-like and fibronectin type-III domains. It is implicated in guidance of migrating neurons as well as axons during development, synaptic plasticity, and neuronal regeneration.
Hierarchical Clustering with Heat Map
PCA treatment plots
Description of Significant Pathways in Common DEGs
Macrophage pathway analysis:
Neuron pathway analysis:
Regulation of cell adhesion pathway
Regulation of cell adhesion pathway is defined as any process that modulates the frequency, rate or extent of attachment of a cell to another cell or to the extracellular matrix [4]. In the brain, cell-cell and cell-matrix adhesions regulate the structure and function of synapses. These proteins are critical for neuronal structure, synaptic vesicle cycle, including maintenance of and transfer between vesicle pools, exocytosis, and vesicle recycling [6]. There have been studies suggesting that issues in cell adhesion may play a role in disease pathology. PD is due to synaptic dysfunction caused by genetic variations in cell adhesion pathways. Adherens Junction (AJ) proteins are involved in the cell-cell adhesions [4]. AJ proteins are involved in maintaining the blood-brain barrier (BBB) [5]. Normally, aging tends to change the BBB, however it has been noted to be more extreme in PD patients. In macrophages, this process is important because without cell adhesion pathway, it would prevent microglia to get to the sites of neuronal damage efficiently.
Regulation of biological quality
Regulation of biological quality is defined as any process that modulates a qualitative or quantitative trait of a physical quality [4]. A biological quality is a measurable attribute of an organism or part of an organism, such as size, mass, shape, color, etc. In neuronal cells, regulation of Parkinson’s disease-associated genes were done by Pumilio proteins and microRNAs. Pumilio proteins post-transcriptionally regulate gene expression through binding conserved motifs in the 3’ untranslated region (UTR) of mRNA targets. miRNAs targeted LRRK2 and SNCA, which are genes responsible for PD [3].
In macrophages, studies have concluded that inflammation plays a role in PD and can be due to resulting immune reactions from surrounding macrophages. However, microglia can act as the macrophage to protect against pathogens and regulate homeostasis in the brain [6]. M1 microglia is the first line of defense to clear the pathogens while M2 inhibits the pro-inflammatory responses [13]. Microglia can be activated to clear pathogens and maintain homeostasis [2].
Discussion
In this study, we looked to identify any overlapping genes and cellular pathways between affected macrophages and neurons which may elucidate the link between how neurons and macrophages are related in Parkinson’s Disease. To accomplish this, we had used the datasets GSE206308 and GSE184956 which are the neuron and macrophage RNA-Seq datasets used from the NCBI’s Gene Expression Omnibus (GEO) database. We have found that the gene tenascin C is present between the Lists Of Differentially Expressed Genes in the Data Sets, which past studies have agreed with. Tenascin C (TnC) is a glycoprotein highly expressed where neuron development is present [12].
Tenascin C modulates cell migration, proliferation and cellular signaling through induction of pro-inflammatory cytokines and oncogenic signaling molecules amongst other mechanisms, and is usually one of the causal roles of inflammation [20]. Previous studies have shown that both macrophage and neurons can cause neuroinflammation which plays an important role in the pathogenesis of Parkinson’s disease [2].
Drugs that target Tenascin-C expression generate hope that increased knowledge about Tenascin-C will improve the management of diseases such as Parkinson’s disease, heart diseases, and cancer. Our study confirms that one gene exists which shows that macrophages and neurons are related in Parkinson’s Disease, however further work needs to be done to help us understand the significant role that macrophages and neurons play in Parkinson’s disease
References
[1] Biju, K. C., Zhou, Q., Li, G., Imam, S. Z., Roberts, J. L., Morgan, W. W., … & Li, S. (2010). Macrophage-mediated GDNF delivery protects against dopaminergic neurodegeneration: a therapeutic strategy for Parkinson’s disease. Molecular Therapy, 18(8), 1536-1544.
[2] Blaylock, R. L. (2017). Parkinson’s disease: microglial/macrophage-induced immunoexcitotoxicity as a central mechanism of neurodegeneration. Surgical Neurology International, 8.
[3] Bogale, T. A., Faustini, G., Longhena, F., Mitola, S., Pizzi, M., & Bellucci, A. (2021). Alpha-Synuclein in the regulation of brain endothelial and perivascular cells: gaps and future perspectives. Frontiers in Immunology, 12, 611761.
[4] Carbon S, Ireland A, Mungall CJ, Shu S, Marshall B, Lewis S, AmiGO Hub, Web Presence Working Group. AmiGO: online access to ontology and annotation data. Bioinformatics. Jan 2009;25(2):288-289.
[5] Chapman, M. A. (2014). Interactions between cell adhesion and the synaptic vesicle cycle in Parkinson’s disease. Medical hypotheses, 83(2), 203-207.
[6] Edwards, Y. J., Beecham, G. W., Scott, W. K., Khuri, S., Bademci, G., Tekin, D., … & Vance, J. M. (2011). Identifying consensus disease pathways in Parkinson’s disease using an integrative systems biology approach. PLoS one, 6(2), e16917.
[7] Mamelak, M. (2018). Parkinson’s disease, the dopaminergic neuron and gammahydroxybutyrate. Neurology and Therapy, 7(1), 5-11.
[8] Parkinson’s Disease: Causes, Symptoms, And Treatments. 2022. NIH National Institute on Aging. https://www.nia.nih.gov/health/parkinsons-disease
[9] Reactome. n.d. Pathway Browser Analysis. https://reactome.org/PathwayBrowser/#/DTAB=AN&ANALYSIS=MjAyMjExMDYxOTMzMjRfNDc4NjM%253D&FILTER=resource:UNIPROT.
[10] Reactome. n.d. Pathway Browser Analysis. https://reactome.org/PathwayBrowser/#/DTAB=AN&ANALYSIS=MjAyMjExMDYxOTM3MzNfNDc4Njg%253D&FILTER=resource:UNIPROT.
[11] Su, R., & Zhou, T. (2021). Alpha-synuclein induced immune cells activation and associated therapy in Parkinson’s disease. Frontiers in Aging Neuroscience, 13.
[12] Tucić, M., Stamenković, V., & Andjus, P. (2021). The Extracellular Matrix Glycoprotein Tenascin C and Adult Neurogenesis. Frontiers in Cell and Developmental Biology, 9, 674199.
[13] Yan, A., Zhang, Y., Lin, J., Song, L., Wang, X., & Liu, Z. (2018). Partial depletion of peripheral M1 macrophages reverses motor deficits in MPTP-treated mouse by suppressing neuroinflammation and dopaminergic neurodegeneration. Frontiers in Aging Neuroscience, 10, 160.
[14] Leinonen, R., Sugawara, H., & Shumway, M. (2010). The Sequence Read Archive. Nucleic Acids Research, 39(Database), D19â D21. https://doi.org/10.1093/nar/gkq1019
[15] Langmead, B., Trapnell, C., Pop, M., & Salzberg, S. L. (2009). Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biology, 10(3), R25. https://doi.org/10.1186/gb-2009-10-3-r25
[16] Liao, Y., Smyth, G. K., & Shi, W. (2013). featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics, 30(7), 923â930. https://doi.org/10.1093/bioinformatics/btt656
[17] Love, M. I., Huber, W., & Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15(12). https://doi.org/10.1186/s13059-014-0550-8
[18] Griss J, Viteri G, Sidiropoulos K, Nguyen V, Fabregat A, Hermjakob H. ReactomeGSA – Efficient Multi-Omics Comparative Pathway Analysis. Mol Cell Proteomics. 2020 Sep 9. doi: 10.1074/mcp. PubMed PMID: 32907876.
[19]Blaylock RL. Parkinson’s disease: Microglial/macrophage-induced immunoexcitotoxicity as a central mechanism of neurodegeneration. Surg Neurol Int. 2017 Apr 26;8:65. doi: 10.4103/sni.sni_441_16. PMID: 28540131; PMCID: PMC5421223.
[20]Midwood KS, Orend G. The role of tenascin-C in tissue injury and tumorigenesis. J Cell Commun Signal. 2009 Dec;3(3-4):287-310. doi: 10.1007/s12079-009-0075-1. Epub 2009 Oct 17. PMID: 19838819; PMCID: PMC2778592.
[21] Perry M (2022). heatmaps: Flexible Heatmaps for Functional Genomics and Sequence Features. R package version 1.22.0.
Biology Computes 