Preprocessing Data

This notebook outlines the data preprocessing steps required for scDynOmics. Before starting, ensure your data meets the following formatting requirements:

Gene Identifiers: The adata.var index must use Ensembl IDs.
ATAC Annotations: Peak annotations must follow the chrN:start-end format.

[1]:

import os
import anndata as ad, hdf5plugin
import scdynomics.utils.pp as pp
from scdynomics.utils.data.reader import scMultiomics_10x_h5

temp_data_name = "temp.h5"

# Download the sample data from 10x Genomics and save it as "temp.h5"
os.system("curl -O https://cf.10xgenomics.com/samples/cell-arc/2.0.0/e18_mouse_brain_fresh_5k/e18_mouse_brain_fresh_5k_filtered_feature_bc_matrix.h5")
os.rename("e18_mouse_brain_fresh_5k_filtered_feature_bc_matrix.h5", temp_data_name)

# Process the data and load it into AnnData objects
raw_gex, raw_peaks = scMultiomics_10x_h5(temp_data_name)
os.remove(temp_data_name)
print("\nData loaded successfully.")
print(f"Raw GEX adata: {raw_gex}")
print(f"Raw PEAKS adata: {raw_peaks}")

# Set the reference data paths
ref_adata = ad.read_h5ad("../../data/gene_dict/mm10.pan_promoter.ind.h5ad")
ref_dict_path = "../../data/gene_dict/mm10.pan_promoter.dict.json.gz"

  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
                                 Dload  Upload   Total   Spent    Left  Speed
100 85.0M  100 85.0M    0     0   106M      0 --:--:-- --:--:-- --:--:--  106M


Data loaded successfully.
Raw GEX adata: View of AnnData object with n_obs × n_vars = 4878 × 32284
Raw PEAKS adata: View of AnnData object with n_obs × n_vars = 4878 × 172193

1. Pretraining Pipeline

The pretraining phase requires aligning your data to a standard reference.

Expected Input: Multimodal samples containing two modalities that represent distinct time points.

1.1 Mapping Raw Inputs to Reference Feature Space

We first align the raw gene expression and peak data to the provided reference genome to ensure consistent feature spaces.

[2]:

# Align the raw data to the reference data
gex_adata = pp.align.align_genes_to_reference(
    raw_gex,
    ref_adata = ref_adata,
    handle_duplicated = 'sum',
)

# Map the peaks to the reference promoters
atac_adata = pp.align.align_peaks_to_reference(
    raw_peaks,
    ref_adata = ref_adata,
    ref_dict_path = ref_dict_path
)

# Check
assert gex_adata.obs.index.equals(atac_adata.obs.index)
assert gex_adata.var.index.equals(atac_adata.var.index)

# Combining into same AnnData object
pt_adata = ad.AnnData(
    X = gex_adata.X,
    obs = gex_adata.obs,
    var = gex_adata.var,
    layers = {
        "pan_promoter": atac_adata.X
    }
)

print(f"Overview of the processed AnnData: {pt_adata}")

