Tissue-specific vulnerabilities of Foxp3-driven transferrin receptor deletion

Analysis performed by: Anton Zhelonkin

Sample Information

Sample ID	WT or KO	Treg or Tcon	Spleen or Lung	Frozen or Fresh	Mouse ID
10041-KV-1	WT	Tcon	Lung	Cryopreserved	Mouse 886-all
10041-KV-2	WT	Treg	Lung	Cryopreserved	Mouse 886-all
10041-KV-3	KO	Tcon	Lung	Cryopreserved	Mouse 5+6
10041-KV-4	KO	Treg	Lung	Cryopreserved	Mouse 5+6
10041-KV-5	WT	Tcon	Lung	Cryopreserved	Mouse 887
10041-KV-6	WT	Treg	Lung	Cryopreserved	Mouse 887
10041-KV-7	KO	Tcon	Lung	Cryopreserved	Mouse 642-2+339-2
10041-KV-8	KO	Treg	Lung	Cryopreserved	Mouse 642-2+339-2
10041-KV-9	KO	Tcon	Lung	Cryopreserved	Mouse 2+3
10041-KV-10	KO	Treg	Lung	Cryopreserved	Mouse 2+3
10041-KV-11	KO	Treg	Lung	Fresh Cells	Mouse 2+4
10041-KV-12	WT	Treg	Lung	Fresh Cells	Mouse 5
10041-KV-13	WT	Treg	Lung	Fresh Cells	Mouse 6
10041-KV-14	WT	Treg	Spleen	Cryopreserved	Mouse 886-1
10041-KV-15	WT	Treg	Spleen	Cryopreserved	Mouse 886-2
10041-KV-16	WT	Treg	Spleen	Cryopreserved	Mouse 886-3
10041-KV-17	KO	Treg	Spleen	Cryopreserved	Mouse 642-2
10041-KV-18	KO	Treg	Spleen	Cryopreserved	Mouse 339-2
10041-KV-19	WT	Treg	Spleen	Fresh Cells	Mouse 5
10041-KV-20	KO	Treg	Spleen	Fresh Cells	Mouse 4

Note: 10041-KV-1 reads were corrupted and removed prior to analysis.

Methods

Reads processing, QC and alignment

Paired-end FASTQ files for each sample were aligned to the mouse genome GRCm39 RefSeq assembly GCF_000001635.27 (Release Feb 2024). Reads were aligned and coordinate-sorted using STAR¹ in two-pass mode (STAR --twopassMode Basic). Each sample underwent pre- and post-alignment quality control with FastQC² 0.12.1, Picard³ 3.2.0, RSeQC⁴ 5.0.4, and MultiQC⁵ 1.25. Adapters and low-quality bases were trimmed with Trimmomatic⁶ 0.39. PCR duplicates were marked with Picard tools, and the downstream impact of deduplication was explicitly evaluated. Gene-level counts were quantified with featureCounts⁷ from the Subread 2.0.2 release.

Downstream Analysis

Bulk RNA-seq count data were analyzed in R⁸ 4.4.2 using Bioconductor packages⁹. Differential expression analysis was performed with edgeR¹⁰,¹¹ 4.2.2 and limma-voom¹² 3.62.2.

Deduplication Assessment

Deduplicated and non-deduplicated count matrices were compared to evaluate the impact of PCR duplicate removal. Deduplication revealed no added benefits, and so we decided, following best practices for bulk RNA-seq without unique molecular identifiers¹³,¹⁴ to proceed with non-deduplicated counts for all downstream analyses.

Statistical Modelling Framework

Analyses were restricted to regulatory T cell (Treg) samples; conventional T cell (Tcon) samples were excluded. Lowly expressed genes were removed with filterByExpr(), and library sizes were normalized using the Trimmed Mean of M-values (TMM) method¹⁵. Linear modeling was performed within the limma framework¹⁶, adopting a two-model comparison to evaluate processing effects (fresh vs cryopreserved).

Model 1: Cryopreserved-only analysis

To eliminate processing-related confounding, Model 1 restricted the dataset to cryopreserved Treg samples and fit a no-intercept organ × genotype interaction model using limma-voom:

