jtleek / sva-devel

28 stars 45 forks source link

combat for epic array data #59

Open bioinfonext opened 1 year ago

bioinfonext commented 1 year ago

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