Single-Cell RNA-seq Walkthrough: How I Mapped Cell-Type-Specific Regulatory Networks Across Cardiac Development

A step-by-step breakdown of how I analysed three scRNA-seq datasets: fetal, infant, and adult cardiac outflow tract, to uncover transcription factor activity and developmental gene regulation.

Why This Project?

The heart doesn’t stop changing after birth. The outflow tract (OFT), the region connecting the heart to the aorta and pulmonary artery, undergoes dramatic structural and cellular remodelling from fetal life through adulthood. Understanding how gene regulation shifts across those developmental stages has real implications for congenital heart disease research and regenerative medicine.

This project gave me three scRNA-seq datasets representing three stages of human cardiac OFT development: fetal, infant, and adult. My job was to go from raw gene expression matrices to biologically meaningful insights, identifying which transcription factors (TFs) are driving cell identity at each stage, and how those regulatory networks differ across development.

Here’s exactly how I did it.

The Dataset

Sample Stage Cells: Sample 1 — Fetal ~17,457 cells Sample 2 — Embryonic ~7,287 cells Sample 3 — Adult ~5,430 cells

Total: ~30,166 cells · 14 clusters · 9 major cell types

Each sample came as a gene expression matrix (genes × cells). The goal: Compare TF activity and cell-type-specific gene regulation across all three stages.

UMAP of our dataset samples

The Workflow at a Glance

A note on this walkthrough: My original analysis was conducted using computational tools: TF2DNA, CELLxGENE, and GOrilla, through their respective interfaces.

(Me) Working on my project

The code snippets throughout this post are provided as a programmatic equivalent for readers who want to replicate or extend this workflow in Python/R. Think of them as your coding companion to the analysis! 🧬

Raw Gene Expression Matrices
        ↓
Quality Control & Preprocessing
        ↓
Dimensionality Reduction (PCA → UMAP)
        ↓
Clustering & Cell Type Annotation
        ↓
TF Activity Scoring (TF2DNA)
        ↓
Differential Gene Expression (CELLxGENE)
        ↓
Pathway Enrichment (GOrilla)
        ↓
Visualisation & Interpretation

Step 1: Quality Control & Preprocessing

Before any analysis, I filtered out low-quality cells, those with too few detected genes (likely empty droplets) or too many (likely doublets), and cells with high mitochondrial gene content (a sign of dying cells).

import scanpy as sc

# Load each sample
adata_fetal = sc.read_h5ad("fetal_OFT.h5ad")
adata_infant = sc.read_h5ad("infant_OFT.h5ad")
adata_adult  = sc.read_h5ad("adult_OFT.h5ad")
# QC metrics
sc.pp.calculate_qc_metrics(adata_fetal, inplace=True)
# Filter
sc.pp.filter_cells(adata_fetal, min_genes=200)
sc.pp.filter_genes(adata_fetal, min_cells=3)
# Remove high mitochondrial content cells
adata_fetal = adata_fetal[adata_fetal.obs['pct_counts_mt'] < 20]

I then normalised each dataset to 10,000 counts per cell and log-transformed to stabilise variance.

sc.pp.normalize_total(adata_fetal, target_sum=1e4)
sc.pp.log1p(adata_fetal)

Step 2: Dimensionality Reduction & Clustering

With ~30,000 cells across three samples, direct visualisation isn’t possible. I used PCA to reduce to 50 components, then UMAP for 2D visualisation, a standard and interpretable approach for scRNA-seq.

sc.pp.highly_variable_genes(adata_fetal, n_top_genes=2000)
sc.tl.pca(adata_fetal, svd_solver='arpack')
sc.pp.neighbors(adata_fetal, n_neighbors=10, n_pcs=40)
sc.tl.umap(adata_fetal)
sc.tl.leiden(adata_fetal, resolution=0.5)  # Community-based clustering

After merging the three datasets and running batch correction (to remove technical variation between samples), I identified 14 distinct clusters corresponding to 9 major cell types, including cardiomyocytes, smooth muscle cells, fibroblasts, endothelial cells, and neural crest-derived cells.

💡 Tip: Always check your UMAP coloured by sample of origin before and after batch correction. If samples are completely separated before correction, your biology will be confounded by batch effects.

Step 3: TF Activity Scoring with TF2DNA

This was the most novel part of the analysis. Rather than just looking at which genes are expressed, I wanted to know which transcription factors are actively regulating those genes, a much more mechanistically informative question.

I used the TF2DNA database to retrieve TF binding site information for 195 TFs. For each TF, a regulon (the set of genes it controls) was defined, and activity scores were calculated per cell using a scoring approach similar to AUCell.

# Pseudocode — adapted for TF2DNA regulons
import decoupler as dc

# Load TF-target network from TF2DNA
network = pd.read_csv("TF2DNA_network.csv")  # TF, target, weight

# Run activity inference
dc.run_aucell(
    mat=adata,
    net=network,
    source='TF',
    target='gene',
    use_raw=True
)

