How I Used Python to Analyse 40,000 Human Gut Cells and Uncover What Makes Crohn's Disease Different from Colitis

A step-by-step walkthrough of my multi-sample single-cell RNA sequencing project — written for anyone curious about how computational biology actually works

Before we start - what is this article actually about?

Imagine being able to take a tiny sample of tissue from a patient's gut, and instead of just knowing "there are cells here", you could read the activity of every single gene inside every single individual cell - thousands of cells at once.

That's what single-cell RNA sequencing (scRNA-seq) does. And in this article, I'll walk you through how I used Python to analyse data from 40,000 human gut cells across 18 patients to understand what makes two similar gut diseases — Crohn's disease and ulcerative colitis - biologically different from each other.

No biology PhD required. I'll explain every term as we go.

The biological problem - two diseases that look the same but aren't

Crohn's disease (CD) and ulcerative colitis (UC) are both types of Inflammatory Bowel Disease (IBD). Both cause chronic inflammation in the gut, both cause pain and discomfort, and both are lifelong conditions. Doctors have known for decades that they're different diseases, but at a surface level they're easy to confuse.

The key difference is where and how they inflame:

• Crohn's disease can affect any part of the digestive tract, tends to cause deep, patchy inflammation, and is driven largely by a type of immune cell called a myeloid cell (think macrophages - the "pac-man" cells of the immune system)
• Ulcerative colitis affects only the colon and rectum, causes continuous surface inflammation, and is driven more by antibody-producing cells called plasma cells

Understanding these differences at the level of individual cells - not just tissue - is critical for developing better, more targeted treatments. That's where single-cell sequencing comes in.

The data - 18 patients, 40,000 cells, one big challenge

I used a publicly available dataset from the CZ CELLxGENE platform containing colonic mucosa (colon lining) biopsies from 18 donors:

• 6 patients with Crohn's disease
• 6 patients with ulcerative colitis
• 6 healthy controls

Each donor contributed thousands of cells, giving us 46,700 cells total (40,084 after quality filtering).

But here's the problem: when you combine data from 18 different people, collected at different times, processed slightly differently in the lab - the data gets messy. The technical differences between samples (called batch effects) can be so strong that they hide the real biological differences you actually care about.

Think of it like trying to compare photos taken by 18 different cameras with different colour settings. The same object might look completely different just because of the camera, not because the object actually changed.

This is where Harmony comes in.

What is Harmony and why does it matter?

Harmony is a batch correction algorithm. In plain English: it's a mathematical tool that looks at all 18 samples, figures out which differences between samples are due to technical variation (the "camera settings"), and removes those differences — leaving only the real biological signal.

Here's what the data looks like before Harmony:

PCA plot coloured by donor - you can see individual donors clustering separately, meaning donor identity is driving the structure more than the actual biology.

And after Harmony:

UMAP plot coloured by donor - all 18 donors are now mixed uniformly within each cell type cluster. The batch effects are gone.

But crucially - when we colour the same post-Harmony UMAP by disease group instead of donor:

UMAP coloured by disease - Crohn's, UC, and normal cells now separate into distinct regions. The biology is preserved, the noise is removed.

This before/after comparison is the core technical contribution of this project. Without Harmony, any findings could be explained by "oh, that's just because donor 3 was processed differently." With Harmony, we can be confident the differences are real.

Setting up the pipeline in Python

The full analysis uses Scanpy - the standard Python library for single-cell analysis - along with harmonypy for the integration step.

import scanpy as sc
import harmonypy

# Load the dataset
adata = sc.read_h5ad('data/ibd_dataset.h5ad')
# 46,700 cells × 32,354 genes

The data is stored in an AnnData object - think of it as a very specialised spreadsheet where rows are cells, columns are genes, and there's extra space to store metadata like disease status, donor ID, and cell type labels.

Step 1 - Quality control

Not every cell in a sequencing experiment is a real, healthy cell. Some are damaged, some are empty droplets that got accidentally captured, and some are doublets (two cells mistakenly counted as one). We filter these out using three metrics:

# Remove low quality cells
adata = adata[adata.obs['n_genes_by_counts'] > 200]   # too few genes = empty droplet
adata = adata[adata.obs['n_genes_by_counts'] < 6000]  # too many genes = doublet
adata = adata[adata.obs['pct_counts_mt'] < 30]        # high mitochondrial % = dying cell

Why mitochondrial genes? When a cell is dying or stressed, the nucleus breaks down and releases its RNA - but mitochondria (the cell's energy factories) have their own separate DNA and RNA that stays intact longer. So a high percentage of mitochondrial gene reads is a reliable sign of a low-quality cell.

After filtering: 40,084 cells remain.

Step 2 - Normalisation

Different cells capture different amounts of RNA simply due to technical variation in the sequencing process. A cell with 10,000 RNA molecules captured will look like it expresses every gene more than a cell with only 2,000 captured — even if their true biology is identical.

We normalise by scaling every cell to have the same total count (10,000), then apply a log transformation to reduce the influence of very highly expressed genes:

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

Step 3 - Finding highly variable genes

Out of 32,354 genes, most are either not expressed at all or expressed at the same level in every cell (housekeeping genes that keep basic cell functions running). These genes add noise without adding information.

We select only the highly variable genes - genes whose expression varies meaningfully between cells - for downstream analysis. We found 2,873 of these, using batch_key='donor_id' to ensure we pick genes that are variable across all donors, not just one:

sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3,
                             min_disp=0.5, batch_key='donor_id')
