Module 11 walkthrough prompts#
Companion to the Module 11 walkthrough. Use each Point to write your own Colab cell-level AI prompt. Run the generated code and check the printed objects before the next step. Open the hint only if you want a short sample prompt.
Setup#
The notebook already installs applied-biostats-helper and gseapy==1.3.2 and imports numpy, pandas, seaborn, matplotlib, scipy.stats, statsmodels BH correction, and gseapy. Do not ask the AI to reinstall those unless an import fails.
Step 1. Load the staged tables#
Point: Load gse199911_samples.csv into samples and gse199911_tpm.csv.gz into tpm. Confirm the design at the sample level (PBMC RNA-seq): print n samples, samples['group'].value_counts(), and tpm.shape, then show tpm.head().
Hint
Load
samplesandtpmfrom the staged CSVs. Print n, group counts, shape, andtpm.head().
Step 2. Abundance filter#
Point: From tpm, set gene symbols as the index and keep only columns in samples['sample_id']. Keep genes with TPM > 1 in at least 10 samples. Store the filtered matrix as expr_f. Print gene counts before and after, then expr_f.head().
Hint
Index by gene, align to sample IDs, keep genes with TPM > 1 in ≥10 samples as
expr_f, and print before/after n genes.
Step 3. Log2 transform#
Point: Create log2 as log2(TPM + 1) from expr_f. Print shape and log2.head().
Hint
Set
log2 = log2(expr_f + 1); print shape and head.
Step 4. Variance filter#
Point: On log2, compute per-gene variance across samples. Drop genes in the bottom 20% of variance. Overwrite log2 with the kept genes. Print the variance cutoff and the number of genes remaining.
Hint
Drop the bottom 20% variance genes from
log2; print cutoff and n remaining.
Step 5. Gene-wise log2 fold change#
Point: From samples, get sample_id lists for hiv and control. For each gene in log2, compute HIV mean minus control mean. Store the Series as log2fc. Print how many genes have a fold change and show the top of log2fc sorted descending.
Hint
Build
log2fcas HIV − control means onlog2; show the top sorted values.
Step 6. Gene-wise Welch t-test#
Point: Using the HIV and control columns of log2, run a vectorized Welch t-test with scipy.stats.ttest_ind(..., equal_var=False) along genes (avoid pingouin loops). Store t_stat and p_val. Print the number of p-values and how many raw p-values are < 0.05.
Hint
Vectorized Welch
ttest_indonlog2→t_stat,p_val; count raw p < 0.05.
Step 7. Benjamini–Hochberg FDR#
Point: Adjust p_val with Benjamini–Hochberg via statsmodels.stats.multitest.multipletests. Build DataFrame results with columns gene_symbol, log2fc, t, p, and fdr, sorted by fdr. Print how many genes have FDR < 0.05 and show results.head(10).
Hint
BH-adjust
p_valintoresults(gene_symbol,log2fc,t,p,fdr); count FDR < 0.05.
Step 8. Build the Enrichr input list#
Point: From results, make list sig of gene symbols with FDR < 0.05 and |log2FC| > 1. Print len(sig) and the first 10 symbols.
Hint
Make
sigfrom FDR < 0.05 and |log2FC| > 1; print length and first 10.
Step 9. Enrichr (Hallmark)#
Point: Run gseapy.enrichr on sig with gene set library MSigDB_Hallmark_2020, organism Human, and outdir=None. Store the table as enrich_df sorted by adjusted p-value and show the top 10 rows.
Hint
gseapy.enrichronsigwith Hallmark 2020 →enrich_df; show top 10 by adj-p.