I am working on TMT-labelled proteomics data (3 conditions in triplicates). Looking at the mean-variance relationship of my data on the peptide-level, it appears to be in good agreement with an additive-multiplicative error model, so I used vsn2 to transform the data. Now that the variance is stabilized, I believe I could directly use limma to asses the significance of changes on the peptide level. However, I would like to have a model for the significance on the protein level, which takes into account the number of detected peptides per protein (which is quite variable) and assigns higher significance to proteins with several peptides showing consistent changes than with just a single peptide showing the same change.  Is there a way to use limma (on the vsn-transformed data) to this end?

Limma fits a separate linear model to each "feature", which in your case is peptides. However, you could use a self-contained gene set test like roast (also part of limma) to combine the results for proteins with multiple peptides, giving a p-value for each protein representing the null hypothesis that none of its peptides are differentially expressed. You would treat each protein as a "gene set" consisting of all the peptides associated with it.

Also, if your proteomics data consists of discrete counts of peptides, you may want to try either voom or edgeR if their assumptions match your data, since both are designed precisely for analyzing count data. (roast is also available for these analyses)