# Result: 2,873 highly variable genes

Step 4 - PCA and Harmony integration

PCA (Principal Component Analysis) reduces our 2,873-gene matrix into 50 dimensions that capture the most important variation. Then Harmony corrects for batch effects within this PCA space:

sc.tl.pca(adata, svd_solver='arpack')
sc.external.pp.harmony_integrate(adata, key='donor_id', max_iter_harmony=20)
# Harmony converged in just 6 iterations

The fact that Harmony converged in only 6 out of 20 possible iterations is a good sign - it means the batch effects were relatively structured and correctable, and the biological signal is strong.

Step 5 - UMAP visualisation

UMAP (Uniform Manifold Approximation and Projection) takes the Harmony-corrected embeddings and projects them into 2D for visualisation. Similar cells end up close together, different cells end up far apart:

sc.pp.neighbors(adata, use_rep='X_pca_harmony', n_neighbors=10)
sc.tl.umap(adata)

What we found - the biology

Finding 1: Five major cell populations in the gut

The UMAP revealed five well-separated clusters corresponding to the major cell types in the colon:

• Plasma cells (15,633 cells) - antibody factories
• Colon epithelial cells (12,347 cells) - the cells lining the gut wall
• T cells (12,128 cells) - immune soldiers
• Myeloid cells (3,771 cells) - macrophages and related immune cells
• Stromal cells (2,821 cells) - structural support cells

Finding 2: UC and CD have completely different cellular makeups

This is where the biology gets interesting. When we look at what proportion of cells belong to each type across disease groups:

Plasma cells:

• Normal: ~27%
• Crohn's disease: ~30%
• Ulcerative colitis: ~52%

UC has nearly double the plasma cell proportion of healthy tissue. Plasma cells make antibodies - this confirms that UC is primarily driven by antibody-mediated (humoral) immunity, a well-established finding that we independently reproduced from raw data.

Myeloid cells:

• Normal: ~6%
• Ulcerative colitis: ~8%
• Crohn's disease: ~13%

CD has more than double the myeloid cell proportion. Myeloid cells include macrophages - the cells responsible for the granulomatous (nodule-forming) inflammation that's the hallmark of Crohn's disease.

Epithelial cells:

• Normal: ~29%
• Crohn's disease: ~19%
• Ulcerative colitis: ~3%

This is striking. UC patients have almost no epithelial cells left - the gut lining is severely disrupted. This explains why UC causes such pronounced mucosal damage and bleeding.

Finding 3: The S100A8/S100A9 signature - a clinical biomarker reproduced at single-cell resolution

To find which specific genes are driving the myeloid difference in Crohn's disease, we ran a differential expression analysis comparing CD myeloid cells against normal myeloid cells using the Wilcoxon rank-sum test.

The top upregulated genes were:

Here's why S100A8 and S100A9 matter: together they form a protein called calprotectin. When gut inflammation is active, immune cells release calprotectin into the stool. Doctors measure this in a routine test called the faecal calprotectin test - one of the most common non-invasive ways to monitor IBD disease activity in clinic.

By analysing raw single-cell data, we independently identified the exact genes behind this clinical test - at the resolution of individual cells. That's the kind of validation that confirms the analysis is biologically meaningful, not just a statistical artefact.

What I learned - technical takeaways

1. Batch correction is non-negotiable for multi-sample studies. Without Harmony, the donor-to-donor variation would swamp the disease signal. Any "finding" could be explained by technical noise. Harmony is now a standard step in any multi-patient single-cell study.

2. Choosing highly variable genes with batch_key matters. If you find HVGs without accounting for batches, you risk selecting genes that are variable only because of technical differences between samples. Using batch_key='donor_id' ensures you're finding genes that are genuinely biologically variable.

3. Clinical relevance is the best validation. When your computational analysis independently reproduces a well-established clinical biomarker (faecal calprotectin), it gives you confidence that the pipeline is working correctly and the findings are real.

4. Public datasets are powerful. This entire analysis used freely available data from CZ CELLxGENE - no lab access required. The tools (Scanpy, harmonypy) are free and open source. Computational biology has an exceptionally low barrier to entry compared to wet lab science.

The full code

The complete annotated notebook is available on GitHub:
👉 github.com/Farhan89082/ibd-harmony-integration

The README includes all figures, biological interpretation, and instructions for reproducing the analysis from scratch.

What's next?

This project is part of a series of three single-cell RNA sequencing analyses I've built for my computational biology portfolio:

1. Alzheimer's Disease - microglial activation and mitochondrial dysfunction in human brain cells → github.com/Farhan89082/alzheimers-scrna-analysis
2. NSCLC Tumour Microenvironment - T cell exhaustion trajectories and macrophage polarisation in lung cancer → github.com/Farhan89082/nsclc-tumour-microenvironment
3. IBD Harmony Integration - this article

If you're a Python developer curious about getting into computational biology, or a biology student learning to code, I hope this walkthrough shows that the barrier is lower than it looks. The tools are excellent, the data is freely available, and the biology is genuinely fascinating.

Have questions about the analysis or want to discuss the methodology? Drop a comment below - I'd love to hear from you.