{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# OP2: 🧠 Biologically-informed dimensionality reduction\n\nMethods like PCA (Principal Component Analysis) and ICA (Independent Component Analysis) are commonly used to reduce the dimensionality of data, creating new features in the process. However, this often comes at the cost of reduced interpretability and the potential neglect of prior biological knowledge. Here we will see how can we use differential expression data to estimate the potential activity of **Transcription Factors (TF)** in the dataset, using latest scientific results on this topic.\n\nTFs are proteins that play a pivotal role in controlling gene expression 🧬. By binding to specific DNA sequences, they act as regulators, determining whether a particular gene is turned on or off. If we can understand the activity of these transcription factors, we gain significant insights into changes in gene expression. Essentially, knowing how these transcription factors operate gives us a clearer picture of why certain genes may change upon perturbation with compounds. By knowing which genes each TF control, we can infer their activity from changes in gene expression, effectively reducing the number of genes by summarizing their activity through TFs.\n\n[*Decoupler*<sup>1</sup>](https://github.com/saezlab/decoupler-py) is a very active python package containing different statistical methods to extract biological activities from omics data within a unified framework. To do so, it needs some biological knowledge about the interactions between TFs and genes. Decoupler uses [*Omnipath*<sup>2</sup>](https://omnipathdb.org/) to obtain data from biological databases. We will use the most recent prior knowledge on TF - Gene interactions called [*CollecTRI*<sup>3</sup>](https://github.com/saezlab/CollecTRI). Combining prior knowledge networks with biological data is a powerful way to reduce the dimensionality of the data while incorporating biological knowledge in the process.\n\nWe will see how can we use it to perform statistical tests and obtain rich and biologically informed features.\n\n**Tool is available on GitHub:** https://github.com/saezlab/decoupler-py","metadata":{}},{"cell_type":"code","source":"%%capture\n!pip install decoupler==1.4.0 # for enrichment\n!pip install omnipath==1.0.7  # for biological prior knowledge networks\n!pip install adjustText       # only for plotting","metadata":{"execution":{"iopub.status.busy":"2023-09-29T08:55:25.128710Z","iopub.execute_input":"2023-09-29T08:55:25.129067Z","iopub.status.idle":"2023-09-29T08:56:01.597673Z","shell.execute_reply.started":"2023-09-29T08:55:25.129040Z","shell.execute_reply":"2023-09-29T08:56:01.596137Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import decoupler as dc # decoupler uses numba, takes some time to jit\nimport pandas as pd","metadata":{"execution":{"iopub.status.busy":"2023-09-29T08:56:01.600213Z","iopub.execute_input":"2023-09-29T08:56:01.600698Z","iopub.status.idle":"2023-09-29T08:58:00.401834Z","shell.execute_reply.started":"2023-09-29T08:56:01.600657Z","shell.execute_reply":"2023-09-29T08:58:00.400629Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = pd.read_parquet(\n    \"/kaggle/input/open-problems-single-cell-perturbations/de_train.parquet\"\n).set_index([\"cell_type\", \"sm_name\"]).iloc[:, 3:]\ndf.index = df.index.map('@'.join)\ndf.head()","metadata":{"execution":{"iopub.status.busy":"2023-09-29T08:58:00.402868Z","iopub.execute_input":"2023-09-29T08:58:00.403526Z","iopub.status.idle":"2023-09-29T08:58:02.947453Z","shell.execute_reply.started":"2023-09-29T08:58:00.403492Z","shell.execute_reply":"2023-09-29T08:58:02.946706Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"try:\n    # Query the DB to obtain the most updated version of the TF - Gene interactions for the enrichment analysis.\n    df_collectri = dc.get_collectri(organism='human', split_complexes=False)\nexcept:\n    df_collectri = pd.read_csv(\"https://github.com/pablormier/omnipath-static/raw/main/op/collectri-26.09.2023.zip\")\ndf_collectri","metadata":{"execution":{"iopub.status.busy":"2023-09-29T08:58:02.949247Z","iopub.execute_input":"2023-09-29T08:58:02.949843Z","iopub.status.idle":"2023-09-29T08:58:05.432010Z","shell.execute_reply.started":"2023-09-29T08:58:02.949816Z","shell.execute_reply":"2023-09-29T08:58:05.430460Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The dataframe above contains the known interactions between TFs and genes that they regulate. Now we will use that knowledge with Decoupler to perform statistical tests to estimate the activity of TFs based on the expression pattern of the genes in control for each TF. There are many enrichment tests in Decoupler, let's use for example ULM, described in the following image ([from the docs](https://decoupler-py.readthedocs.io/en/latest/notebooks/bulk.html)):\n\n\n<img src=\"https://github.com/saezlab/decoupler-py/blob/main/docs/source/ulm.png?raw=true\" width=\"800px\"/>","metadata":{}},{"cell_type":"code","source":"tf_activities, p_values = dc.run_ulm(\n    mat=df,\n    net=df_collectri,\n    source='source',\n    target='target',\n    weight='weight',\n    verbose=True,\n)\ntf_activities","metadata":{"execution":{"iopub.status.busy":"2023-09-29T08:58:05.434719Z","iopub.execute_input":"2023-09-29T08:58:05.435247Z","iopub.status.idle":"2023-09-29T08:58:51.382514Z","shell.execute_reply.started":"2023-09-29T08:58:05.435198Z","shell.execute_reply":"2023-09-29T08:58:51.381456Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As a result, we summarized the activity of thousands of genes into 632 transcription factors.","metadata":{}},{"cell_type":"code","source":"tf_activities.shape","metadata":{"execution":{"iopub.status.busy":"2023-09-29T08:58:51.383897Z","iopub.execute_input":"2023-09-29T08:58:51.384850Z","iopub.status.idle":"2023-09-29T08:58:51.392617Z","shell.execute_reply.started":"2023-09-29T08:58:51.384811Z","shell.execute_reply":"2023-09-29T08:58:51.391655Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Visualizing TF activities\n\nA very interesting feature for biological interpretation is the use of the different plots available in Decoupler to have a look at the data. Let's see for example which are the top 25 enriched TFs in NK cells treated with Dabrafenib","metadata":{}},{"cell_type":"code","source":"dc.plot_barplot(tf_activities, 'NK cells@Dabrafenib', top = 25, vertical = True)","metadata":{"execution":{"iopub.status.busy":"2023-09-29T08:58:51.393777Z","iopub.execute_input":"2023-09-29T08:58:51.394811Z","iopub.status.idle":"2023-09-29T08:58:52.691536Z","shell.execute_reply.started":"2023-09-29T08:58:51.394775Z","shell.execute_reply":"2023-09-29T08:58:52.690345Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The plot shows a strong inhibition in activity of the transcription factor MYC. This comes from the fact that many genes under the control of MYC are downregulated after treatment. This downregulation could potentially dampen NK cell proliferation, metabolism, and effector functions, given MYC's role in promoting these cellular activities.\n\nBut why MYC appears downregulated? we can explore the values of the target genes\n\n> NOTE: Remeber that in the dataset, what we have is -log(p_values) of the differential expression analysis. This does not give us information about the effect size. For better interpretation, these plots should be used with log fold changes","metadata":{}},{"cell_type":"code","source":"df.loc[[\"NK cells@Dabrafenib\"], :]","metadata":{"execution":{"iopub.status.busy":"2023-09-29T08:58:52.692852Z","iopub.execute_input":"2023-09-29T08:58:52.693232Z","iopub.status.idle":"2023-09-29T08:58:52.718998Z","shell.execute_reply.started":"2023-09-29T08:58:52.693201Z","shell.execute_reply":"2023-09-29T08:58:52.717761Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dc.plot_targets(\n    df.loc[[\"NK cells@Dabrafenib\"],:].T, \n    stat = 'NK cells@Dabrafenib',\n    source_name = 'MYC', \n    net = df_collectri, \n    top = 15\n)","metadata":{"execution":{"iopub.status.busy":"2023-09-29T08:58:52.720980Z","iopub.execute_input":"2023-09-29T08:58:52.721333Z","iopub.status.idle":"2023-09-29T08:58:56.211531Z","shell.execute_reply.started":"2023-09-29T08:58:52.721282Z","shell.execute_reply":"2023-09-29T08:58:56.210289Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Since Decoupler returns two dataframes, one with effect sizes and the other with p-values, we can use for example the p-values to generate features that will look very similar to the original ones used in the competition:","metadata":{}},{"cell_type":"code","source":"import numpy as np\n\ndf_tfpval = -np.log(p_values+1e-10) * np.sign(tf_activities)\ndf_tfpval","metadata":{"execution":{"iopub.status.busy":"2023-09-29T08:58:56.214013Z","iopub.execute_input":"2023-09-29T08:58:56.214342Z","iopub.status.idle":"2023-09-29T08:58:56.254277Z","shell.execute_reply.started":"2023-09-29T08:58:56.214314Z","shell.execute_reply":"2023-09-29T08:58:56.253139Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Reducing the dimensionality even further\n\nIf the idea is to use only a bunch of TFs as features for prediction, a good strategy can be to pick only the top variable genes across conditions:","metadata":{}},{"cell_type":"code","source":"# Using the TF activities\ntf_activities.std(axis=0).sort_values(ascending=False).head(10)","metadata":{"execution":{"iopub.status.busy":"2023-09-29T08:58:56.255446Z","iopub.execute_input":"2023-09-29T08:58:56.255757Z","iopub.status.idle":"2023-09-29T08:58:56.269116Z","shell.execute_reply.started":"2023-09-29T08:58:56.255731Z","shell.execute_reply":"2023-09-29T08:58:56.268037Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Using the -log(p vals)\ndf_tfpval.std(axis=0).sort_values(ascending=False).head(10)","metadata":{"execution":{"iopub.status.busy":"2023-09-29T08:58:56.270524Z","iopub.execute_input":"2023-09-29T08:58:56.270882Z","iopub.status.idle":"2023-09-29T08:58:56.286689Z","shell.execute_reply.started":"2023-09-29T08:58:56.270854Z","shell.execute_reply":"2023-09-29T08:58:56.285539Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"There are multiple methods to calculate TF activities and P-values, as outlined in Decoupler's documentation https://decoupler-py.readthedocs.io/en/latest/api.html. The different assumptions made by these statistical enrichment tests can provide diverse insights for models. These insights can be leveraged, for instance, by using ensembles of models with features computed using different statistical tests.\n\n## References\n\n1. Badia-i-Mompel, Pau, et al. \"decoupleR: ensemble of computational methods to infer biological activities from omics data.\" Bioinformatics Advances 2.1 (2022): vbac016.\n2. Türei, Dénes, Tamás Korcsmáros, and Julio Saez-Rodriguez. \"OmniPath: guidelines and gateway for literature-curated signaling pathway resources.\" Nature methods 13.12 (2016): 966-967.\n3. Mueller-Dott, Sophia, et al. \"Expanding the coverage of regulons from high-confidence prior knowledge for accurate estimation of transcription factor activities.\" bioRxiv (2023): 2023-03.","metadata":{}}]}