Walkthrough (INTENDED)

Walkthrough (INTENDED)#

Instructor-complete version. Students receive a blank notebook plus the walkthrough prompts page for Colab cell-level AI. Each analysis step below is meant to be one prompt.

Learning objectives#

  • Filter low-abundance genes, then log2(TPM + 1)

  • Drop low-variance genes on the log2 scale

  • Compute gene-wise log2 fold changes, Welch t-tests, and Benjamini–Hochberg FDR

  • Summarize significant genes with Enrichr (MSigDB_Hallmark_2020)

Biological framing#

Acute HIV-1 infection triggers a large host transcriptional response in circulating immune cells. This walkthrough uses PBMC RNA-seq from GSE199911 (Parker et al., 2023): people sampled about one month into infection versus contemporaneous uninfected controls. The scientific goal is to find genes that differ between groups on a comparable expression scale, control the false-discovery rate across thousands of tests, and summarize the significant genes as Hallmark programs with Enrichr—the interferon and cell-cycle themes reported in the paper.

# Setting up the Colab environment. DO NOT EDIT!
try:
  from applied_biostats import setup_environment
except ImportError:
  !pip -q install applied-biostats-helper
  from applied_biostats import setup_environment
finally:
  grader = setup_environment('Module11_walkthrough')

!pip -q install gseapy==1.3.2

import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from scipy import stats
from statsmodels.stats.multitest import multipletests
import gseapy

Step 1. Load the staged tables#

Load the gse199911_tpm.csv.gz and gse199911_samples.csv sample sheet.

samples = pd.read_csv('gse199911_samples.csv')
tpm = pd.read_csv('gse199911_tpm.csv.gz')

print('Unit of analysis: sample (person-level PBMC RNA-seq)')
print('n samples:', len(samples))
print(samples['group'].value_counts())
print('tpm shape (genes x samples + symbol col):', tpm.shape)
tpm.head()

Step 2. Abundance filter#

Keep genes with TPM > 1 in at least 10 samples.

sample_cols = samples['sample_id'].tolist()
expr = tpm.set_index('gene_symbol')[sample_cols]

keep_abund = (expr > 1).sum(axis=1) >= 10
expr_f = expr.loc[keep_abund]
print('genes before abundance filter:', expr.shape[0])
print('genes after abundance filter:', expr_f.shape[0])
expr_f.head()

Step 3. Log2 transform#

Compute log2(TPM + 1) for the abundance-filtered matrix.

log2 = np.log2(expr_f + 1.0)
print('log2 shape:', log2.shape)
log2.head()

Step 4. Variance filter#

On the log2 matrix, drop genes in the bottom 20% of variance across samples.

gene_var = log2.var(axis=1)
var_cut = gene_var.quantile(0.20)
keep_var = gene_var > var_cut
log2 = log2.loc[keep_var]
print('variance cutoff (20th pct):', float(var_cut))
print('genes after variance filter:', log2.shape[0])

Step 5. Gene-wise log2 fold change#

For each gene, HIV mean minus control mean on the log2 scale.

hiv_cols = samples.loc[samples['group'] == 'hiv', 'sample_id'].tolist()
ctrl_cols = samples.loc[samples['group'] == 'control', 'sample_id'].tolist()

log2fc = log2[hiv_cols].mean(axis=1) - log2[ctrl_cols].mean(axis=1)
print('n genes with fold change:', len(log2fc))
log2fc.sort_values(ascending=False).head()

Step 6. Gene-wise Welch t-test#

Use a vectorized scipy.stats.ttest_ind down the genes to calculate a p-value for each gene.

hiv = log2[hiv_cols].to_numpy()
ctrl = log2[ctrl_cols].to_numpy()
t_stat, p_val = stats.ttest_ind(hiv, ctrl, axis=1, equal_var=False, nan_policy='omit')
print('n p-values:', len(p_val))
print('raw p < 0.05:', int(np.sum(p_val < 0.05)))

Step 7. Benjamini–Hochberg FDR#

Adjust the gene-wise p-values for multiple testing using the BH method and build the results table.

_reject, fdr, _, _ = multipletests(p_val, method='fdr_bh')

results = pd.DataFrame({
    'gene_symbol': log2.index,
    'log2fc': log2fc.to_numpy(),
    't': t_stat,
    'p': p_val,
    'fdr': fdr,
}).sort_values('fdr')

print('genes FDR < 0.05:', int((results['fdr'] < 0.05).sum()))
results.head(10)

Step 8. Build the Enrichr input list#

Keep genes with FDR < 0.05 and absolute log2 fold change > 1.

sig = results.loc[(results['fdr'] < 0.05) & (results['log2fc'].abs() > 1), 'gene_symbol'].tolist()
print('genes for Enrichr:', len(sig))
print(sig[:10])

Step 9. Enrichr (Hallmark)#

Run Enrichr on sig with MSigDB_Hallmark_2020. Use outdir=None.

enr = gseapy.enrichr(
    gene_list=sig,
    gene_sets='MSigDB_Hallmark_2020',
    organism='Human',
    outdir=None,
)
enrich_df = enr.results.sort_values('Adjusted P-value')
enrich_df.head(10)

Checkpoint#

You should now have:

  • results — one row per tested gene with log2fc, p, and fdr

  • enrich_df — Hallmark enrichment for the significant symbol list

The lab loads a pre-saved copy of the DE table (same abundance + variance filters) so it does not depend on this notebook having been run.

Contents