Identify proteins with significantly different abundance between experimental conditions. The central decisions are how to handle missing values (model the left-censored MNAR dropout, do not impute it), how to moderate per-protein variance at the small sample sizes proteomics uses, and at what level to test (protein summary vs feature/peptide). The skill covers limma, DEqMS, proDA, msqrob2, MSstats, and a Python Welch+BH fallback, plus minimum-fold-change testing and fold-change shrinkage.
pip install numpy pandas scipy statsmodels
BiocManager::install(c("limma", "DEqMS", "proDA", "msqrob2", "MSstats", "ashr"))
Tell your AI agent what you want to do: - "Find differentially abundant proteins between treatment and control in my intensity matrix" - "Run limma with empirical-Bayes moderation on my small-n protein data" - "Use DEqMS because I have PSM counts per protein from a TMT experiment" - "My label-free data has 30% missing values and some on/off proteins -- test it without imputing" - "Test for at least a 1.5-fold change instead of just nonzero, without inflating FDR"
"I have 4 replicates per group of protein-level LFQ intensities. Pick and run the right differential test."
"Analyze my TMT proteomics data for differential abundance. I have PSM counts per protein, so use DEqMS, and remember it is multi-batch."
"My label-free data has many proteins detected in one group but missing in the other. Test these honestly instead of imputing."
"Set up a limma model with condition and batch as covariates -- do not remove the batch effect before testing."
"Some proteins are missing because they are below the detection limit. Use a method that models this dropout."
"Test whether fold changes exceed 1.5-fold using treat and topTreat, not a post-hoc filter."
"Report raw fold changes for GSEA and shrunken estimates for the figure."
removeBatchEffect before testing) and define contrasts.trend=TRUE, robust=TRUE; DEqMS count prior when available).treat()+topTreat() rather than a FC+significance double filter.| Method | Best for | Key advantage |
|---|---|---|
| limma | Small n (3-5), protein-level summaries | Borrows variance across proteins via empirical Bayes |
| DEqMS | PSM/peptide counts available | Prior keyed on quantification depth; dominates limma-trend |
| proDA | Label-free with extensive MNAR missing values | Models dropout in the likelihood; no imputation |
| msqrob2 | Outlier-peptide / unbalanced coverage | Peptide-level robust ridge; keeps feature df |
| MSstats | Technical replicates, nested/labeled designs | Feature-level mixed models |
| Welch t-test + BH | Large n (>10), Python-only | Simple; no moderation, unsuitable at small n |
The dominant statistical problem is missingness, and in label-free MS it is left-censored MNAR -- missing because the intensity is low. Imputing it injects directional bias: - Perseus/MaxQuant downshift manufactures systematic false positives (the volcano "anchor/wing" artifact) by giving on/off proteins an inflated mean offset and a collapsed within-group variance. - kNN/mean imputation is mean-reverting and compresses real down-regulation. - The correct approach models the dropout: proDA (probabilistic dropout), msqrob2, or MSstats with AFT censoring. Report on/off proteins as "undetected in group X", not as a giant fold change.
eBayes(trend = TRUE, robust = TRUE) -- the trend is effectively mandatory for label-free intensity data.equal_var=False in scipy.stats.ttest_ind (the default is Student's), and pass method='fdr_bh' to multipletests (the default is Holm-Sidak).topTable/topTreat return adj.P.Val (there is no $FDR column); topTreat omits B.treat()+topTreat(); never topTable(lfc=...) nor a post-hoc abs(logFC)>1 & adj.P.Val<0.05 double filter (it inflates realized FDR above 50%).removeBatchEffect() only for PCA/visualization, never as input to lmFit.