Building Genome Reference Files

Note: Building custom reference files and token dictionaries is completely optional. If you are working with standard human or mouse datasets, you can skip this tutorial and directly use the pre-compiled reference files provided at data\gene_dict\.

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.

We will walk through setting up the inputs, initializing the reference maker, and exporting the indices.

import anndata as ad, hdf5plugin
from scdynomics.utils.ref import Maker

1. Initializing the Reference Maker

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.

For the transcription factor (TF) list, we are using reference files retrieved from the SCENIC+ project.

---

# Define paths to the GFF and TF list files
mm10_gff3_path = "../../data/gff/Mus_musculus.GRCm38.102.gff3.gz"
mm10_tf_list_path = "../../data/tf_list/allTFs_mm.txt"

# Instantiate the Maker class
reference_maker = Maker(
    gff_path=mm10_gff3_path,
    tf_list_path=mm10_tf_list_path,
    remove_non_coding=True,
    promoter_tss_upstream=1500,
    promoter_tss_downstream=500
)

print(f"Loaded {len(reference_maker.gene_dict)} genes from the mm10 GFF.")

Loaded 22452 genes from the mm10 GFF.

2. Generating and Saving References

With the data parsed, we can now output our reference formats.

1. JSON Dictionary: Contains the mapped genes, transcripts, proteins, and promoters.

2. Gene Expression AnnData Index: Creates a structural index for expression matrices, sorting transcription factors to the top.

---

output_json_path = "ref.json"
output_h5ad_path = "ref.h5ad"

# 1. Output the dictionary
ref_dict = reference_maker.output_dict(save_path=output_json_path)
print("Saved JSON reference.")

# 2. Output the Gene Expression AnnData index
gex_adata = reference_maker.output_gex_index(
    skip_mt=False, # Keeping mitochondrial genes in the reference
    save_path=output_h5ad_path
)
print("Saved AnnData reference.")

Saved JSON reference.
Saved AnnData reference.

3. Inspecting the Result

Let’s look at the generated AnnData object to verify our gene indexing and TF tagging worked correctly.

---

# View the AnnData structure
print(gex_adata)

# View the variable (gene) metadata
display(gex_adata.var.head())

AnnData object with n_obs × n_vars = 0 × 22452
    var: 'name', 'type', 'chr', 'start', 'end', 'strand', 'promoters', 'mt'
    uns: 'n_tf'

	name	type	chr	start	end	strand	promoters	mt
id
ENSMUSG00000025902	Sox17	tf	1	4490931	4497354	-	chr1:4495913-4497913;chr1:4496257-4498257;chr1...	False
ENSMUSG00000025912	Mybl1	tf	1	9667415	9700209	-	chr1:9699709-9701709;chr1:9699524-9701524;chr1...	False
ENSMUSG00000099032	Tcf24	tf	1	9960163	9967932	-	chr1:9967432-9969432;chr1:9962309-9964309	False
ENSMUSG00000042414	Prdm14	tf	1	13113457	13127163	-	chr1:13126663-13128663	False
ENSMUSG00000005886	Ncoa2	tf	1	13139105	13374083	-	chr1:13371939-13373939;chr1:13373583-13375583;...	False

4. Building Input Token Dictionary

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.

Here is a breakdown of what the code accomplishes:

• Common Tokens: Assigns baseline tokens (<MASK> = 0, <PAD> = 1) used for sequence padding and masking during model training.

• Species Integration: Loads the pre-computed AnnData references for both mouse (mm10) and human (hg38), explicitly filtering out mitochondrial genes.

• Sequential Numbering: Iterates through the genes of each species, assigning a unique, sequential integer ID to each gene.

• 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.

• TF Counts: Extracts and stores the total number of TFs (n_tf) for each species, directly indcating the boundary between TFs and other genes.

---

# Initialization
token_dict = {
    'common': dict(),
    'mm10': dict(),
    'hg38': dict(),
    'reverse': dict(),
    'n_tf':{
        'mm10': None,
        'hg38': None
    }
}

# Read in references
mm10_ref = ad.read_h5ad(
    '../../data/gene_dict/mm10.pan_promoter.ind.h5ad',
    backed='r'
)
hg38_ref = ad.read_h5ad(
    '../../data/gene_dict/hg38.pan_promoter.ind.h5ad',
    backed='r'
)

# Remove MT genes
mm10_ref = mm10_ref[:, mm10_ref.var['mt'] != True]
hg38_ref = hg38_ref[:, hg38_ref.var['mt'] != True]

# common tokens no matter of species
token_dict['common']['<MASK>'] = 0
token_dict['common']['<PAD>'] = 1
used_token = len(token_dict['common'])

# mm10 tokens
token_dict['mm10'] = {
    gene: i+used_token
    for i, gene in enumerate(mm10_ref.var.index)
}
used_token += len(token_dict['mm10'])

# hg38 tokens
token_dict['hg38'] = {
    gene: i+used_token
    for i, gene in enumerate(hg38_ref.var.index)
}
used_token += len(token_dict['hg38'])

# Reverse dictionary
for k,v in token_dict['common'].items():
    assert v not in token_dict['reverse']
    token_dict['reverse'] [v] = {
        'gene': k,
        'type': 'common'
    }
for k,v in token_dict['mm10'].items():
    assert v not in token_dict['reverse']
    token_dict['reverse'] [v] = {
        'gene': k,
        'type': 'mm10'
    }
for k,v in token_dict['hg38'].items():
    assert v not in token_dict['reverse']
    token_dict['reverse'] [v] = {
        'gene': k,
        'type': 'hg38'
    }

# n_tf
token_dict['n_tf']['mm10'] = int(mm10_ref.uns['n_tf'])
token_dict['n_tf']['hg38'] = int(hg38_ref.uns['n_tf'])