{ "cells": [ { "cell_type": "markdown", "id": "82321cbb-e1b9-49f3-b826-32c4fafd96f4", "metadata": {}, "source": [ "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/lucascamillomd/pyaging/blob/main/tutorials/tutorial_atacseq.ipynb) [![Open In nbviewer](https://img.shields.io/badge/View%20in-nbviewer-orange)](https://nbviewer.jupyter.org/github/lucascamillomd/pyaging/blob/main/tutorials/tutorial_atacseq.ipynb)" ] }, { "cell_type": "markdown", "id": "cf837dcf-de19-46f7-9c81-c063a45b14b6", "metadata": {}, "source": [ "# Bulk ATAC-Seq" ] }, { "cell_type": "markdown", "id": "133e64f6-a0d2-4a70-84a0-33a4c2bea725", "metadata": {}, "source": [ "This tutorial is a brief guide for the implementation of the two ATAC clocks developed by Morandini et al. Link to [paper](https://link.springer.com/article/10.1007/s11357-023-00986-0)." ] }, { "cell_type": "markdown", "id": "a0469770-c9af-4d69-a055-bf3c312286db", "metadata": {}, "source": [ "We just need two packages for this tutorial." ] }, { "cell_type": "code", "execution_count": 1, "id": "e64f6bb6-7c95-4b9e-b37a-4ae811dc088d", "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "import pyaging as pya " ] }, { "cell_type": "markdown", "id": "5c367e47-ceca-4c73-ac14-b6a9a781c66a", "metadata": {}, "source": [ "## Download and load example data" ] }, { "cell_type": "markdown", "id": "52b6e9c2-b7b7-4655-ab9c-84afe3185b78", "metadata": {}, "source": [ "If you have your own ATAC-Seq data, please follow the recommendations in the Ocampo paper. Specifically, one needs to count the number of reads for each of the peak regions from the paper (file [here](https://static-content.springer.com/esm/art%3A10.1007%2Fs11357-023-00986-0/MediaObjects/11357_2023_986_MOESM9_ESM.tsv)). This can be done through the code found on their [GitHub](https://github.com/SunScript0/ATAC-clock/blob/main/pipeline_atac/02_peakset_and_counts.sh) using featureCounts." ] }, { "cell_type": "markdown", "id": "76f0a1a6-ad70-4a40-8fd0-63de208c7ad5", "metadata": {}, "source": [ "For testing purposes, let's download an example of input for the ATAC clocks. For instructions on how to go from raw sequencing reads to the data table, please refer to the paper. " ] }, { "cell_type": "code", "execution_count": 2, "id": "992db2fb-2b58-4f8b-92ac-f760df3758dc", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "|-----> 🏗️ Starting download_example_data function\n", "|-----------> Data found in pyaging_data/GSE193140.pkl\n", "|-----> 🎉 Done! [0.4942s]\n" ] } ], "source": [ "pya.data.download_example_data('GSE193140')" ] }, { "cell_type": "code", "execution_count": 3, "id": "078b2b96-7317-4f84-a8c4-16276fb76137", "metadata": {}, "outputs": [], "source": [ "df = pd.read_pickle('pyaging_data/GSE193140.pkl')" ] }, { "cell_type": "code", "execution_count": 4, "id": "cd1a8e9b-53a5-4cfc-ac6d-22d9ed09784a", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
chr1:817100-817691chr1:826742-828191chr1:841908-843021chr1:844055-844921chr1:857908-859108chr1:869571-870271chr1:898378-899076chr1:904303-905702chr1:906675-907111chr1:912617-913368...chrY:21073148-21074236chrY:21174455-21175401chrY:21177324-21177828chrY:21180682-21181317chrY:21239902-21241040chrY:21248553-21249961chrY:21256824-21257260chrY:21259823-21260874chrY:22086084-22086722chrY:22499696-22500344
CR_124182265215119843271455637...6210465319020502122
CR_12296268827254010971378616712...11133125270372918912
CR_12113727854246691297863835124...0000000000
CR_120169281929354613732093130110...798331514750183214
CR_1192053005184537102533113824136...1518171257257877
\n", "

5 rows × 80400 columns

\n", "
" ], "text/plain": [ " chr1:817100-817691 chr1:826742-828191 chr1:841908-843021 \\\n", "CR_124 182 2652 15 \n", "CR_122 96 2688 27 \n", "CR_121 137 2785 42 \n", "CR_120 169 2819 29 \n", "CR_119 205 3005 18 \n", "\n", " chr1:844055-844921 chr1:857908-859108 chr1:869571-870271 \\\n", "CR_124 11 9 843 \n", "CR_122 25 40 1097 \n", "CR_121 46 69 1297 \n", "CR_120 35 46 1373 \n", "CR_119 45 37 1025 \n", "\n", " chr1:898378-899076 chr1:904303-905702 chr1:906675-907111 \\\n", "CR_124 2 714 556 \n", "CR_122 13 786 167 \n", "CR_121 8 638 351 \n", "CR_120 20 931 301 \n", "CR_119 33 1138 241 \n", "\n", " chr1:912617-913368 ... chrY:21073148-21074236 \\\n", "CR_124 37 ... 62 \n", "CR_122 12 ... 11 \n", "CR_121 24 ... 0 \n", "CR_120 10 ... 7 \n", "CR_119 36 ... 15 \n", "\n", " chrY:21174455-21175401 chrY:21177324-21177828 \\\n", "CR_124 104 65 \n", "CR_122 13 31 \n", "CR_121 0 0 \n", "CR_120 9 8 \n", "CR_119 18 17 \n", "\n", " chrY:21180682-21181317 chrY:21239902-21241040 \\\n", "CR_124 31 90 \n", "CR_122 25 270 \n", "CR_121 0 0 \n", "CR_120 33 151 \n", "CR_119 12 57 \n", "\n", " chrY:21248553-21249961 chrY:21256824-21257260 \\\n", "CR_124 20 50 \n", "CR_122 37 29 \n", "CR_121 0 0 \n", "CR_120 47 50 \n", "CR_119 25 7 \n", "\n", " chrY:21259823-21260874 chrY:22086084-22086722 chrY:22499696-22500344 \n", "CR_124 21 2 2 \n", "CR_122 18 9 12 \n", "CR_121 0 0 0 \n", "CR_120 18 32 14 \n", "CR_119 8 7 7 \n", "\n", "[5 rows x 80400 columns]" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df.head()" ] }, { "cell_type": "markdown", "id": "251495e7-082f-45ae-841c-a2dd86a3cb15", "metadata": {}, "source": [ "## Convert data to AnnData object" ] }, { "cell_type": "markdown", "id": "73ba54e0-4292-4d85-b208-e56e267d6797", "metadata": {}, "source": [ "AnnData objects are highly flexible and are thus our preferred method of organizing data for age prediction." ] }, { "cell_type": "code", "execution_count": 5, "id": "4dab5019-9f54-4e32-be19-abbb1c71a2d6", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "|-----> 🏗️ Starting df_to_adata function\n", "|-----> ⚙️ Create anndata object started\n", "|-----> ✅ Create anndata object finished [0.0289s]\n", "|-----> ⚙️ Add metadata to anndata started\n", "|-----------? No metadata provided. Leaving adata.obs empty\n", "|-----> ⚠️ Add metadata to anndata finished [0.0004s]\n", "|-----> ⚙️ Log data statistics started\n", "|-----------> There are 157 observations\n", "|-----------> There are 80400 features\n", "|-----------> Total missing values: 0\n", "|-----------> Percentage of missing values: 0.00%\n", "|-----> ✅ Log data statistics finished [0.0049s]\n", "|-----> ⚙️ Impute missing values started\n", "|-----------> No missing values found. No imputation necessary\n", "|-----> ✅ Impute missing values finished [0.0053s]\n", "|-----> 🎉 Done! [0.0419s]\n" ] } ], "source": [ "adata = pya.preprocess.df_to_adata(df)" ] }, { "cell_type": "markdown", "id": "5042e04f-17c0-4eb2-8c5d-2c2fc5d6d2d6", "metadata": {}, "source": [ "Note that the original DataFrame is stored in `X_original` under layers. is This is what the `adata` object looks like:" ] }, { "cell_type": "code", "execution_count": 6, "id": "503da312-2256-4e67-9747-107f5c4587ec", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "AnnData object with n_obs × n_vars = 157 × 80400\n", " var: 'percent_na'\n", " layers: 'X_original'" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adata" ] }, { "cell_type": "markdown", "id": "c072990d-0f54-49b3-bb7a-7bbd13301e2a", "metadata": {}, "source": [ "## Predict age" ] }, { "cell_type": "markdown", "id": "5fe08978-f1ba-49b3-a0d8-52df4b6efb4e", "metadata": {}, "source": [ "We can either predict one clock at once or all at the same time. For convenience, let's simply input all two clocks of interest at once. The function is invariant to the capitalization of the clock name. " ] }, { "cell_type": "code", "execution_count": 7, "id": "96e008fe-9f8c-45fb-8dc6-6a39f1ecb7ac", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "|-----> 🏗️ Starting predict_age function\n", "|-----> ⚙️ Set PyTorch device started\n", "|-----------> Using device: cpu\n", "|-----> ✅ Set PyTorch device finished [0.0006s]\n", "|-----> 🕒 Processing clock: ocampoatac1\n", "|-----------> ⚙️ Load clock started\n", "|-----------------> Data found in pyaging_data/ocampoatac1.pt\n", "|-----------> ✅ Load clock finished [0.5113s]\n", "|-----------> ⚙️ Check features in adata started\n", "|-----------------> All features are present in adata.var_names.\n", "|-----------------> Added prepared input matrix to adata.obsm[X_ocampoatac1]\n", "|-----------> ✅ Check features in adata finished [3.8480s]\n", "|-----------> ⚙️ Predict ages with model started\n", "|-----------------> The preprocessing method is tpm_norm_log1p\n", "|-----------------> There is no postprocessing necessary\n", "|-----------------> in progress: 100.0000%\n", "|-----------> ✅ Predict ages with model finished [0.1635s]\n", "|-----------> ⚙️ Add predicted ages and clock metadata to adata started\n", "|-----------> ✅ Add predicted ages and clock metadata to adata finished [0.0007s]\n", "|-----> 🕒 Processing clock: ocampoatac2\n", "|-----------> ⚙️ Load clock started\n", "|-----------------> Data found in pyaging_data/ocampoatac2.pt\n", "|-----------> ✅ Load clock finished [0.4514s]\n", "|-----------> ⚙️ Check features in adata started\n", "|-----------------> All features are present in adata.var_names.\n", "|-----------------> Added prepared input matrix to adata.obsm[X_ocampoatac2]\n", "|-----------> ✅ Check features in adata finished [4.9598s]\n", "|-----------> ⚙️ Predict ages with model started\n", "|-----------------> The preprocessing method is tpm_norm_log1p\n", "|-----------------> There is no postprocessing necessary\n", "|-----------------> in progress: 100.0000%\n", "|-----------> ✅ Predict ages with model finished [0.0690s]\n", "|-----------> ⚙️ Add predicted ages and clock metadata to adata started\n", "|-----------> ✅ Add predicted ages and clock metadata to adata finished [0.0007s]\n", "|-----> 🎉 Done! [10.1175s]\n" ] } ], "source": [ "pya.pred.predict_age(adata, ['OcampoATAC1', 'OcampoATAC2'])" ] }, { "cell_type": "code", "execution_count": 8, "id": "b83a10e2-7984-4427-9e85-05329a16feb3", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
ocampoatac1ocampoatac2
CR_12429.52712428.114206
CR_12239.00309740.061162
CR_12140.71600843.095199
CR_12032.38037233.033456
CR_11936.44071138.301516
\n", "
" ], "text/plain": [ " ocampoatac1 ocampoatac2\n", "CR_124 29.527124 28.114206\n", "CR_122 39.003097 40.061162\n", "CR_121 40.716008 43.095199\n", "CR_120 32.380372 33.033456\n", "CR_119 36.440711 38.301516" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adata.obs.head()" ] }, { "cell_type": "markdown", "id": "a2cfa1e5-a7f4-4157-8c66-4afcc2323ef7", "metadata": {}, "source": [ "Having so much information printed can be overwhelming, particularly when running several clocks at once. In such cases, just set verbose to False." ] }, { "cell_type": "code", "execution_count": 9, "id": "055761d9-7e22-49f3-a1db-31c3ed3749ba", "metadata": {}, "outputs": [], "source": [ "pya.data.download_example_data('GSE193140', verbose=False)\n", "df = pd.read_pickle('pyaging_data/GSE193140.pkl')\n", "adata = pya.preprocess.df_to_adata(df, verbose=False)\n", "pya.pred.predict_age(adata, ['OcampoATAC1', 'OcampoATAC2'], verbose=False)" ] }, { "cell_type": "code", "execution_count": 10, "id": "fdd9d6c2-7f0a-4f96-a095-4a492ed73f8d", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
ocampoatac1ocampoatac2
CR_12429.52712428.114206
CR_12239.00309740.061162
CR_12140.71600843.095199
CR_12032.38037233.033456
CR_11936.44071138.301516
\n", "
" ], "text/plain": [ " ocampoatac1 ocampoatac2\n", "CR_124 29.527124 28.114206\n", "CR_122 39.003097 40.061162\n", "CR_121 40.716008 43.095199\n", "CR_120 32.380372 33.033456\n", "CR_119 36.440711 38.301516" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adata.obs.head()" ] }, { "cell_type": "markdown", "id": "4bad3df8-f868-4cf5-be74-00ffd02c18f5", "metadata": {}, "source": [ "After age prediction, the clocks are added to `adata.obs`. Moreover, the percent of missing values for each clock and other metadata are included in `adata.uns`." ] }, { "cell_type": "code", "execution_count": 11, "id": "0d13fb55-8a12-4d28-83e9-ec7c9fbbe30c", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "AnnData object with n_obs × n_vars = 157 × 80400\n", " obs: 'ocampoatac1', 'ocampoatac2'\n", " var: 'percent_na'\n", " uns: 'ocampoatac1_percent_na', 'ocampoatac1_missing_features', 'ocampoatac1_metadata', 'ocampoatac2_percent_na', 'ocampoatac2_missing_features', 'ocampoatac2_metadata'\n", " layers: 'X_original'" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adata" ] }, { "cell_type": "markdown", "id": "a4e7ad8d-44ae-4ced-a626-f9e3b2d04114", "metadata": {}, "source": [ "## Get citation" ] }, { "cell_type": "markdown", "id": "d7d089b0-5433-47a9-b031-bc4504c6b55d", "metadata": {}, "source": [ "The doi, citation, and some metadata are automatically added to the AnnData object under `adata.uns[CLOCKNAME_metadata]`." ] }, { "cell_type": "code", "execution_count": 12, "id": "6b368506-55d1-4b74-be61-817bcf575ade", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'clock_name': 'ocampoatac1',\n", " 'data_type': 'atac',\n", " 'species': 'Homo sapiens',\n", " 'year': 2023,\n", " 'approved_by_author': '⌛',\n", " 'citation': 'Morandini, Francesco, et al. \"ATAC-clock: An aging clock based on chromatin accessibility.\" GeroScience (2023): 1-18.',\n", " 'doi': 'https://doi.org/10.1007/s11357-023-00986-0',\n", " 'notes': None,\n", " 'version': None}" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adata.uns['ocampoatac1_metadata']" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.17" } }, "nbformat": 4, "nbformat_minor": 5 }