y_g = X β_g + ε_g,
X = [Lung_KO, Lung_WT, Spleen_KO, Spleen_WT]

Each indicator column equals 1 for samples in the specified organ × genotype group (0 otherwise), enabling KO vs WT comparisons within organ without processing covariates.

Model 2: Batch-corrected full dataset

Model 2 leveraged all Treg samples while accounting for processing differences via a no-intercept design with four organ × genotype indicators and a binary fresh covariate:

y_g = X β_g + ε_g,
X = [Lung_KO, Lung_WT, Spleen_KO, Spleen_WT, fresh]

The fresh term equals 1 for fresh samples and 0 for cryopreserved samples, corresponding in R to ~ 0 + organ:genotype + fresh.

Model Fitting and Comparison

Linear models were fit with voomLmFit() without sample weights. Empirical sample quality weighting¹⁷ was evaluated (voomLmFit(..., sample.weights = TRUE)) to capture increased variability in KO samples. However, weighting algorithm downweighted KO samples (KV-0017, KV-0018) showing the largest genotype effects, suppressing biologically meaningful fold changes. Fits without weights preserved the signal of interest, so we procedeed with sample.weights = FALSE. Empirical Bayes moderation was performed using eBayes() with robust = TRUE¹⁸ for improved variance estimation, and removeBatchEffect()¹⁶ was used only for visualization (PCA/MDS).

Model diagnostics included Bland-Altman plots²⁰, MA plots, mean-variance trends, multidimensional scaling, log fold-change correlations, and differential expression overlap. Key comparisons between models highlighted:

• Fold-change agreement: Bland-Altman analysis showed agreement for |mean logFC| < 2, while Model 1 inflated fold-changes for strongly dysregulated genes (|mean logFC| > 5).
• Variance stabilization: Model 1 exhibited extreme variance inflation at high expression levels; Model 2 stabilized variance by modeling processing effects.
• Statistical power: Model 1 seemed to suffer reduced power with many genes untestable (p-value = 1) probably due to zero estimated variance; Model 2 estimated accurate p-values for all genes, with no genes falling back to p=1.
• Differential expression robustness: Only one-third of Model 1 hits remained significant in Model 2, probably indicating that throwing out replicates inflated estimates in the cryopreserved-only analysis. No genes changed direction between models, demonstrating that batch correction attenuated but did not invert biological signals. For the KO_vs_WT in (i)the spleen, (ii) in the lung, and (iii) for the difference contrasts between (KO_vs_WT in the spleen) vs (KO_vs_WT in the lung) there were uniquely identified model 2 specific DEGs in the amount of (i) 38, (ii) 111, (iii) 18, respectively.

See ./3_Results/Followup/Model_fit for the results of model comparison.

Contrasts Design

Both models tested the same pre-specified contrasts:

• Tissue-specific KO effects
  • Spleen: KO_vs_WT_Spleen = Spleen_KO − Spleen_WT
  • Lung: KO_vs_WT_Lung = Lung_KO − Lung_WT
• Genotype × tissue interaction
  • Diff_Effect = (Spleen_KO − Spleen_WT) − (Lung_KO − Lung_WT)
• Overall KO effect
  • KO_Overall = (Lung_KO + Spleen_KO)/2 − (Lung_WT + Spleen_WT)/2

Multiple testing was controlled using the Benjamini–Hochberg FDR procedure¹⁹ with p < 0.05.

Gene Set Enrichment Analysis

Gene set enrichment analysis²¹ was performed with clusterProfiler²²⁻²⁴ 4.14.6 using fgsea²⁵ for fast preranked enrichment. All genes passing the low-expression filter were ranked by moderated t-statistic to minimize pathway bias²⁶,²⁷. Pathways from MSigDB²⁸ (via msigdbr 10.0.2) were tested, including Hallmark²⁹ (H), Gene Ontology Biological Processes³⁰,³¹ (C5:GO:BP), Gene Ontology Molecular Functions (C5:GO:MF), KEGG³² (C2:CP:KEGG), Reactome³³ (C2:CP:REACTOME), WikiPathways³⁴ (C2:CP:WIKIPATHWAYS), and BioCarta³⁵ (C2:CP:BIOCARTA). GSEA used 100,000 permutations for robust p-value estimation. Pathways with FDR-adjusted p-value < 0.05 were deemed significant; dot plots outline significant pathways in black. Custom visualisation scripts were used, and were called from the main analysis .R file, the scripts are available at their dedicated versioned GitHub repo

