{ "cells": [ { "cell_type": "markdown", "id": "4e690b3c-4dec-450e-a7f8-f63987e60cdb", "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_dnam.ipynb) [![Open In nbviewer](https://img.shields.io/badge/View%20in-nbviewer-orange)](https://nbviewer.jupyter.org/github/lucascamillomd/pyaging/blob/main/tutorials/tutorial_dnam.ipynb)" ] }, { "cell_type": "markdown", "id": "62e1ac68-927d-4ca8-a2ab-bd99a7ee52ab", "metadata": {}, "source": [ "# RRBS DNA methylation" ] }, { "cell_type": "markdown", "id": "9552602a-777c-42a5-900a-41c85096c3d8", "metadata": {}, "source": [ "This tutorial focuses on predicting age from Mus musculus reduced-representation bisulfite sequencing (RRBS) data. There are a few clocks available that were trained on RRBS data. Moreover, it is possible to use Horvath's mammalian clocks by converting the genomic location to the probes in the Horvath methylation array." ] }, { "cell_type": "code", "execution_count": 1, "id": "dc347ae0-41ae-46ac-ba50-08cacd4c9241", "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "import pyaging as pya\n", "import os\n", "import numpy as np" ] }, { "cell_type": "markdown", "id": "5e620ab9-4837-4a7a-83f1-726be9c9f7bf", "metadata": {}, "source": [ "## Download and load example data" ] }, { "cell_type": "markdown", "id": "5ff79235-46fb-4c59-a629-1f479f9f13a3", "metadata": {}, "source": [ "Let's download the publicly available dataset GSE130735 with RRBS samples from mouse. Given it is RRBS, there are millions of CpG sites." ] }, { "cell_type": "code", "execution_count": 2, "id": "332d96c6-6b12-4cd1-b216-c32ce21673b9", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "|-----> 🏗️ Starting download_example_data function\n", "|-----------> Data found in pyaging_data/GSE130735_subset.pkl\n", "|-----> 🎉 Done! [0.5425s]\n" ] } ], "source": [ "pya.data.download_example_data('GSE130735')" ] }, { "cell_type": "code", "execution_count": 3, "id": "2a6fcd90-8ded-40d5-a606-e32e21816ebf", "metadata": {}, "outputs": [], "source": [ "df = pd.read_pickle('pyaging_data/GSE130735_subset.pkl')" ] }, { "cell_type": "markdown", "id": "df8ea1a3-313f-42bc-aeef-ec5349975b80", "metadata": {}, "source": [ "It is important to note that the features for RRBS clocks are the genomic coordinates in the format below." ] }, { "cell_type": "code", "execution_count": 4, "id": "c4d8245b-8d04-4ae0-945d-3aed4956a3bb", "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:3020814chr1:3020842chr1:3020877chr1:3020891chr1:3020945chr1:3020971chr1:3020987chr1:3021012chr1:3037802chr1:3037820...chrY:1825397chrY:4682362chrY:32122892chrY:85867071chrY:85867083chrY:85867117chrY:85867137chrY:85867139chrY:85867178chrY:88224179
GSM37526310.6090.250.4080.1890.0680.3730.5710.2520.3330.158...NaNNaNNaNNaNNaNNaNNaNNaNNaNNaN
GSM3752625NaNNaN0.9730.9840.9120.9150.9870.9740.9910.932...NaNNaNNaNNaNNaNNaNNaNNaNNaNNaN
GSM3752634NaNNaN0.5260.1310.0000.0380.4690.7690.7720.146...NaNNaNNaNNaNNaNNaNNaNNaNNaNNaN
GSM37526200.9310.920.9880.9490.8970.9210.9070.9581.0000.867...NaNNaNNaNNaNNaNNaNNaNNaNNaNNaN
GSM3752622NaNNaN0.2050.3820.0910.1320.1740.2270.1080.053...NaNNaNNaNNaNNaNNaNNaNNaNNaNNaN
\n", "

5 rows × 1778324 columns

\n", "
" ], "text/plain": [ " chr1:3020814 chr1:3020842 chr1:3020877 chr1:3020891 \\\n", "GSM3752631 0.609 0.25 0.408 0.189 \n", "GSM3752625 NaN NaN 0.973 0.984 \n", "GSM3752634 NaN NaN 0.526 0.131 \n", "GSM3752620 0.931 0.92 0.988 0.949 \n", "GSM3752622 NaN NaN 0.205 0.382 \n", "\n", " chr1:3020945 chr1:3020971 chr1:3020987 chr1:3021012 \\\n", "GSM3752631 0.068 0.373 0.571 0.252 \n", "GSM3752625 0.912 0.915 0.987 0.974 \n", "GSM3752634 0.000 0.038 0.469 0.769 \n", "GSM3752620 0.897 0.921 0.907 0.958 \n", "GSM3752622 0.091 0.132 0.174 0.227 \n", "\n", " chr1:3037802 chr1:3037820 ... chrY:1825397 chrY:4682362 \\\n", "GSM3752631 0.333 0.158 ... NaN NaN \n", "GSM3752625 0.991 0.932 ... NaN NaN \n", "GSM3752634 0.772 0.146 ... NaN NaN \n", "GSM3752620 1.000 0.867 ... NaN NaN \n", "GSM3752622 0.108 0.053 ... NaN NaN \n", "\n", " chrY:32122892 chrY:85867071 chrY:85867083 chrY:85867117 \\\n", "GSM3752631 NaN NaN NaN NaN \n", "GSM3752625 NaN NaN NaN NaN \n", "GSM3752634 NaN NaN NaN NaN \n", "GSM3752620 NaN NaN NaN NaN \n", "GSM3752622 NaN NaN NaN NaN \n", "\n", " chrY:85867137 chrY:85867139 chrY:85867178 chrY:88224179 \n", "GSM3752631 NaN NaN NaN NaN \n", "GSM3752625 NaN NaN NaN NaN \n", "GSM3752634 NaN NaN NaN NaN \n", "GSM3752620 NaN NaN NaN NaN \n", "GSM3752622 NaN NaN NaN NaN \n", "\n", "[5 rows x 1778324 columns]" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df.head()" ] }, { "cell_type": "markdown", "id": "ea7c44d7-73c3-4cd7-844d-bab34aa2dcee", "metadata": {}, "source": [ "## Convert data to AnnData object" ] }, { "cell_type": "markdown", "id": "04f2758f-fb8c-4a52-983a-29ec826dba6c", "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": "682f7e04-22a6-4561-b389-c8f336f19862", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "|-----> 🏗️ Starting df_to_adata function\n", "|-----> ⚙️ Create anndata object started\n", "|-----> ✅ Create anndata object finished [0.9882s]\n", "|-----> ⚙️ Add metadata to anndata started\n", "|-----------? No metadata provided. Leaving adata.obs empty\n", "|-----> ⚠️ Add metadata to anndata finished [0.0006s]\n", "|-----> ⚙️ Log data statistics started\n", "|-----------> There are 14 observations\n", "|-----------> There are 1778324 features\n", "|-----------> Total missing values: 6322346\n", "|-----------> Percentage of missing values: 25.39%\n", "|-----> ✅ Log data statistics finished [0.0205s]\n", "|-----> ⚙️ Impute missing values started\n", "|-----------> Imputing missing values using mean strategy\n", "|-----> ✅ Impute missing values finished [0.4631s]\n", "|-----> ⚙️ Add imputer strategy to adata.uns started\n", "|-----> ✅ Add imputer strategy to adata.uns finished [0.0087s]\n", "|-----> 🎉 Done! [1.4897s]\n" ] } ], "source": [ "adata = pya.pp.df_to_adata(df, imputer_strategy='mean') # knn might be a bit slow" ] }, { "cell_type": "markdown", "id": "7349164c-f28b-4222-bf41-6f80d8b79c3b", "metadata": {}, "source": [ "This is what the `adata` object looks like:" ] }, { "cell_type": "markdown", "id": "4b5ff1ef-e724-407a-b6d4-9907558f21ba", "metadata": {}, "source": [ "## Predict age with RRBS clocks" ] }, { "cell_type": "markdown", "id": "eb197ded-91dd-4319-8dbb-a635d09c8367", "metadata": {}, "source": [ "We can either predict one clock at once or all at the same time. For convenience, let's simply input all four available mammalian clocks at once. The function is invariant to the capitalization of the clock name." ] }, { "cell_type": "code", "execution_count": 6, "id": "cdd18ee5-9af1-404e-80f6-42a83685273e", "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.0033s]\n", "|-----> 🕒 Processing clock: thompson\n", "|-----------> ⚙️ Load clock started\n", "|-----------------> Data found in pyaging_data/thompson.pt\n", "|-----------> ✅ Load clock finished [0.5324s]\n", "|-----------> ⚙️ Check features in adata started\n", "|-----------------? 1 out of 582 features (0.17%) are missing: ['chr4:91376687'], etc.\n", "|-----------------> Filling missing features entirely with 0\n", "|-----------------> Added prepared input matrix to adata.obsm[X_thompson]\n", "|-----------> ⚠️ Check features in adata finished [0.0654s]\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.0013s]\n", "|-----------> ⚙️ Add predicted ages and clock metadata to adata started\n", "|-----------> ✅ Add predicted ages and clock metadata to adata finished [0.0008s]\n", "|-----> 🕒 Processing clock: meer\n", "|-----------> ⚙️ Load clock started\n", "|-----------------> Data found in pyaging_data/meer.pt\n", "|-----------> ✅ Load clock finished [0.4402s]\n", "|-----------> ⚙️ Check features in adata started\n", "|-----------------? 225 out of 435 features (51.72%) are missing: ['chr10:111559529', 'chr10:115250413', 'chr10:127620127'], etc.\n", "|-----------------> Filling missing features entirely with 0\n", "|-----------------> Added prepared input matrix to adata.obsm[X_meer]\n", "|-----------> ⚠️ Check features in adata finished [0.0412s]\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.0010s]\n", "|-----------> ⚙️ Add predicted ages and clock metadata to adata started\n", "|-----------> ✅ Add predicted ages and clock metadata to adata finished [0.0006s]\n", "|-----> 🕒 Processing clock: petkovich\n", "|-----------> ⚙️ Load clock started\n", "|-----------------> Data found in pyaging_data/petkovich.pt\n", "|-----------> ✅ Load clock finished [0.5167s]\n", "|-----------> ⚙️ Check features in adata started\n", "|-----------------? 58 out of 90 features (64.44%) are missing: ['chr19:23893237', 'chr18:45589182', 'chr16:10502162'], etc.\n", "|-----------------> Filling missing features entirely with 0\n", "|-----------------> Added prepared input matrix to adata.obsm[X_petkovich]\n", "|-----------> ⚠️ Check features in adata finished [0.0161s]\n", "|-----------> ⚙️ Predict ages with model started\n", "|-----------------> There is no preprocessing necessary\n", "|-----------------> The postprocessing method is petkovich\n", "|-----------------> in progress: 100.0000%\n", "|-----------> ✅ Predict ages with model finished [0.0033s]\n", "|-----------> ⚙️ Add predicted ages and clock metadata to adata started\n", "|-----------> ✅ Add predicted ages and clock metadata to adata finished [0.0020s]\n", "|-----> 🕒 Processing clock: stubbs\n", "|-----------> ⚙️ Load clock started\n", "|-----------------> Data found in pyaging_data/stubbs.pt\n", "|-----------> ✅ Load clock finished [0.4679s]\n", "|-----------> ⚙️ Check features in adata started\n", "|-----------------? 8889 out of 17992 features (49.41%) are missing: ['chr1:10038066', 'chr1:106173313', 'chr1:106759301'], etc.\n", "|-----------------> Using reference feature values for stubbs\n", "|-----------------> Added prepared input matrix to adata.obsm[X_stubbs]\n", "|-----------> ⚠️ Check features in adata finished [0.8672s]\n", "|-----------> ⚙️ Predict ages with model started\n", "|-----------------> The preprocessing method is quantile_normalization_and_scale_with_gold_standard\n", "|-----------------> The postprocessing method is stubbs\n", "|-----------------> in progress: 100.0000%\n", "|-----------> ✅ Predict ages with model finished [0.0263s]\n", "|-----------> ⚙️ Add predicted ages and clock metadata to adata started\n", "|-----------> ✅ Add predicted ages and clock metadata to adata finished [0.0014s]\n", "|-----> 🎉 Done! [3.2757s]\n" ] } ], "source": [ "pya.pred.predict_age(adata, ['Thompson', 'Meer', 'Petkovich', 'Stubbs'])" ] }, { "cell_type": "markdown", "id": "3978afec-40a0-4e1f-8ff6-1d048da8a894", "metadata": {}, "source": [ "All of the age predictions are in unit of months." ] }, { "cell_type": "code", "execution_count": 7, "id": "b04572f9-23dd-4eb1-8e84-16a9b25c2d6a", "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", "
thompsonmeerpetkovichstubbs
GSM375263119.6341137.3151838.0751770.957770
GSM3752625-1.4104610.0282212.953822-0.074265
GSM375263461.05878321.3221789.6404891.389193
GSM3752620-2.6638151.6119473.019351-0.092710
GSM375262220.5941147.5921457.1047660.667168
\n", "
" ], "text/plain": [ " thompson meer petkovich stubbs\n", "GSM3752631 19.634113 7.315183 8.075177 0.957770\n", "GSM3752625 -1.410461 0.028221 2.953822 -0.074265\n", "GSM3752634 61.058783 21.322178 9.640489 1.389193\n", "GSM3752620 -2.663815 1.611947 3.019351 -0.092710\n", "GSM3752622 20.594114 7.592145 7.104766 0.667168" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adata.obs.head()" ] }, { "cell_type": "code", "execution_count": 8, "id": "4bb259c5-2cba-4dc1-b123-2387a5bb7749", "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", "
thompsonmeerpetkovichstubbs
GSM375263119.6341137.3151838.0751770.957770
GSM3752625-1.4104610.0282212.953822-0.074265
GSM375263461.05878321.3221789.6404891.389193
GSM3752620-2.6638151.6119473.019351-0.092710
GSM375262220.5941147.5921457.1047660.667168
\n", "
" ], "text/plain": [ " thompson meer petkovich stubbs\n", "GSM3752631 19.634113 7.315183 8.075177 0.957770\n", "GSM3752625 -1.410461 0.028221 2.953822 -0.074265\n", "GSM3752634 61.058783 21.322178 9.640489 1.389193\n", "GSM3752620 -2.663815 1.611947 3.019351 -0.092710\n", "GSM3752622 20.594114 7.592145 7.104766 0.667168" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adata.obs.head()" ] }, { "cell_type": "markdown", "id": "8519affc-ffc8-4904-ad7a-bd6a6d6458cf", "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": "18b44cfa-36d5-49c9-badf-7ba9e189bbc0", "metadata": {}, "outputs": [], "source": [ "pya.data.download_example_data('GSE130735', verbose=False)\n", "df = pd.read_pickle('pyaging_data/GSE130735_subset.pkl')\n", "adata = pya.preprocess.df_to_adata(df, imputer_strategy='mean', verbose=False)\n", "pya.pred.predict_age(adata, ['Thompson', 'Meer', 'Petkovich', 'Stubbs'], verbose=False)" ] }, { "cell_type": "code", "execution_count": 10, "id": "a2520978-b693-474f-88cf-91bcde1a5d95", "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", "
thompsonmeerpetkovichstubbs
GSM375263119.6341137.3151838.0751770.957770
GSM3752625-1.4104610.0282212.953822-0.074265
GSM375263461.05878321.3221789.6404891.389193
GSM3752620-2.6638151.6119473.019351-0.092710
GSM375262220.5941147.5921457.1047660.667168
\n", "
" ], "text/plain": [ " thompson meer petkovich stubbs\n", "GSM3752631 19.634113 7.315183 8.075177 0.957770\n", "GSM3752625 -1.410461 0.028221 2.953822 -0.074265\n", "GSM3752634 61.058783 21.322178 9.640489 1.389193\n", "GSM3752620 -2.663815 1.611947 3.019351 -0.092710\n", "GSM3752622 20.594114 7.592145 7.104766 0.667168" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adata.obs.head()" ] }, { "cell_type": "markdown", "id": "33119798-f1b3-4c4c-9f18-e4e4b7ca21e8", "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": "92cfc16e-71ff-4767-9c75-04e52455eb6c", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "AnnData object with n_obs × n_vars = 14 × 1778324\n", " obs: 'thompson', 'meer', 'petkovich', 'stubbs'\n", " var: 'percent_na'\n", " uns: 'imputer_strategy', 'thompson_percent_na', 'thompson_missing_features', 'thompson_metadata', 'meer_percent_na', 'meer_missing_features', 'meer_metadata', 'petkovich_percent_na', 'petkovich_missing_features', 'petkovich_metadata', 'stubbs_percent_na', 'stubbs_missing_features', 'stubbs_metadata'\n", " layers: 'X_original', 'X_imputed'" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adata" ] }, { "cell_type": "markdown", "id": "569b1ef8-dd55-45e4-8792-e919dc207808", "metadata": {}, "source": [ "## Predict age with mammalian clocks" ] }, { "cell_type": "markdown", "id": "b9266699-881f-41c7-9a13-c38307527bfa", "metadata": {}, "source": [ "We can predict age by converting the genomic locations directly into the probes from Horvath's methylation array. " ] }, { "cell_type": "code", "execution_count": 12, "id": "a5f71a02-cbc0-4cd0-a123-e7226307ae84", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "os.system('git clone https://github.com/shorvath/MammalianMethylationConsortium.git')\n", "\n", "# Let's read the manifest from the mammalian consortium\n", "annotation_df = pd.read_csv('MammalianMethylationConsortium/Annotations, Amin Haghani/Mammals/Mus_musculus.grcm38.100.HorvathMammalMethylChip40.v1.csv', index_col=0)\n", "annotation_df = annotation_df[~annotation_df.seqnames.isna()]\n", "mm_genomic_locations = 'chr' + annotation_df['seqnames'].astype(str) + ':' + annotation_df['CGstart'].astype(int).astype(str)\n", "mm_genomic_locations = mm_genomic_locations.tolist()\n", "mammalian_probes = annotation_df['CGid'].tolist()\n", "mm_loc_to_probe = dict(zip(mm_genomic_locations, mammalian_probes))\n", "\n", "# Let's get the previous RRBS dataset and filter only for the genomic locations in the manifest file\n", "df_columns_set = set(df.columns)\n", "mm_loc_to_probe_set = set(mm_loc_to_probe.keys())\n", "common_columns = df_columns_set.intersection(mm_loc_to_probe_set)\n", "df_converted = df[list(common_columns)].copy()\n", "\n", "# Then, convert the genomic location to the probe name\n", "df_converted.columns = [mm_loc_to_probe[col] for col in df_converted.columns]\n", "\n", "# Let's clean the GitHub\n", "os.system('rm -r MammalianMethylationConsortium')" ] }, { "cell_type": "code", "execution_count": 13, "id": "6989844c-736d-4637-acfc-8f8f7ae58108", "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", "
cg05347424cg26718996cg07727941cg16852837cg12870762cg26080798cg02899039cg12839061cg05267150cg13170453...cg02179016cg20836420cg18831685cg08992395cg13679010cg12982463cg17146242cg13649253cg07588415cg14814195
GSM3752631NaNNaN0.0000.00.0150.0000.005NaN0.0230.000...NaN0.0000.028NaN0.0000.0000.0180.0000.021NaN
GSM37526250.938NaN0.000NaN0.0000.0000.000NaN0.596NaN...NaN0.000NaNNaN0.0000.8950.2270.1560.025NaN
GSM37526340.125NaN0.000NaN0.6270.0170.033NaN0.745NaN...NaN0.495NaNNaN0.0140.2780.5190.7860.012NaN
GSM37526200.769NaN0.0910.00.0700.0060.012NaN0.6070.092...NaN0.0100.054NaN0.0000.9330.2770.1480.000NaN
GSM3752622NaNNaN0.000NaN0.0000.0000.000NaN0.052NaN...NaN0.000NaNNaN0.0000.0000.0640.0220.000NaN
\n", "

5 rows × 5149 columns

\n", "
" ], "text/plain": [ " cg05347424 cg26718996 cg07727941 cg16852837 cg12870762 \\\n", "GSM3752631 NaN NaN 0.000 0.0 0.015 \n", "GSM3752625 0.938 NaN 0.000 NaN 0.000 \n", "GSM3752634 0.125 NaN 0.000 NaN 0.627 \n", "GSM3752620 0.769 NaN 0.091 0.0 0.070 \n", "GSM3752622 NaN NaN 0.000 NaN 0.000 \n", "\n", " cg26080798 cg02899039 cg12839061 cg05267150 cg13170453 ... \\\n", "GSM3752631 0.000 0.005 NaN 0.023 0.000 ... \n", "GSM3752625 0.000 0.000 NaN 0.596 NaN ... \n", "GSM3752634 0.017 0.033 NaN 0.745 NaN ... \n", "GSM3752620 0.006 0.012 NaN 0.607 0.092 ... \n", "GSM3752622 0.000 0.000 NaN 0.052 NaN ... \n", "\n", " cg02179016 cg20836420 cg18831685 cg08992395 cg13679010 \\\n", "GSM3752631 NaN 0.000 0.028 NaN 0.000 \n", "GSM3752625 NaN 0.000 NaN NaN 0.000 \n", "GSM3752634 NaN 0.495 NaN NaN 0.014 \n", "GSM3752620 NaN 0.010 0.054 NaN 0.000 \n", "GSM3752622 NaN 0.000 NaN NaN 0.000 \n", "\n", " cg12982463 cg17146242 cg13649253 cg07588415 cg14814195 \n", "GSM3752631 0.000 0.018 0.000 0.021 NaN \n", "GSM3752625 0.895 0.227 0.156 0.025 NaN \n", "GSM3752634 0.278 0.519 0.786 0.012 NaN \n", "GSM3752620 0.933 0.277 0.148 0.000 NaN \n", "GSM3752622 0.000 0.064 0.022 0.000 NaN \n", "\n", "[5 rows x 5149 columns]" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df_converted.head()" ] }, { "cell_type": "markdown", "id": "e46ca085-026e-4e02-b316-97b880125507", "metadata": {}, "source": [ "Now we can finally put the dataframe into pyaging after defining the species as Mus musculus." ] }, { "cell_type": "code", "execution_count": 14, "id": "10bbdf6e-a63a-4c18-bc97-e0872fb9895f", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "|-----> 🏗️ Starting df_to_adata function\n", "|-----> ⚙️ Create anndata object started\n", "|-----> ✅ Create anndata object finished [0.0057s]\n", "|-----> ⚙️ Add metadata to anndata started\n", "|-----------? No metadata provided. Leaving adata.obs empty\n", "|-----> ⚠️ Add metadata to anndata finished [0.0006s]\n", "|-----> ⚙️ Log data statistics started\n", "|-----------> There are 14 observations\n", "|-----------> There are 5150 features\n", "|-----------> Total missing values: 17862\n", "|-----------> Percentage of missing values: 24.77%\n", "|-----> ✅ Log data statistics finished [0.0013s]\n", "|-----> ⚙️ Impute missing values started\n", "|-----------> Imputing missing values using mean strategy\n", "|-----> ✅ Impute missing values finished [0.0060s]\n", "|-----> ⚙️ Add imputer strategy to adata.uns started\n", "|-----> ✅ Add imputer strategy to adata.uns finished [0.0004s]\n", "|-----> 🎉 Done! [0.0174s]\n" ] } ], "source": [ "df_converted['Mus musculus'] = 1\n", "adata_mammalian = pya.pp.df_to_adata(df_converted, imputer_strategy='mean')" ] }, { "cell_type": "markdown", "id": "7edae9ce-a0bc-4b78-8575-23745714b42b", "metadata": {}, "source": [ "Let's use these five mammalian predictors." ] }, { "cell_type": "code", "execution_count": 15, "id": "56dc5c27-e793-4343-85e2-9ce30a365d64", "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.0013s]\n", "|-----> 🕒 Processing clock: mammalian1\n", "|-----------> ⚙️ Load clock started\n", "|-----------------> Data found in pyaging_data/mammalian1.pt\n", "|-----------> ✅ Load clock finished [0.4780s]\n", "|-----------> ⚙️ Check features in adata started\n", "|-----------------? 274 out of 335 features (81.79%) are missing: ['cg00249943', 'cg00250826', 'cg00292639'], etc.\n", "|-----------------> Filling missing features entirely with 0\n", "|-----------------> Added prepared input matrix to adata.obsm[X_mammalian1]\n", "|-----------> ⚠️ Check features in adata finished [0.0173s]\n", "|-----------> ⚙️ Predict ages with model started\n", "|-----------------> There is no preprocessing necessary\n", "|-----------------> The postprocessing method is anti_logp2\n", "|-----------------> in progress: 100.0000%\n", "|-----------> ✅ Predict ages with model finished [0.0083s]\n", "|-----------> ⚙️ Add predicted ages and clock metadata to adata started\n", "|-----------> ✅ Add predicted ages and clock metadata to adata finished [0.0017s]\n", "|-----> 🕒 Processing clock: mammalian2\n", "|-----------> ⚙️ Load clock started\n", "|-----------------> Data found in pyaging_data/mammalian2.pt\n", "|-----------> ✅ Load clock finished [0.4544s]\n", "|-----------> ⚙️ Check features in adata started\n", "|-----------------? 2406 out of 2572 features (93.55%) are missing: ['cg00020468', 'cg00096922', 'cg00098422'], etc.\n", "|-----------------> Using reference feature values for mammalian2\n", "|-----------------> Added prepared input matrix to adata.obsm[X_mammalian2]\n", "|-----------> ⚠️ Check features in adata finished [0.0407s]\n", "|-----------> ⚙️ Predict ages with model started\n", "|-----------------> There is no preprocessing necessary\n", "|-----------------> The postprocessing method is mammalian2\n", "|-----------------> in progress: 100.0000%\n", "|-----------> ✅ Predict ages with model finished [0.0196s]\n", "|-----------> ⚙️ Add predicted ages and clock metadata to adata started\n", "|-----------> ✅ Add predicted ages and clock metadata to adata finished [0.0006s]\n", "|-----> 🕒 Processing clock: mammalian3\n", "|-----------> ⚙️ Load clock started\n", "|-----------------> Data found in pyaging_data/mammalian3.pt\n", "|-----------> ✅ Load clock finished [0.5081s]\n", "|-----------> ⚙️ Check features in adata started\n", "|-----------------? 2299 out of 2467 features (93.19%) are missing: ['cg00101675', 'cg06259996', 'cg15168457'], etc.\n", "|-----------------> Using reference feature values for mammalian3\n", "|-----------------> Added prepared input matrix to adata.obsm[X_mammalian3]\n", "|-----------> ⚠️ Check features in adata finished [0.0222s]\n", "|-----------> ⚙️ Predict ages with model started\n", "|-----------------> There is no preprocessing necessary\n", "|-----------------> The postprocessing method is mammalian3\n", "|-----------------> in progress: 100.0000%\n", "|-----------> ✅ Predict ages with model finished [0.0095s]\n", "|-----------> ⚙️ Add predicted ages and clock metadata to adata started\n", "|-----------> ✅ Add predicted ages and clock metadata to adata finished [0.0006s]\n", "|-----> 🕒 Processing clock: mammalianlifespan\n", "|-----------> ⚙️ Load clock started\n", "|-----------------> Data found in pyaging_data/mammalianlifespan.pt\n", "|-----------> ✅ Load clock finished [0.4420s]\n", "|-----------> ⚙️ Check features in adata started\n", "|-----------------? 133 out of 152 features (87.50%) are missing: ['cg00039845', 'cg00300233', 'cg00810217'], etc.\n", "|-----------------> Using reference feature values for mammalianlifespan\n", "|-----------------> Added prepared input matrix to adata.obsm[X_mammalianlifespan]\n", "|-----------> ⚠️ Check features in adata finished [0.0043s]\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.0018s]\n", "|-----------> ⚙️ Add predicted ages and clock metadata to adata started\n", "|-----------> ✅ Add predicted ages and clock metadata to adata finished [0.0009s]\n", "|-----> 🕒 Processing clock: mammalianfemale\n", "|-----------> ⚙️ Load clock started\n", "|-----------------> Data found in pyaging_data/mammalianfemale.pt\n", "|-----------> ✅ Load clock finished [0.4532s]\n", "|-----------> ⚙️ Check features in adata started\n", "|-----------------? 73 out of 101 features (72.28%) are missing: ['cg01145947', 'cg02053792', 'cg02407848'], etc.\n", "|-----------------> Filling missing features entirely with 0\n", "|-----------------> Added prepared input matrix to adata.obsm[X_mammalianfemale]\n", "|-----------> ⚠️ Check features in adata finished [0.0135s]\n", "|-----------> ⚙️ Predict ages with model started\n", "|-----------------> There is no preprocessing necessary\n", "|-----------------> The postprocessing method is sigmoid\n", "|-----------------> in progress: 100.0000%\n", "|-----------> ✅ Predict ages with model finished [0.0051s]\n", "|-----------> ⚙️ Add predicted ages and clock metadata to adata started\n", "|-----------> ✅ Add predicted ages and clock metadata to adata finished [0.0024s]\n", "|-----> 🎉 Done! [3.0755s]\n" ] } ], "source": [ "pya.pred.predict_age(adata_mammalian, ['Mammalian1', 'Mammalian2', 'Mammalian3', \"MammalianLifespan\", \"MammalianFemale\"])" ] }, { "cell_type": "markdown", "id": "de8ed8d2-213a-4ca9-a0d0-e141ba7d1789", "metadata": {}, "source": [ "Note that RRBS clocks are in units of months whereas the mammalian clocks are in units of years." ] }, { "cell_type": "code", "execution_count": 16, "id": "32c3c7fe-acbb-4fe8-b438-5c664a533c41", "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", "
mammalian1mammalian2mammalian3mammalianlifespanmammalianfemale
GSM37526312.537895-0.009800-0.0484141.2021340.732238
GSM37526253.3539350.064448-0.0482031.4800800.952105
GSM37526344.4906100.813899-0.0350591.4002780.978554
GSM37526203.6038020.122934-0.0469581.6447010.948646
GSM37526222.9512630.005617-0.0477251.3852820.741336
GSM37526375.7185150.895781-0.0376491.3997690.975114
GSM45582167.4562450.684752-0.0121951.4481190.785830
GSM37526435.8819430.880053-0.0376561.4084830.969979
GSM45582136.7200800.855574-0.0262111.4772360.821308
GSM37526406.4529340.766701-0.0317631.3718310.938162
GSM45582225.0492470.120607-0.0451091.3667290.789795
GSM45582196.0987100.857243-0.0298781.4221440.796568
GSM37526282.7469490.098284-0.0479891.4805990.794259
GSM37526172.7398680.078172-0.0482491.4840910.807197
\n", "
" ], "text/plain": [ " mammalian1 mammalian2 mammalian3 mammalianlifespan \\\n", "GSM3752631 2.537895 -0.009800 -0.048414 1.202134 \n", "GSM3752625 3.353935 0.064448 -0.048203 1.480080 \n", "GSM3752634 4.490610 0.813899 -0.035059 1.400278 \n", "GSM3752620 3.603802 0.122934 -0.046958 1.644701 \n", "GSM3752622 2.951263 0.005617 -0.047725 1.385282 \n", "GSM3752637 5.718515 0.895781 -0.037649 1.399769 \n", "GSM4558216 7.456245 0.684752 -0.012195 1.448119 \n", "GSM3752643 5.881943 0.880053 -0.037656 1.408483 \n", "GSM4558213 6.720080 0.855574 -0.026211 1.477236 \n", "GSM3752640 6.452934 0.766701 -0.031763 1.371831 \n", "GSM4558222 5.049247 0.120607 -0.045109 1.366729 \n", "GSM4558219 6.098710 0.857243 -0.029878 1.422144 \n", "GSM3752628 2.746949 0.098284 -0.047989 1.480599 \n", "GSM3752617 2.739868 0.078172 -0.048249 1.484091 \n", "\n", " mammalianfemale \n", "GSM3752631 0.732238 \n", "GSM3752625 0.952105 \n", "GSM3752634 0.978554 \n", "GSM3752620 0.948646 \n", "GSM3752622 0.741336 \n", "GSM3752637 0.975114 \n", "GSM4558216 0.785830 \n", "GSM3752643 0.969979 \n", "GSM4558213 0.821308 \n", "GSM3752640 0.938162 \n", "GSM4558222 0.789795 \n", "GSM4558219 0.796568 \n", "GSM3752628 0.794259 \n", "GSM3752617 0.807197 " ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adata_mammalian.obs" ] }, { "cell_type": "markdown", "id": "615f8fbf-f5e1-4af9-a2a0-5f4f781001fe", "metadata": {}, "source": [ "## Get citation" ] }, { "cell_type": "markdown", "id": "0fe55edd-9271-4b41-857d-ef3fceafc2a6", "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": 17, "id": "cdf4c609-7a24-4c3f-a891-647315b77d54", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'clock_name': 'thompson',\n", " 'data_type': 'methylation',\n", " 'species': 'Mus musculus',\n", " 'year': 2018,\n", " 'approved_by_author': '✅',\n", " 'citation': 'Thompson, Michael J., et al. \"A multi-tissue full lifespan epigenetic clock for mice.\" Aging (Albany NY) 10.10 (2018): 2832.',\n", " 'doi': 'https://doi.org/10.18632/aging.101590',\n", " 'notes': None,\n", " 'version': None}" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adata.uns['thompson_metadata']" ] }, { "cell_type": "code", "execution_count": 18, "id": "466b1bc5-6207-469e-b479-260bbf55f2a7", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'clock_name': 'meer',\n", " 'data_type': 'methylation',\n", " 'species': 'Mus musculus',\n", " 'year': 2018,\n", " 'approved_by_author': '⌛',\n", " 'citation': 'Meer, Margarita V., et al. \"A whole lifespan mouse multi-tissue DNA methylation clock.\" Elife 7 (2018): e40675.',\n", " 'doi': 'https://doi.org/10.7554/eLife.40675',\n", " 'notes': None,\n", " 'version': None}" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adata.uns['meer_metadata']" ] }, { "cell_type": "code", "execution_count": 19, "id": "2808782d-04d3-4527-8328-b18a583cf15b", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'clock_name': 'petkovich',\n", " 'data_type': 'methylation',\n", " 'species': 'Mus musculus',\n", " 'year': 2017,\n", " 'approved_by_author': '⌛',\n", " 'citation': 'Petkovich, Daniel A., et al. \"Using DNA methylation profiling to evaluate biological age and longevity interventions.\" Cell metabolism 25.4 (2017): 954-960.',\n", " 'doi': 'https://doi.org/10.1016/j.cmet.2017.03.016',\n", " 'notes': None,\n", " 'version': None}" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adata.uns['petkovich_metadata']" ] }, { "cell_type": "code", "execution_count": 20, "id": "ec210b8b-9ba1-45df-9b16-80a5f5ac86f4", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'clock_name': 'stubbs',\n", " 'data_type': 'methylation',\n", " 'species': 'Mus musculus',\n", " 'year': 2017,\n", " 'approved_by_author': '⌛',\n", " 'citation': 'Stubbs, Thomas M., et al. \"Multi-tissue DNA methylation age predictor in mouse.\" Genome biology 18 (2017): 1-14.',\n", " 'doi': 'https://doi.org/10.1186/s13059-017-1203-5',\n", " 'notes': None,\n", " 'version': None}" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adata.uns['stubbs_metadata']" ] }, { "cell_type": "code", "execution_count": 21, "id": "78ed8018-bfdc-42b7-b814-d3a55184fc05", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'clock_name': 'mammalian1',\n", " 'data_type': 'methylation',\n", " 'species': 'multi',\n", " 'year': 2023,\n", " 'approved_by_author': '⌛',\n", " 'citation': 'Lu, A. T., et al. \"Universal DNA methylation age across mammalian tissues.\" Nature aging 3.9 (2023): 1144-1166.',\n", " 'doi': 'https://doi.org/10.1038/s43587-023-00462-6',\n", " 'notes': None,\n", " 'version': None}" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adata_mammalian.uns['mammalian1_metadata']" ] }, { "cell_type": "code", "execution_count": 22, "id": "f4b9c4c7-d23b-4850-abf2-5f9a092c7bc4", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'clock_name': 'mammalianlifespan',\n", " 'data_type': 'methylation',\n", " 'species': 'multi',\n", " 'year': 2023,\n", " 'approved_by_author': '⌛',\n", " 'citation': 'Li, Caesar Z., et al. \"Epigenetic predictors of species maximum lifespan and other life history traits in mammals.\" bioRxiv (2023): 2023-11.',\n", " 'doi': 'https://doi.org/10.1101/2023.11.02.565286',\n", " 'notes': None,\n", " 'version': None}" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adata_mammalian.uns['mammalianlifespan_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 }