ChAMP Methylation Pipeline

library(ChAMP)

idat_folder ="./Idats/"

working_dir = "./Working_Dir"

array_type = "450K" # CAN BE EITHER EPIC OF 450K

core_num = 4 # Number of Cores used for a few functions

group_column = "Sample_Group" # NAME OF COLUMN FROM SAMPLE SHEET CONTAINING THE GROUP LABEL (0 - 1 or 2 levels)

# LOAD

#?champ.load

setwd(working_dir)

myLoad <- champ.load(directory = idat_folder,
                     method="ChAMP",
                     methValue="B",
                     autoimpute=FALSE,
                     filterDetP=TRUE,
                     ProbeCutoff=0,
                     SampleCutoff=0.1,
                     detPcut=0.01,
                     filterBeads=TRUE,
                     beadCutoff=0.05,
                     filterNoCG=TRUE,
                     filterSNPs=FALSE,
                     population=NULL,
                     filterMultiHit=TRUE,
                     filterXY=FALSE,
                     force=FALSE,
                     arraytype=array_type)
#?champ.QC

setwd(working_dir)

champ.QC(beta = myLoad$beta,
  pheno=myLoad$pd[,group_column],
  mdsPlot=TRUE,
  densityPlot=TRUE,
  dendrogram=TRUE,
  PDFplot=TRUE,
  Rplot=TRUE,
  Feature.sel="None",
  resultsDir="./CHAMP_QCimages/")

#?champ.norm
# USING BMIQ Normalisation (good enough for most cases)
myNorm <- champ.norm(beta=myLoad$beta,
                     rgSet=myLoad$rgSet,
                     mset=myLoad$mset,
                     resultsDir="./CHAMP_Normalization/",
                     method="BMIQ",
                     plotBMIQ=FALSE,
                     arraytype=array_type,
                     cores=core_num)
#?champ.SVD
champ.SVD(beta = myNorm,
          rgSet=NULL,
          pd=myLoad$pd,
          RGEffect=FALSE,
          PDFplot=TRUE,
          Rplot=TRUE,
          resultsDir="./CHAMP_SVDimages/")

#?champ.runCombat
myCombat <- champ.runCombat(beta=myNorm,
                            pd=myLoad$pd,
                            variablename=group_column,
                            batchname=c("Slide", "Array"),
                            logitTrans=TRUE)

champ.SVD(beta = myCombat,
          rgSet=NULL,
          pd=myLoad$pd,
          RGEffect=FALSE,
          PDFplot=TRUE,
          Rplot=TRUE,
          resultsDir="./CHAMP_SVDimages_after_Combat/")

# Replace myNorm by myCombat

myNorm <- myCombat

champ.QC(beta = myNorm,
         pheno=myLoad$pd[, group_column],
         mdsPlot=TRUE,
         densityPlot=TRUE,
         dendrogram=TRUE,
         PDFplot=TRUE,
         Rplot=TRUE,
         Feature.sel="None",
         resultsDir="./CHAMP_QCimages_after_Norm_and_Combat/")
# EXPORT 

write.csv(myNorm, "beta_matrix.csv")
saveRDS(myNorm, "beta_matrix.RDS")