Why TF activity scoring instead of just differential expression? Gene expression tells you what is happening. TF activity tells you who is in charge. Two cells can express similar genes but be driven by completely different master regulators, and that distinction matters enormously in developmental biology.

TF expression visualisation using TF2DNA

Step 4: Visualisation & Exploration with CELLxGENE

For interactive exploration of gene expression and cluster annotation, I used the CELLxGENE app. This browser-based tool lets you query any gene across any cluster in real time without writing code.

This was particularly useful for:

• Validating cluster annotations against known marker genes
• Exploring TF expression patterns across developmental stages
• Comparing fetal vs adult expression of specific regulators visually

The dot plot was especially informative. It showed both the fraction of cells expressing a gene and the mean expression level, giving a cleaner picture than a standard heatmap.

Dot Plot of the genes of interest

Step 5: Differential Gene Expression Across Developmental Stages

With clusters annotated, I ran pairwise differential expression between developmental stages within each cell type using a Wilcoxon rank-sum test.

sc.tl.rank_genes_groups(
    adata,
    groupby='stage',        # fetal / infant / adult
    method='wilcoxon',
    key_added='stage_DEGs'
)

sc.pl.rank_genes_groups_dotplot(
    adata,
    key='stage_DEGs',
    n_genes=5,
    groupby='stage'
)

Key findings from this step:

• Several TFs (including NFKB1, SP1, SP4, FP1) were highly expressed in the fetal sample but showed markedly reduced activity in the adult
• GTFIA, ZNF226, USF2, and ZNF189 were present in the fetal but absent or low in adult tissue
• MEF2C showed major behavioural differences. Highly expressed in fetal cardiomyocytes, absent in adults
• SPAS showed consistent expression across NR (HC, SP1, ZNF226, ZNF238, ZNP189), suggesting a conserved regulatory role

Comparative Analysis of the datasets on CELLxGENE app

Step 6: Pathway Enrichment with GOrilla

To understand the biological meaning of the differentially expressed genes, I ran pathway enrichment analyses using GOrilla (Gene Ontology enrichment analysis and visualisation tool).

I submitted ranked gene lists for each developmental stage and each major cell type, asking: What biological processes are enriched in genes that are upregulated in fetal vs adult cardiac OFT cells?

Key enriched pathways included:

• Heart morphogenesis and OFT development: expected, and reassuring as a validation
• Neural crest cell migration: consistent with the known role of neural crest cells in OFT formation
• Transcription factor binding and DNA-templated regulation: confirming TF activity is the dominant regulatory mechanism in this tissue
• Hematopoietic lineage specification: interesting, suggesting shared regulatory logic between cardiac and blood progenitors

Pathway Enrichment (DAG) of the TFs representing the processes they are involved in

Step 7: Publication-Standard Visualisation

The final outputs were produced to Expression Atlas-style publication standards, including:

• UMAP plots coloured by cell type, developmental stage, and TF activity score
• Dot plots showing TF expression across clusters and stages
• Heatmaps of top DEGs per cell type
• Pathway enrichment bubble plots from GOrilla outputs

# Example: UMAP coloured by TF activity
sc.pl.umap(
    adata,
    color=['MEF2C_activity', 'NFKB1_activity', 'stage'],
    cmap='RdYlBu_r',
    size=5,
    save='TF_activity_UMAP.pdf'
)

💡 Always export figures as vector formats (PDF/SVG) for publication. Raster images (PNG/JPG) lose quality when scaled.

Key Takeaways

1. TF activity scoring > raw expression for developmental questions. Knowing a gene is expressed tells you little. Knowing it’s actively driving a regulatory programme tells you a lot.

2. Comparative scRNA-seq requires careful batch correction. Fetal, infant, and adult tissues are processed differently and at different times. Without batch correction, you’re comparing technical noise, not biology.

3. CELLxGENE is underused in academic workflows. It’s not just for browsing public datasets; it’s a powerful exploratory tool for your own data.

4. GOrilla is simple but effective. It doesn’t require a custom R environment or complex input formats. For quick, clean GO enrichment on a ranked list, it’s hard to beat.

Tools Used

Scanpy (Python), TF2DNA (TF-target network database), CELLxGENE, GOrilla (Gene Ontology pathway enrichment), UMAP, Leiden algorithm, Community detection/clustering

What’s Next?

This analysis opens up several follow-up questions I’d love to explore:

• Can we use pseudotime analysis (Monocle/scVelo) to trace developmental trajectories through the OFT cell types?
• Do the TF regulatory networks identified here overlap with known congenital heart disease risk genes?
• Could multi-omics integration (scATAC-seq + scRNA-seq) validate the TF binding activity we inferred computationally?

If you’ve worked on similar cardiac or developmental scRNA-seq datasets, I’d love to compare notes. Drop a comment or find me on LinkedIn or GitHub.

This analysis was conducted as part of my MSc in Bioinformatics and Systems Biology at the University of Manchester, supervised by Prof. Nicoletta Bobola.

Poster presentation of my thesis

Tags: Single-Cell RNA-seq · Bioinformatics · Genomics · Computational Biology · Transcription Factors · Heart Development · Python · Scanpy