{ "cells": [ { "cell_type": "markdown", "id": "0cb8e51a", "metadata": {}, "source": [ "# Building Genome Reference Files\n", "*Note: \n", "Building custom reference files and token dictionaries is completely optional.\n", "If you are working with standard human or mouse datasets,\n", " you can skip this tutorial and directly use the pre-compiled reference files provided at* `data\\gene_dict\\`.\n", "\n", "This tutorial demonstrates how to parse GFF annotation files, annotate transcription factors (TFs), and generate the required JSON and AnnData (`.h5ad`) reference files for scDynOmics.\n", "\n", "We will walk through setting up the inputs, initializing the reference maker, and exporting the indices." ] }, { "cell_type": "code", "execution_count": 7, "id": "f4387c26", "metadata": {}, "outputs": [], "source": [ "import anndata as ad, hdf5plugin\n", "from scdynomics.utils.ref import Maker" ] }, { "cell_type": "markdown", "id": "75fc56da", "metadata": {}, "source": [ "## 1. Initializing the Reference Maker\n", "\n", "To ensure this tutorial is easy to follow and fully reproducible, let's use the default `GFF3` file already provided in the `data` directory. The included `GFF3` annotations are sourced directly from [Ensembl](https://www.ensembl.org/Mus_musculus/Info/Index).\n", "\n", "For the transcription factor (TF) list, we are using reference files retrieved from the [SCENIC+](https://github.com/aertslab/scenicplus/tree/main/resources) project.\n", "\n", "---" ] }, { "cell_type": "code", "execution_count": 2, "id": "790e2650", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Loaded 22452 genes from the mm10 GFF.\n" ] } ], "source": [ "# Define paths to the GFF and TF list files\n", "mm10_gff3_path = \"../../data/gff/Mus_musculus.GRCm38.102.gff3.gz\"\n", "mm10_tf_list_path = \"../../data/tf_list/allTFs_mm.txt\"\n", "\n", "# Instantiate the Maker class\n", "reference_maker = Maker(\n", " gff_path=mm10_gff3_path,\n", " tf_list_path=mm10_tf_list_path,\n", " remove_non_coding=True,\n", " promoter_tss_upstream=1500,\n", " promoter_tss_downstream=500\n", ")\n", "\n", "print(f\"Loaded {len(reference_maker.gene_dict)} genes from the mm10 GFF.\")" ] }, { "cell_type": "markdown", "id": "73f1a5d2", "metadata": {}, "source": [ "## 2. Generating and Saving References\n", "\n", "With the data parsed, we can now output our reference formats. \n", "\n", "1. **JSON Dictionary:** Contains the mapped genes, transcripts, proteins, and promoters.\n", "2. **Gene Expression AnnData Index:** Creates a structural index for expression matrices, sorting transcription factors to the top.\n", "\n", "---" ] }, { "cell_type": "code", "execution_count": 3, "id": "c31eba11", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Saved JSON reference.\n", "Saved AnnData reference.\n" ] } ], "source": [ "output_json_path = \"ref.json\"\n", "output_h5ad_path = \"ref.h5ad\"\n", "\n", "# 1. Output the dictionary\n", "ref_dict = reference_maker.output_dict(save_path=output_json_path)\n", "print(\"Saved JSON reference.\")\n", "\n", "# 2. Output the Gene Expression AnnData index\n", "gex_adata = reference_maker.output_gex_index(\n", " skip_mt=False, # Keeping mitochondrial genes in the reference\n", " save_path=output_h5ad_path\n", ")\n", "print(\"Saved AnnData reference.\")" ] }, { "cell_type": "markdown", "id": "91d3e25c", "metadata": {}, "source": [ "## 3. Inspecting the Result\n", "\n", "Let's look at the generated AnnData object to verify our gene indexing and TF tagging worked correctly.\n", "\n", "---" ] }, { "cell_type": "code", "execution_count": 4, "id": "e7cc5502", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "AnnData object with n_obs × n_vars = 0 × 22452\n", " var: 'name', 'type', 'chr', 'start', 'end', 'strand', 'promoters', 'mt'\n", " uns: 'n_tf'\n" ] }, { "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", "
nametypechrstartendstrandpromotersmt
id
ENSMUSG00000025902Sox17tf144909314497354-chr1:4495913-4497913;chr1:4496257-4498257;chr1...False
ENSMUSG00000025912Mybl1tf196674159700209-chr1:9699709-9701709;chr1:9699524-9701524;chr1...False
ENSMUSG00000099032Tcf24tf199601639967932-chr1:9967432-9969432;chr1:9962309-9964309False
ENSMUSG00000042414Prdm14tf11311345713127163-chr1:13126663-13128663False
ENSMUSG00000005886Ncoa2tf11313910513374083-chr1:13371939-13373939;chr1:13373583-13375583;...False
\n", "
" ], "text/plain": [ " name type chr start end strand \\\n", "id \n", "ENSMUSG00000025902 Sox17 tf 1 4490931 4497354 - \n", "ENSMUSG00000025912 Mybl1 tf 1 9667415 9700209 - \n", "ENSMUSG00000099032 Tcf24 tf 1 9960163 9967932 - \n", "ENSMUSG00000042414 Prdm14 tf 1 13113457 13127163 - \n", "ENSMUSG00000005886 Ncoa2 tf 1 13139105 13374083 - \n", "\n", " promoters mt \n", "id \n", "ENSMUSG00000025902 chr1:4495913-4497913;chr1:4496257-4498257;chr1... False \n", "ENSMUSG00000025912 chr1:9699709-9701709;chr1:9699524-9701524;chr1... False \n", "ENSMUSG00000099032 chr1:9967432-9969432;chr1:9962309-9964309 False \n", "ENSMUSG00000042414 chr1:13126663-13128663 False \n", "ENSMUSG00000005886 chr1:13371939-13373939;chr1:13373583-13375583;... False " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# View the AnnData structure\n", "print(gex_adata)\n", "\n", "# View the variable (gene) metadata\n", "display(gex_adata.var.head())" ] }, { "cell_type": "markdown", "id": "0d39eca4", "metadata": {}, "source": [ "## 4. Building Input Token Dictionary\n", "\n", "In this final step, we construct a unified token dictionary that maps genes from different species (like mouse `mm10` and human `hg38`) to unique integer IDs. This is crucial for downstream machine learning tasks where text-based gene names must be converted into numerical tokens. \n", "\n", "Here is a breakdown of what the code accomplishes:\n", "* **Common Tokens:** Assigns baseline tokens (`` = 0, `` = 1) used for sequence padding and masking during model training.\n", "* **Species Integration:** Loads the pre-computed AnnData references for both mouse (`mm10`) and human (`hg38`), explicitly filtering out mitochondrial genes.\n", "* **Sequential Numbering:** Iterates through the genes of each species, assigning a unique, sequential integer ID to each gene.\n", "* **Reverse Mapping:** Generates a `reverse` lookup dictionary that allows us to easily map any given integer ID back to its original gene name and source species.\n", "* **TF Counts:** Extracts and stores the total number of TFs (`n_tf`) for each species, directly indcating the boundary between TFs and other genes.\n", "\n", "---" ] }, { "cell_type": "code", "execution_count": null, "id": "5f758bee", "metadata": {}, "outputs": [], "source": [ "# Initialization\n", "token_dict = {\n", " 'common': dict(),\n", " 'mm10': dict(),\n", " 'hg38': dict(),\n", " 'reverse': dict(),\n", " 'n_tf':{\n", " 'mm10': None,\n", " 'hg38': None\n", " }\n", "}\n", "\n", "# Read in references\n", "mm10_ref = ad.read_h5ad(\n", " '../../data/gene_dict/mm10.pan_promoter.ind.h5ad',\n", " backed='r'\n", ")\n", "hg38_ref = ad.read_h5ad(\n", " '../../data/gene_dict/hg38.pan_promoter.ind.h5ad',\n", " backed='r'\n", ")\n", "\n", "# Remove MT genes\n", "mm10_ref = mm10_ref[:, mm10_ref.var['mt'] != True]\n", "hg38_ref = hg38_ref[:, hg38_ref.var['mt'] != True]\n", "\n", "# common tokens no matter of species\n", "token_dict['common'][''] = 0\n", "token_dict['common'][''] = 1\n", "used_token = len(token_dict['common'])\n", "\n", "# mm10 tokens\n", "token_dict['mm10'] = {\n", " gene: i+used_token\n", " for i, gene in enumerate(mm10_ref.var.index)\n", "}\n", "used_token += len(token_dict['mm10'])\n", "\n", "# hg38 tokens\n", "token_dict['hg38'] = {\n", " gene: i+used_token\n", " for i, gene in enumerate(hg38_ref.var.index)\n", "}\n", "used_token += len(token_dict['hg38'])\n", "\n", "# Reverse dictionary\n", "for k,v in token_dict['common'].items():\n", " assert v not in token_dict['reverse'] \n", " token_dict['reverse'] [v] = {\n", " 'gene': k,\n", " 'type': 'common'\n", " }\n", "for k,v in token_dict['mm10'].items():\n", " assert v not in token_dict['reverse'] \n", " token_dict['reverse'] [v] = {\n", " 'gene': k,\n", " 'type': 'mm10'\n", " }\n", "for k,v in token_dict['hg38'].items():\n", " assert v not in token_dict['reverse'] \n", " token_dict['reverse'] [v] = {\n", " 'gene': k,\n", " 'type': 'hg38'\n", " }\n", "\n", "# n_tf\n", "token_dict['n_tf']['mm10'] = int(mm10_ref.uns['n_tf'])\n", "token_dict['n_tf']['hg38'] = int(hg38_ref.uns['n_tf'])" ] } ], "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.11.6" } }, "nbformat": 4, "nbformat_minor": 5 }