Entering edit mode
Hi all,
I am trying to correct batch effect using combat, could you please help if this is the correct approach or do I need to change the code? I also want to include Sentrix_position with Batch but not sure how combine two variable here?
library(sva)
> #phenotype file
> sample.info <-read.table("Sample.Info.csv", header=TRUE, sep=",")
> # first column as row names
>sample.info2<-sample.info[,-1]
> rownames(sample.info2)<-sample.info[,1]
>
> head(sample.info2)
Sentrix_ID Sentrix_Position Batch Recruitment Category Gender eGFR
DC541 2.03e+11 R06C01 Batch 15 US Control Male 141.3943
DC485 2.03e+11 R04C01 Batch 8 US Control Male 133.6376
DC490 2.03e+11 R08C01 Batch 8 US Control Male 133.1413
DC131 2.03e+11 R08C01 Batch 7 US Control Female 131.9288
DC574 2.03e+11 R03C01 Batch 16 US Control Female 130.6548
DC411 2.03e+11 R02C01 Batch 18 Belfast Control Male 127.5505
Age RRT CD8T CD4T NK Bcell Mono
DC541 29 FALSE 0.09559439 0.22490564 0.00000000 0.087392632 0.03744009
DC485 42 FALSE 0.04212238 0.12661890 0.04028556 0.024276752 0.09722832
DC490 39 FALSE 0.03568210 0.15196063 0.03766905 0.031379030 0.07066799
DC131 58 FALSE 0.10451144 0.04749711 0.03571969 0.004003534 0.06539260
DC574 24 FALSE 0.07358490 0.17676019 0.07706284 0.023125340 0.08663073
DC411 31 FALSE 0.01487569 0.11073860 0.01740834 0.065593039 0.05207325
Gran
DC541 0.4991278
DC485 0.6401493
DC490 0.6507270
DC131 0.7140364
DC574 0.4966720
DC411 0.7113963
> #read mehtylation data
>
> betaval<-read.table("methylation.data.txt", header=TRUE, stringsAsFactors=FALSE)
> # first column as row names
>
> betaval2<-betaval[,-1]
> rownames(betaval2)<-betaval[,1]
>
> #convert beta value into m value
> mval<-log2(betaval2/(1-betaval2))
>
modcombat <- model.matrix(~1, data=sample.info2)
data.cb <- ComBat(dat=as.matrix(mval), batch=sample.info2$Batch, mod=modcombat, par.prior=TRUE, prior.plots=FALSE)
Many thanks