{ "cells": [ { "cell_type": "markdown", "id": "a76ae282-3b11-4246-8292-a9276267832d", "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_histonemarkchipseq.ipynb) [![Open In nbviewer](https://img.shields.io/badge/View%20in-nbviewer-orange)](https://nbviewer.jupyter.org/github/lucascamillomd/pyaging/blob/main/tutorials/tutorial_histonemarkchipseq.ipynb)" ] }, { "cell_type": "markdown", "id": "d444a24e-6a98-4db1-8688-7f3f80ed2876", "metadata": {}, "source": [ "# Bulk histone mark ChIP-Seq" ] }, { "cell_type": "markdown", "id": "186154f3-1c8d-4284-a5a4-01f28d4db533", "metadata": {}, "source": [ "This tutorial is a brief guide for the implementation of the seven histone-mark-specific clocks and the pan-histone-mark clock developed ourselves. See the [published article](https://doi.org/10.1126/sciadv.adk9373)." ] }, { "cell_type": "markdown", "id": "270379c1-9159-4677-92fa-10b08aa9f703", "metadata": {}, "source": [ "We just need two packages for this tutorial." ] }, { "cell_type": "code", "execution_count": null, "id": "f816b2da-312d-4b51-b2bb-a64fd0ff5b4e", "metadata": {}, "outputs": [], "source": [ "!pip install pybigwig -q" ] }, { "cell_type": "code", "execution_count": 1, "id": "dd281360-7e16-45d9-ae2b-8f8f3fff809d", "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "import pyaging as pya" ] }, { "cell_type": "markdown", "id": "b6893601-615e-449b-829b-c144276f402f", "metadata": {}, "source": [ "## Download and load example data" ] }, { "cell_type": "markdown", "id": "fd3e80a9-5361-40f0-bf3e-6f6057181594", "metadata": {}, "source": [ "Let's download an example of H3K4me3 ChIP-Seq bigWig file from the ENCODE project." ] }, { "cell_type": "code", "execution_count": 2, "id": "85c15bf3-6cf1-4f71-abf2-d0d7ee81b86b", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "|-----> 🏗️ Starting download_example_data function\n", "|-----------> Downloading data to pyaging_data/ENCFF386QWG.bigWig\n", "|-----------> in progress: 100.0000%\n", "|-----> 🎉 Done! [2075.2977s]\n" ] } ], "source": [ "pya.data.download_example_data('ENCFF386QWG')" ] }, { "cell_type": "markdown", "id": "3880246a-471e-4f75-bd2f-ed2623458a48", "metadata": {}, "source": [ "To exemplify that multiple bigWigs can be turned into a df object at once, let's just repeat the file path." ] }, { "cell_type": "code", "execution_count": 3, "id": "f65f5cc7-4c42-45a5-a04e-83e0520eccff", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "|-----> 🏗️ Starting bigwig_to_df function\n", "|-----> ⚙️ Load Ensembl genome metadata started\n", "|-----------> Downloading data to pyaging_data/Ensembl-105-EnsDb-for-Homo-sapiens-genes.csv\n", "|-----------> in progress: 100.0000%\n", "|-----> ✅ Load Ensembl genome metadata finished [20.4846s]\n", "|-----> ⚙️ Processing bigWig files started\n", "|-----------> Processing file: pyaging_data/ENCFF386QWG.bigWig\n", "|-----------> in progress: 100.0000%\n", "|-----------> Processing file: pyaging_data/ENCFF386QWG.bigWig\n", "|-----------> in progress: 100.0000%\n", "|-----> ✅ Processing bigWig files finished [12.4227s]\n", "|-----> 🎉 Done! [44.8066s]\n" ] } ], "source": [ "df = pya.pp.bigwig_to_df(['pyaging_data/ENCFF386QWG.bigWig', 'pyaging_data/ENCFF386QWG.bigWig'])" ] }, { "cell_type": "code", "execution_count": 4, "id": "1a24e0a5-f97f-4f01-95a7-dd96246d9eb2", "metadata": {}, "outputs": [], "source": [ "df.index = ['sample1', 'sample2'] # just to avoid an annoying anndata warning that samples have same names" ] }, { "cell_type": "code", "execution_count": 5, "id": "769858ac-9d6d-43f8-9c53-0f4a88c5484c", "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", "
ENSG00000223972ENSG00000227232ENSG00000278267ENSG00000243485ENSG00000284332ENSG00000237613ENSG00000268020ENSG00000240361ENSG00000186092ENSG00000238009...ENSG00000237801ENSG00000237040ENSG00000124333ENSG00000228410ENSG00000223484ENSG00000124334ENSG00000270726ENSG00000185203ENSG00000182484ENSG00000227159
sample10.0286160.0304150.0277830.0286160.0286160.0286160.0441710.0364740.0307840.03181...0.0344350.0068221.4131190.0294240.1400050.0497860.0692960.3321260.0285960.028616
sample20.0286160.0304150.0277830.0286160.0286160.0286160.0441710.0364740.0307840.03181...0.0344350.0068221.4131190.0294240.1400050.0497860.0692960.3321260.0285960.028616
\n", "

2 rows × 62241 columns

\n", "
" ], "text/plain": [ " ENSG00000223972 ENSG00000227232 ENSG00000278267 ENSG00000243485 \\\n", "sample1 0.028616 0.030415 0.027783 0.028616 \n", "sample2 0.028616 0.030415 0.027783 0.028616 \n", "\n", " ENSG00000284332 ENSG00000237613 ENSG00000268020 ENSG00000240361 \\\n", "sample1 0.028616 0.028616 0.044171 0.036474 \n", "sample2 0.028616 0.028616 0.044171 0.036474 \n", "\n", " ENSG00000186092 ENSG00000238009 ... ENSG00000237801 \\\n", "sample1 0.030784 0.03181 ... 0.034435 \n", "sample2 0.030784 0.03181 ... 0.034435 \n", "\n", " ENSG00000237040 ENSG00000124333 ENSG00000228410 ENSG00000223484 \\\n", "sample1 0.006822 1.413119 0.029424 0.140005 \n", "sample2 0.006822 1.413119 0.029424 0.140005 \n", "\n", " ENSG00000124334 ENSG00000270726 ENSG00000185203 ENSG00000182484 \\\n", "sample1 0.049786 0.069296 0.332126 0.028596 \n", "sample2 0.049786 0.069296 0.332126 0.028596 \n", "\n", " ENSG00000227159 \n", "sample1 0.028616 \n", "sample2 0.028616 \n", "\n", "[2 rows x 62241 columns]" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df.head()" ] }, { "cell_type": "markdown", "id": "e303dc0f-9e77-4524-9c04-90540e9ee75d", "metadata": {}, "source": [ "## Convert data to AnnData object" ] }, { "cell_type": "markdown", "id": "ae8e44bc-67fc-4508-9623-faea44301fa8", "metadata": {}, "source": [ "AnnData objects are highly flexible and are thus our preferred method of organizing data for age prediction." ] }, { "cell_type": "code", "execution_count": 6, "id": "6c167be6-1bd3-407c-ae12-771739189c3c", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "|-----> 🏗️ Starting df_to_adata function\n", "|-----> ⚙️ Create anndata object started\n", "|-----> ✅ Create anndata object finished [0.0282s]\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 2 observations\n", "|-----------> There are 62241 features\n", "|-----------> Total missing values: 0\n", "|-----------> Percentage of missing values: 0.00%\n", "|-----> ✅ Log data statistics finished [0.0010s]\n", "|-----> ⚙️ Impute missing values started\n", "|-----------> No missing values found. No imputation necessary\n", "|-----> ✅ Impute missing values finished [0.0016s]\n", "|-----> 🎉 Done! [0.0333s]\n" ] } ], "source": [ "adata = pya.preprocess.df_to_adata(df)" ] }, { "cell_type": "markdown", "id": "3f82813b-3db2-4570-9e4c-3dce08dc5108", "metadata": {}, "source": [ "Note that the original DataFrame is stored in `X_original` under layers. This is what the `adata` object looks like:" ] }, { "cell_type": "code", "execution_count": 7, "id": "641a61a6-46fc-4d47-b176-eb39524ce94f", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "AnnData object with n_obs × n_vars = 2 × 62241\n", " var: 'percent_na'\n", " layers: 'X_original'" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adata" ] }, { "cell_type": "markdown", "id": "c72aa719-efd3-4094-90f5-bffcaea76a34", "metadata": {}, "source": [ "## Predict age" ] }, { "cell_type": "markdown", "id": "aff9395b-4954-4148-9cbb-6681e7217cf3", "metadata": {}, "source": [ "We can either predict one clock at once or all at the same time. For convenience, let's simply input a few clocks of interest at once. The function is invariant to the capitalization of the clock name. " ] }, { "cell_type": "code", "execution_count": 8, "id": "c02455b4-06dd-44c2-b4b3-a2bb434eae7d", "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.0031s]\n", "|-----> 🕒 Processing clock: camilloh3k4me3\n", "|-----------> ⚙️ Load clock started\n", "|-----------------> Downloading data to pyaging_data/camilloh3k4me3.pt\n", "|-----------------> in progress: 100.0000%\n", "|-----------> ✅ Load clock finished [5.0816s]\n", "|-----------> ⚙️ Check features in adata started\n", "|-----------------> All features are present in adata.var_names.\n", "|-----------------> Added prepared input matrix to adata.obsm[X_camilloh3k4me3]\n", "|-----------> ✅ Check features in adata finished [0.0076s]\n", "|-----------> ⚙️ Predict ages with model started\n", "|-----------------> There is no preprocessing necessary\n", "|-----------------> There is no postprocessing necessary\n", "|-----------------> in progress: 100.0000%\n", "|-----------> ✅ Predict ages with model finished [0.0250s]\n", "|-----------> ⚙️ Add predicted ages and clock metadata to adata started\n", "|-----------> ✅ Add predicted ages and clock metadata to adata finished [0.0013s]\n", "|-----> 🕒 Processing clock: camilloh3k9me3\n", "|-----------> ⚙️ Load clock started\n", "|-----------------> Downloading data to pyaging_data/camilloh3k9me3.pt\n", "|-----------------> in progress: 100.0000%\n", "|-----------> ✅ Load clock finished [12.2932s]\n", "|-----------> ⚙️ Check features in adata started\n", "|-----------------> All features are present in adata.var_names.\n", "|-----------------> Added prepared input matrix to adata.obsm[X_camilloh3k9me3]\n", "|-----------> ✅ Check features in adata finished [0.0187s]\n", "|-----------> ⚙️ Predict ages with model started\n", "|-----------------> There is no preprocessing necessary\n", "|-----------------> There is no postprocessing necessary\n", "|-----------------> in progress: 100.0000%\n", "|-----------> ✅ Predict ages with model finished [0.0030s]\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: camillopanhistone\n", "|-----------> ⚙️ Load clock started\n", "|-----------------> Downloading data to pyaging_data/camillopanhistone.pt\n", "|-----------------> in progress: 100.0000%\n", "|-----------> ✅ Load clock finished [222.1535s]\n", "|-----------> ⚙️ Check features in adata started\n", "|-----------------> All features are present in adata.var_names.\n", "|-----------------> Added prepared input matrix to adata.obsm[X_camillopanhistone]\n", "|-----------> ✅ Check features in adata finished [0.0201s]\n", "|-----------> ⚙️ Predict ages with model started\n", "|-----------------> There is no preprocessing necessary\n", "|-----------------> There is no postprocessing necessary\n", "|-----------------> in progress: 100.0000%\n", "|-----------> ✅ Predict ages with model finished [0.0109s]\n", "|-----------> ⚙️ Add predicted ages and clock metadata to adata started\n", "|-----------> ✅ Add predicted ages and clock metadata to adata finished [0.0008s]\n", "|-----> 🎉 Done! [239.8831s]\n" ] } ], "source": [ "pya.pred.predict_age(adata, ['CamilloH3K4me3', 'CamilloH3K9me3', 'CamilloPanHistone'])" ] }, { "cell_type": "code", "execution_count": 9, "id": "f64fb182-937b-4f67-b58e-5fffb0e2fad0", "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", "
camilloh3k4me3camilloh3k9me3camillopanhistone
sample153.99854444.322954.021884
sample253.99854444.322954.021884
\n", "
" ], "text/plain": [ " camilloh3k4me3 camilloh3k9me3 camillopanhistone\n", "sample1 53.998544 44.3229 54.021884\n", "sample2 53.998544 44.3229 54.021884" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adata.obs.head()" ] }, { "cell_type": "markdown", "id": "bbaa2243-e380-4020-bf04-f7aa7da83cd4", "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": 10, "id": "e8dd3457-8983-41a4-aaab-41563b91a866", "metadata": {}, "outputs": [], "source": [ "pya.data.download_example_data('ENCFF386QWG', verbose=False)\n", "df = pya.pp.bigwig_to_df(['pyaging_data/ENCFF386QWG.bigWig', 'pyaging_data/ENCFF386QWG.bigWig'], verbose=False)\n", "df.index = ['sample1', 'sample2']\n", "adata = pya.preprocess.df_to_adata(df, verbose=False)\n", "pya.pred.predict_age(adata, ['CamilloH3K4me3', 'CamilloH3K9me3', 'CamilloPanHistone'], verbose=False)" ] }, { "cell_type": "code", "execution_count": 11, "id": "8192ab67-a1cc-4728-8ca0-f81a56940fbf", "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", "
camilloh3k4me3camilloh3k9me3camillopanhistone
sample153.99854444.322954.021884
sample253.99854444.322954.021884
\n", "
" ], "text/plain": [ " camilloh3k4me3 camilloh3k9me3 camillopanhistone\n", "sample1 53.998544 44.3229 54.021884\n", "sample2 53.998544 44.3229 54.021884" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adata.obs.head()" ] }, { "cell_type": "markdown", "id": "9832aa0b-99a8-4938-a2a2-5e9b484a3353", "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": 12, "id": "a4b22bf1-116f-456f-82d2-58b300f863f1", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "AnnData object with n_obs × n_vars = 2 × 62241\n", " obs: 'camilloh3k4me3', 'camilloh3k9me3', 'camillopanhistone'\n", " var: 'percent_na'\n", " uns: 'camilloh3k4me3_percent_na', 'camilloh3k4me3_missing_features', 'camilloh3k4me3_metadata', 'camilloh3k9me3_percent_na', 'camilloh3k9me3_missing_features', 'camilloh3k9me3_metadata', 'camillopanhistone_percent_na', 'camillopanhistone_missing_features', 'camillopanhistone_metadata'\n", " layers: 'X_original'" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adata" ] }, { "cell_type": "markdown", "id": "c08ff758-675c-4136-9fb8-c19f0e05fefd", "metadata": {}, "source": [ "## Get citation" ] }, { "cell_type": "markdown", "id": "8407c418-6251-4b08-9d29-166f9a4339d2", "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": 13, "id": "2946393e-a199-46ba-a9dd-80bc8fa88787", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'clock_name': 'camilloh3k4me3',\n", " 'data_type': 'histone mark',\n", " 'species': 'Homo sapiens',\n", " 'year': 2023,\n", " 'approved_by_author': '✅',\n", " 'citation': 'de Lima Camillo, Lucas Paulo, et al. \"Histone mark age of human tissues and cell types.\" Science Advances 11.1 (2025): eadk9373.',\n", " 'doi': 'https://doi.org/10.1126/sciadv.adk9373',\n", " 'notes': None,\n", " 'version': None}" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adata.uns['camilloh3k4me3_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 }