limpa case study: two groups with Spectronaut quantification
First version 17 December 2025, Last revised 29 December 2025
Introduction
This case study analyses the HYE100 hybrid proteome samples from Navarro et al (2016). Two samples, A and B, were prepared by mixing human, yeast and E. coli tryptic peptides in specified proportions. The human peptides was obtained from HeLa cells, which is a cervical cancer cell line. Sample A was composed of 67% human, 30% yeast, and 3% E. coli, while sample B was composed of 67% human, 3% yeast, and 30% E. coli. The human proteins should not be differentially expressed between A and B, while yeast proteins should be 10-fold up-regulated in A, and E. coli proteins should be 10-fold up-regulated in B. Three technical replicates were prepared of A and B, to make a two group experiment with \(n=3\) samples in each group.
Demichev et al (2020) re-processed the mass spectrometry data using Spectronaut, and the Spectronaut report file with the precursor quantifications can be downloaded from https://doi.org/10.17605/OSF.IO/6G3UX.
We use this data for two purposes. First, we use it as an example analysis using Spectronaut data. Second, we explore the data a little more to illustrate the advantages and limitations of a mixed-species dataset such as this.
Example analysis
The peptide precursor intensities are provided in a tab-delimited report file output by Spectronaut, which we can read into limpa using readDIANN():
> library(limpa)Loading required package: limma> y.peptide <- readSpectronaut("Fig2-spectronaut.csv", sep = ",",
+ run.column = "R.FileName",
+ precursor.column = "EG.ModifiedPeptide",
+ qty.column = "FG.Quantity",
+ q.columns = "EG.Qvalue",
+ extra.columns = "PG.ProteinAccessions")
> dim(y.peptide)[1] 39338 6The samples are in two groups.
> Replicate <- factor(c(3,3,2,2,1,1))
> Group <- factor(c("B","A","B","A","B","A"))
> colnames(y.peptide) <- paste(Group, Replicate, sep=".")Add a species annotation column:
> PG <- y.peptide$genes$PG.ProteinAccessions
> n <- nchar(PG)
> y.peptide$genes$Species <- substring(PG,n-4L,n)
> table(y.peptide$genes$Species)
ECOLI HUMAN YEAS8
11384 16125 11829 The detection probability slope is estimated to be about 0.7.
Here we estimate DPC from the complete-normal model using dpcCN because it is robust to very large fold-changes in the data.
> dpcest <- dpcCN(y.peptide)Iter 1: dpc = -3.3195 0.7289
Iter 2: dpc = -3.1723 0.6979> plotDPC(dpcest)
Quantify protein expression:
> y <- dpcQuant(y.peptide, dpc=dpcest, protein.id="PG.ProteinAccessions")Estimating hyperparameters ...Quantifying proteins ...Proteins: 1000 Peptides: 8417Proteins: 2000 Peptides: 16265Proteins: 3000 Peptides: 22929Proteins: 4000 Peptides: 27590Proteins: 5000 Peptides: 34177Proteins: 6000 Peptides: 38047Proteins: 6507 Peptides: 39338The A and B samples separate very clearly, showing that there is considerable DE in the data:
> plotMDSUsingSEs(y)
Now we form the design matrix and proceed to differential expression analysis. First, limpa identifies a clear variance trend that will be used in the empirical Bayes modelling:
> design <- model.matrix(~Group)
> fit <- dpcDE(y, design, plot=TRUE)
Differential expression analysis identifies Yeast and E. coli samples as highly signficant:
> fit <- eBayes(fit)
> summary(dt <- decideTests(fit[,2])) GroupB
Down 1695
NotSig 3395
Up 1417> topTable(fit, coef=2, n=10) Species NPeptides PropObs logFC AveExpr t P.Value adj.P.Val B
1/tr|C8ZD01|C8ZD01_YEAS8 YEAS8 2 1.0000 -3.651 9.676 -98.29 1.446e-15 5.728e-12 26.43
1/tr|C8Z4U5|C8Z4U5_YEAS8 YEAS8 9 0.8889 -3.549 7.750 -91.87 2.740e-15 5.728e-12 25.78
2/tr|C8ZDM8|C8ZDM8_YEAS8/tr|C8ZGS8|C8ZGS8_YEAS8 YEAS8 2 1.0000 -3.439 9.826 -90.97 3.006e-15 5.728e-12 25.77
1/sp|P69908|DCEA_ECOLI ECOLI 2 1.0000 3.655 10.362 87.18 4.496e-15 5.728e-12 25.39
1/sp|P0A910|OMPA_ECOLI ECOLI 23 0.8696 3.122 8.528 84.68 5.916e-15 5.728e-12 25.14
1/sp|P63284|CLPB_ECOLI ECOLI 57 0.8450 3.238 6.711 83.26 6.945e-15 5.728e-12 24.93
1/sp|P0A853|TNAA_ECOLI ECOLI 55 0.8061 3.096 7.796 81.98 8.040e-15 5.728e-12 24.81
2/sp|P33898|G3P3_ECOLI/sp|P0A9B2|G3P1_ECOLI ECOLI 1 1.0000 3.509 12.890 81.71 8.296e-15 5.728e-12 24.75
1/sp|P0A9G6|ACEA_ECOLI ECOLI 39 0.7692 3.192 7.478 82.06 7.969e-15 5.728e-12 24.74
1/tr|C8Z3S1|C8Z3S1_YEAS8 YEAS8 8 0.8542 -3.579 8.255 -81.20 8.802e-15 5.728e-12 24.70A mean-difference plot shows that limpa recovers correct fold-changes for most protein, even single-peptides proteins and those at very low intensities:
> plotMD(fit, coef=2, status=fit$genes$Species, main="B vs A")
> abline(h=log2(10), lty=2)
> abline(h=0, lty=2)
> abline(h=-log2(10), lty=2)
The plot becomes clearer if we restrict to proteins with two or more precursors:
> i <- which(fit$genes$NPeptides>1)
> plotMD(fit[i,], coef=2, status=fit$genes$Species[i], main="B vs A")
> abline(h=log2(10),lty=2)
> abline(h=0, lty=2)
> abline(h=-log2(10),lty=2)
Proteins that are entirely missing in one group
One of the advantages of limpa is its ability to assess DE even for proteins that are entirely missing in one of the groups. We can compute the number of non-missing observations for each protein in each group:
> NObs.A <- rowSums(y$other$n.observations[,Group=="A"])
> NObs.B <- rowSums(y$other$n.observations[,Group=="B"])For this data, 499 of the proteins that are significantly up are missing entirely in the A samples, and 801 of the proteins that are significantly down are missing entirely in the B samples:
> table(DE=dt, pmin(NObs.A,5))
DE 0 1 2 3 4 5
-1 0 6 9 245 33 1402
0 314 204 195 968 97 1617
1 302 160 83 176 51 645> table(DE=dt, pmin(NObs.B,5))
DE 0 1 2 3 4 5
-1 450 178 98 253 71 645
0 467 223 169 887 67 1582
1 0 8 16 161 17 1215As as example, we pick out most significant of the latter proteins:
> PValue <- fit$p.value[,2]
> PGB <- names(which.min(PValue[NObs.B==0]))Plot the protein quantifications with standard errors:
> plotProtein(y, PGB, main=PGB)
Now plot the corresponding precursor log-intensities.
This yeast protein has 15 precursors, all of which are missing in all the B samples.
In plot, closed circles indicate observations and open circles indicate missing values.
For the purposes of the plot, the missing values are plotted 0.5 units below the minimum observation for the same precursor:
> plotPeptides(y.peptide[y.peptide$genes$PG.ProteinAccessions==PGB,],
+ main=paste(PGB,"precursors"))
Normalization seems unnecessary
We did not do any extra normalization on top of that already done by Spectronaut. The human proteins seem equally distributed between samples, so no extra normalization seems necessary:
> isH <- which(y$genes$Species == "HUMAN")
> boxplot(as.data.frame(y$E[isH,]),
+ ylab="log2-Intensity", main="Human proteins")
As expected, the yeast proteins are lower in B samples:
> isY <- which(y$genes$Species == "YEAS8")
> boxplot(as.data.frame(y$E[isY,]),
+ ylab="log2-Intensity", main="Yeast proteins")
while the E. coli proteins are lower in A samples:
> isE <- which(y$genes$Species == "ECOLI")
> boxplot(as.data.frame(y$E[isE,]),
+ ylab="log2-Intensity", main="E coli proteins")
Confirming that missingness is intensity-dependent
Now we explore the pattern of missing values at the precursor level:
> NMissing.A <- rowSums(is.na(y.peptide$E[,Group=="A"]))
> NMissing.B <- rowSums(is.na(y.peptide$E[,Group=="B"]))
> table(NMissing.A, y.peptide$genes$Species)
NMissing.A ECOLI HUMAN YEAS8
0 2243 14462 9097
1 802 722 1380
2 1531 601 1241
3 6808 340 111> table(NMissing.B, y.peptide$genes$Species)
NMissing.B ECOLI HUMAN YEAS8
0 9211 14464 1476
1 1124 780 648
2 973 574 1395
3 76 307 8310We see that human precursors have similar numbers of missing values in A and B samples, while E. coli precursors have far more missing values in A and yeast precursors have far more missing values in B samples. This confirms that the missing values are intensity-dependent, because the missingness pattern correlates with global changes in precursor intensities.
Limitations of this dataset
This dataset was created by Navarro et al (2016) and Demichev et al (2020) in order to benchmark peptide quantification algorithms. The dataset is, however, different from real biological datasets in a number of ways that make it unsuitable for benchmarking differential expression statistical methods. One issue that is that the DE proteins and peptides have systematically lower intensities (and consequently more missing values) than the non-DE proteins and peptides. Another issue is that log-fold-changes are very large, while the variation between replicates is very small and not representative of biological variation. Another issue is that the dataset contains false positives, whereby some of the human proteins are clearly differentially expressed, although the experimental design should have ensured that they were not. A striking example of a DE human protein is albumin, for which both precursors are clearly down in the B samples:
> i <- grep("ALBU_HUMAN",y.peptide$genes$PG.ProteinAccessions)
> plotPeptides(y.peptide[i,], main="ALBU_HUMAN precursors")
The fact that human albumin appears DE is an artefact of the dataset rather than an error of the statistical testing method.
Cited references
Demichev V, Messner CB, Vernardis SI, Lilley KS, Ralser M (2020). Spectronaut: neural networks and interference correction enable deep proteome coverage in high throughput. Nature Methods 17(1), 41-44.
Navarro P, Kuharev J, Gillet LC, Bernhardt OM, MacLean B, Rost HL, Tate SA, Tsou CC, Reiter L, Distler U, Rosenberger G, Perez-Riverol Y, Nesvizhskii AI, Aebersold R, Tenzer S (2016). A multicenter study benchmarks software tools for label-free proteome quantification. Nature Biotechnology, 34(11), 1130-1136.