Dear all,

I have a complex bulk RNA-seq dataset (from a human tissue) with multiple covariates (sex, city, age), different cell type composition, repeated measures (same patient sequenced 4 times --> 4 different locations of the tissue). I am applying sva to correct for different cell type compositions and other unwanted variation due to other factors. limma-voom for normalization and differential expression across groups.This is the code:

y=DGEList(counts=counts)
A<-rowSums(y$counts)
isexpr<-A>500
y=y[isexpr,keep.lib.size=FALSE]
dim(y)
y=calcNormFactors(y)
group=factor(info_new$group)
sex=factor(info_new$sex)
city=factor(info_new$city)
age=cut(info_new$age,4)
levels(age)=c(age1,age2,age3,age4)
mod <- model.matrix(~sex+age+city+group)
mod0 <- model.matrix(~sex+age+city)
v=voom(y,mod)
n.sv = num.sv(v$E,mod,method="leek")
sva_obj <- sva(v$E, mod, mod0, n.sv=n.sv)
mod1 <- model.matrix(~sex+age+city+group+sva_obj$sv)
v <- voom(counts=y, design = mod1)
fit <- lmFit(v, mod1) 
fit2<- eBayes(fit)
summary(decideTests(fit2))

My question is, does sva remove the subject effect as well and it is not necessary to apply duplicateCorrelation?

Here a simple version of the sample information (I have 92 RNA-seq samples, 63 patients, 7 cities, 28 groups, control group only from loc4):

patient

city

sex

age

group

1

city1

M

52

disease1_tissue_loc1

1

city1

M

52

disease1_tissue_loc2

1

city1

M

52

disease1_tissue_loc3

1

city1

M

52

disease1_tissue_loc4

2

city1

F

28

disease2_tissue_loc1

2

city1

F

28

disease2_tissue_loc2

2

city1

F

28

disease2_tissue_loc3

2

city1

F

28

disease2_tissue_loc4

3

city1

M

38

disease2_tissue_loc1

3

city1

M

38

disease2_tissue_loc2

3

city1

M

38

disease2_tissue_loc3

3

city1

M

38

disease2_tissue_loc4

4

city1

F

61

disease1_tissue_loc1

4

city2

F

61

disease1_tissue_loc2

4

city2

F

61

disease1_tissue_loc3

4

city2

F

61

disease1_tissue_loc4

5

city2

F

29

control_tissue_loc4

6

city2

F

53

control_tissue_loc4

7

city2

M

43

control_tissue_loc4

In principle, if you don't include 'patient' in the full model (mod) it will be captured by SVA (assuming that it has a significant contribution to the overall variation observed). However, there are a few points to consider:

1- You have already included sex, age, and city in the null model during estimation of surrogate variables, therefore including them in the design formula alongside surrogate variables is redundant. What you want is:

mod1 <- model.matrix(~group+sva_obj$sv)

2- The goal of SVA is to estimate unknown batch that you have no knowledge of, but have reasons to believe they exist, for example after exploratory analysis using PCA or hierarchical clustering. If your only goal is to adjust for known factors in differential expression analysis, you can directly include them in your design formula as covariates.

3- The magic of voom lies in precision weights that are used to model the mean-variance trend, which you are not taking any advantage of in your estimation of surrogate variables. By using v$E you're only taking log-transformed counts. Whether a simple log-transformation is sufficient to adjust for the mean-variance trend is another point that you can check. If the answer is no, you may want to give sva-seq a shot, which directly takes RNAseq counts as input.

Hope these help.