Blocked design and Surrogate Variables Analsysis(SVA)
1
2
Entering edit mode
kodream ▴ 20
@kodream-6952
Last seen 7.6 years ago
United States

I would like to try to use SVA with a block design.  I've been using SVA in separate non-blocked experiments with similar data, with known artifacts which SVA removes quite well, however I am not quite sure how to use it with blocking.  SVA doesn't seem to be aware of blocking, so I just run it first, then added the surrogate vectors to the design, before I run duplicate correlation.  This is what I have so far, it seems like it will miss something:

svobj<-sva(exprs(eset),modMMR, modMMR0, method="irw", n.sv=n.sv)

modsMMR <-cbind(modMMR,svobj$sv)

#my version of duplicateCorrelation that uses the parallel foreach instead of for
corfit <- parallel_dupCor(eset,modsMMR,block=pData(eset)$Pair)
corfit$consensus

lfit <- lmFit(exprs(eset), design=modsMMR, block = pData(eset)$Pair, cor = corfit$consensus)

fit<-eBayes(lfit)
limma sva block • 2.2k views
ADD COMMENT
0
Entering edit mode

I know this post is very old but could you post your code that runs duplicateCorrelation in parallel? I have a very big dataset and can't quite figure out how to update the rho variable in parallel... Thanks!

ADD REPLY
1
Entering edit mode
Jeff Leek ▴ 650
@jeff-leek-5015
Last seen 3.7 years ago
United States

To handle the blocking, you need to include the blocking variables in both mod and mod0, then you should be able to calculate the surrogate variables and include them in the downstream analysis at the lmFit stage in the model matrix. If you show modMMR and modMMR0, I may be able to help you construct the exact model matrices. 

ADD COMMENT
1
Entering edit mode

Hi Jeff,

I have the same question regarding a blocked design. I am applying SVA to a large EPIC methylation array dataset with paired samples (patients/controls with enrollment and follow up). I am currently incorporating the block at the lmfit stage as the original poster, as I wasn't quite sure how to incorporate it in the mod correctly. My code is provided below. If you incorporate the blocking variable in the mod, then do you exclude corfit$consensus and block at the lmfit stage?  Any help you could provide me with would be very greatly appreciated!

Best,

Katie

 

#Followup and Enrollment Pooled

############################

mval <- edata

mod <- model.matrix(~as.factor(SubjectType), data=pheno)

mod0 <- model.matrix(~1, data=pheno)

sva.results <- sva(mval, mod, mod0)

 

sv = sva.results$sv

modSv = cbind(mod,sv)

 

fit <- lmFit(edata,modSv, block=pheno$patientID,correlation=corfit$consensus)

 

contrast.matrix <- cbind("PDvsControl"=c(0,1,rep(0,sva.results$n.sv)))

fitContrasts = contrasts.fit(fit,contrast.matrix)

eb = eBayes(fitContrasts)

topTableF(eb, adjust="BH")

summary(decideTests(eb))

 

 

#Followup and enrollment separated by condition

########################################

mval <- edata

condition_visit <- paste(pheno$SubjectType, pheno$Visit_Name2, sep="_" ) #4 levels: PD_Enrollment, PD_Followup, HC_Enrollment, HC_Followup

 

mod <- model.matrix(~0+as.factor(condition_visit), data=pheno)

mod0 <- model.matrix(~1, data=pheno)

sva.results <- sva(mval, mod, mod0)

 

sv = sva.results$sv

modSv = cbind(mod,sv)

 

fit <- lmFit(edata,modSv, block=pheno$patientID,correlation=corfit$consensus)

 

contrast.matrix <- cbind("PDvsHC_Enrollment"=c(-1,0,1,0,rep(0,sva.results$n.sv)),

                        "PDvsHC_Followup"=c(0,-1,0,1,rep(0,sva.results$n.sv)),

                         "PD_FollowupvsEnrollment" = c(0,0,-1,1,rep(0,sva.results$n.sv)),

                         "HC_FollowupvsEnrollment" =c(-1,1,0,0,rep(0,sva.results$n.sv)))

fitContrasts = contrasts.fit(fit,contrast.matrix)

eb = eBayes(fitContrasts)

topTableF(eb, adjust="BH")

summary(decideTests(eb))

ADD REPLY

Login before adding your answer.

Traffic: 668 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6