Data and Code Availability

All raw and processed data, analysis scripts, and results are available at GitHub repo.

References

1. Dobin A, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29:15-21.
2. Andrews S. FastQC: a quality control tool for high throughput sequence data. 2010. Available online at https://www.bioinformatics.babraham.ac.uk/projects/fastqc/.
3. Broad Institute. Picard Toolkit. http://broadinstitute.github.io/picard/.
4. Wang L, et al. RSeQC: quality control of RNA-seq experiments. Bioinformatics. 2012;28:2184-2185.
5. Ewels P, et al. MultiQC: summarize analysis results for multiple tools and samples. Bioinformatics. 2016;32:3047-3048.
6. Bolger AM, et al. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114-2120.
7. Liao Y, et al. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30:923-930.
8. R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing; 2024.
9. Huber W, et al. Orchestrating high-throughput genomic analysis with Bioconductor. Nat Methods. 2015;12:115-121.
10. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139-140.
11. McCarthy DJ, Chen Y, Smyth GK. Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Res. 2012;40:4288-4297.
12. Law CW, et al. voom: precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biol. 2014;15:R29.
13. Parekh S, et al. The impact of amplification on differential expression analyses by RNA-seq. Sci Rep. 2016;6:25533.
14. Fu Y, et al. Elimination of PCR duplicates in RNA-seq and small RNA-seq using unique molecular identifiers. BMC Genomics. 2018;19:531.
15. Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 2010;11:R25.
16. Smyth GK. Linear models and empirical Bayes methods for assessing differential expression in microarray experiments. Stat Appl Genet Mol Biol. 2004;3:Article3.
17. Liu R, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43:e47.
18. Phipson B, et al. Robust hyperparameter estimation protects against hypervariable genes and improves power to detect differential expression. Biostatistics. 2016;17:364-376.
19. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B. 1995;57:289-300.
20. Bland JM, Altman DG. Statistical methods for assessing agreement between two methods of clinical measurement. Lancet. 1986;1:307-310.
21. Subramanian A, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA. 2005;102:15545-15550.
22. Yu G, et al. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16:284-287.
23. Wu T, et al. clusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innovations. 2021;2:100141.
24. Yu G, He QY. ReactomePA: an R/Bioconductor package for reactome pathway analysis and visualization. Mol Biosyst. 2016;12:477-479.
25. Korotkevich G, et al. Fast gene set enrichment analysis. bioRxiv. 2021. doi:10.1101/060012.
26. Goeman JJ, Bühlmann P. Analyzing gene expression data in terms of gene sets: methodological issues. Bioinformatics. 2007;23:980-987.
27. Tamayo P, et al. The limitations of simple gene set analysis for gene expression data. PLoS Comput Biol. 2016;12:e1005147.
28. Liberzon A, et al. Molecular signatures database (MSigDB) 3.0. Bioinformatics. 2011;27:1739-1740.
29. Liberzon A, et al. The Molecular Signatures Database Hallmark gene set collection. Cell Syst. 2015;1:417-425.
30. Ashburner M, et al. Gene ontology: tool for the unification of biology. Nat Genet. 2000;25:25-29.
31. The Gene Ontology Consortium. The Gene Ontology resource: enriching a GOld mine. Nucleic Acids Res. 2021;49:D325-D334.
32. Kanehisa M, et al. KEGG: integrating viruses and cellular organisms. Nucleic Acids Res. 2023;51:D678-D684.
33. Gillespie M, et al. The reactome pathway knowledgebase 2022. Nucleic Acids Res. 2022;50:D687-D692.
34. Slenter DN, et al. WikiPathways: a multifaceted pathway database bridging metabolomics to other omics research. Nucleic Acids Res. 2018;46:D661-D667.
35. Nishimura D. BioCarta. Biotech Software & Internet Report. 2001;2:117-120.