Yes, we have several biomarkers, and a proteomic set, so we are characterizing each biomarker by it's proteomic association(s). The use for surrogate variables in this case would be to account for unobserved batch effects across samples. So we can generate svas for each biomarker independently (with the whole proteomic set as the expression set).

But the thought occurred that in having multiple biomarkers tested on the same samples, presumably true unobserved sampling batch effects might be those that appear across multiple independently measured biomarkers. ie for i in [biomarker[1], ..., biomarker[n]] --> biomarker[i] ~ covs + shared_svas + protein

To generate the svas, typically it would be

mod0 <- model.matrix('~ covs', data = pheno_data) 
mod <- model.matrix('~ covs + biomarker[i]', data = pheno_data) 
svobj <- sva(prot_data, mod, mod0) # svas for i only

but putting mod <- model.matrix('~ covs + biomarker[1] + ... + biomarker[n]', data = pheno_data) would... generate svas that adjust for latent variables that occur between any individual or combo of biomarkers and the prots (like a union), without telling me which biomarker it varies across? Or generate svas that adjust for latent variables represented in the association of all biomarkers and the prots (like an intersection), i.e only latent variation if it shows up in all the biomarkers, not just one or a few?

I would prefer the latter, as you would just have a set of svas that could be shared for all biomarker ~ covs + shared_svas + prot_data comparisons instead of sva[i][1], sva[i][2], ... for each biomarker[i].