/home/ag_klughammer/gyu/.conda/envs/pg_0/lib/python3.11/site-packages/scdynomics/utils/pp/align.py:108: UserWarning: Some chromosomes were not found in the reference and were ignored: {'GL456233.1': 1, 'JH584304.1': 1, 'GL456216.1': 1, 'JH584292.1': 1, 'JH584295.1': 1}
  warnings.warn(

Overview of the processed AnnData: AnnData object with n_obs × n_vars = 4878 × 22452
    var: 'name', 'type', 'chr', 'start', 'end', 'strand', 'promoters', 'mt'
    layers: 'pan_promoter'

1.2 Quality Control (QC) and Normalization

Next, we filter out low-quality cells, remove mitochondrial genes, and normalize the data.

[ ]:

# WARNING: adata.X seems to be already log-transformed. -> it is not, so this can be ignored.
# https://github.com/scverse/scanpy/issues/1333
pped_adata = pp.pipeline(
    pt_adata,
    ref_adata = ref_adata,
    layers = ['x', 'pan_promoter'],
    matching_strategy = 'align',
    filtering = True,
    min_genes = 100, # Only filtering cells
    remove_outliers = True,
    log1p = True,
    normalizing = True,
    target_sum = 1e4,
    methods = ['tpm', 'cpm_promoter'],
    remove_mt = True, # Remove mitochondrial genes
)

print(f"\nOverview of the adata object: {pped_adata}")

Genes are already aligned to the reference.
WARNING: adata.X seems to be already log-transformed.

Overview of the adata object: AnnData object with n_obs × n_vars = 4850 × 22439
    var: 'name', 'type', 'chr', 'start', 'end', 'strand', 'promoters'
    uns: 'log1p'
    layers: 'pan_promoter'

2. Fine-Tuning Pipeline

Preparing data for downstream fine-tuning requires a slightly different, more flexible approach compared to strict pretraining:

Feature Selection (Optional): To reduce dimensionality, you may optionally subset the dataset to only include Highly Variable Genes (HVGs) or peaks. If you prefer to retain the full feature space, this step can be bypassed.
Filtering: Strict cell and gene filtering is highly dependent on your specific downstream task. However, to maintain computational stability, establishing a baseline filter of min_genes = 1 and min_cells = 1 is strongly recommended.
Flexibility: The fine-tuning preprocessing pipeline is highly modular. It allows you to customize or skip steps based on the state of your input data. For example, if your dataset has already been normalized using a custom external method, you could skip the normalization step.
Hyperparameters: Please ensure you adjust the functional pipeline arguments below to best fit the distribution and specific needs of your dataset.

2.1: Multimodal Processing

This example uses a 10-cell sample from the Gastrulation dataset, formatted as processed RNA velocity with scVelo.

Note: For standard scMultiomics data, you must map the feature space first, as demonstrated in the Pretraining section 1.1 .

[4]:

# Load the raw data and set the index as the Ensembl gene ID (Accession)
ft_adata = ad.read_h5ad("../../data/sample/raw.gastrulation_e75.h5ad") # example data included in the scDynOmics github repository relative to this notebook
ft_adata.var['name'] = ft_adata.var.index
ft_adata.var.index = ft_adata.var['Accession']
ft_adata.var = ft_adata.var.drop(columns=['Accession'])

print(f"Overview of the adata object: {ft_adata} \n")

Overview of the adata object: AnnData object with n_obs × n_vars = 10 × 53801
    obs: 'barcode', 'sample', 'stage', 'sequencing.batch', 'theiler', 'doub.density', 'doublet', 'cluster', 'cluster.sub', 'cluster.stage', 'cluster.theiler', 'stripped', 'celltype', 'colour', 'umapX', 'umapY', 'haem_gephiX', 'haem_gephiY', 'haem_subclust', 'endo_gephiX', 'endo_gephiY', 'endo_trajectoryName', 'endo_trajectoryDPT', 'endo_gutX', 'endo_gutY', 'endo_gutDPT', 'endo_gutCluster', 'cell_velocyto_loom'
    var: 'Chromosome', 'End', 'Start', 'Strand', 'name'
    obsm: 'X_pca', 'X_umap'
    layers: 'spliced', 'unspliced'

[5]:

# WARNING: adata.X seems to be already log-transformed. -> it is not, so this can be ignored.
# https://github.com/scverse/scanpy/issues/1333
pped_adata = pp.pipeline(
    ft_adata,
    ref_adata = ref_adata,
    layers = ['spliced', 'unspliced'],
    matching_strategy = 'stratify',
    filtering = True,
    min_genes = 1,
    min_cells = 1,
    remove_outliers = True,
    log1p = True,
    normalizing = True,
    target_sum = 1e4,
    methods = ['tpm', 'tpm'],
    remove_mt = True,
)

print(f"\nOverview of the adata object: {pped_adata} \n")

WARNING: adata.X seems to be already log-transformed.

Overview of the adata object: AnnData object with n_obs × n_vars = 9 × 7870
    obs: 'barcode', 'sample', 'stage', 'sequencing.batch', 'theiler', 'doub.density', 'doublet', 'cluster', 'cluster.sub', 'cluster.stage', 'cluster.theiler', 'stripped', 'celltype', 'colour', 'umapX', 'umapY', 'haem_gephiX', 'haem_gephiY', 'haem_subclust', 'endo_gephiX', 'endo_gephiY', 'endo_trajectoryName', 'endo_trajectoryDPT', 'endo_gutX', 'endo_gutY', 'endo_gutDPT', 'endo_gutCluster', 'cell_velocyto_loom'
    var: 'Chromosome', 'End', 'Start', 'Strand', 'name.original', 'name', 'type', 'chr', 'start', 'end', 'strand', 'promoters'
    uns: 'log1p'
    obsm: 'X_pca', 'X_umap'
    layers: 'spliced', 'unspliced'

2.2: Monomodal Processing

If you are not using multimodal data, you can process raw monomodal read counts directly. This approach also supports Highly Variable Gene (HVG) selection.

[6]:

# Default to use adata.X. No need to specify if you don't have any other layers.
pped_adata = pp.pipeline(
    ft_adata,
    ref_adata = ref_adata,
    matching_strategy = 'stratify',
    filtering = True,
    min_genes = 1,
    min_cells = 1,
    remove_outliers = True,
    log1p = True,
    normalizing = True,
    target_sum = 1e4,
    methods = ['tpm'],
    remove_mt = True,
    n_top_genes = 2000,
    hvg_flavor = 'seurat',
)

print(f"Overview of the adata object: {pped_adata} \n")

/home/ag_klughammer/gyu/.conda/envs/pg_0/lib/python3.11/site-packages/scanpy/preprocessing/_simple.py:293: ImplicitModificationWarning: Trying to modify attribute `.var` of view, initializing view as actual.
  adata.var["n_cells"] = number

Overview of the adata object: AnnData object with n_obs × n_vars = 8 × 2000
    obs: 'barcode', 'sample', 'stage', 'sequencing.batch', 'theiler', 'doub.density', 'doublet', 'cluster', 'cluster.sub', 'cluster.stage', 'cluster.theiler', 'stripped', 'celltype', 'colour', 'umapX', 'umapY', 'haem_gephiX', 'haem_gephiY', 'haem_subclust', 'endo_gephiX', 'endo_gephiY', 'endo_trajectoryName', 'endo_trajectoryDPT', 'endo_gutX', 'endo_gutY', 'endo_gutDPT', 'endo_gutCluster', 'cell_velocyto_loom', 'n_genes', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_20_genes', 'outlier'
    var: 'Chromosome', 'End', 'Start', 'Strand', 'name.original', 'name', 'type', 'chr', 'start', 'end', 'strand', 'promoters', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts', 'n_cells', 'highly_variable', 'means', 'dispersions', 'dispersions_norm'
    uns: 'log1p', 'hvg'
    obsm: 'X_pca', 'X_umap'
    layers: 'spliced', 'unspliced'