EPIC/codescript/EPIC.raw.R

suppressPackageStartupMessages(library("optparse"))
suppressPackageStartupMessages(library("stats"))

# specify our desired options in a list
# by default OptionParser will add an help option equivalent to 
# make_option(c("-h", "--help"), action="store_true", default=FALSE, 
#               help="Show this help message and exit")

option_list <- list( 
  make_option(c("-p", "--prefix"), type="character",default="./",
              help="The output files prefix [default ./]"),
  make_option(c("-i", "--input"),type="character", default=NULL,
              help="The directory of input EPIC files. required!")
)

# get command line options, if help option encountered print help and exit,
# otherwise if options not found on command line then set defaults, 
opt <- parse_args(OptionParser(option_list=option_list))

# load libraries

library("minfi")
library("IlluminaHumanMethylationEPICmanifest")
library("IlluminaHumanMethylationEPICanno.ilm10b4.hg19")

## 1. data import
targets <- read.metharray.sheet(opt$input)
rgSet <- read.metharray.exp(targets = targets)
targets$ID <- paste(targets$Sample_Group,targets$Sample_Name,sep=".")
sampleNames(rgSet) <- targets$ID
#phenoData <- pData(rgSet)
message("data import:finished")

## 2. Quality control
# 2.1 qc report by minfi
# qcReport(rgSet, sampNames=targets$ID, sampGroups=targets$Sample_Group,pdf="qcReport.pdf")
          
# 2.2 data filtering
# get detected p value
# a. remove samples with average p value less than 0.05
detP <- detectionP(rgSet)
message("calculate p value:finished")

# raw
mSetRaw <- preprocessRaw(rgSet)
raw_mVals <- getM(mSetRaw)
raw_bVals <- getBeta(mSetRaw)
message("m value and beta value output:finished")

# normalization
MSet.noob <- preprocessNoob(rgSet)
RSet <- ratioConvert(MSet.noob, what = "both", keepCN = TRUE)
bVals <- getBeta(RSet)
mVals <- getM(RSet)

# write output
message("saving R data")
rdata_filename = paste(opt$prefix, '.RData',sep="")
save(rgSet, targets, detP, file = rdata_filename)
message("writing raw table")
m_raw_filename = paste(opt$prefix,'.raw.mVal.txt',sep="")
write.table(raw_mVals,m_raw_filename,col.names = T,row.names = T,sep="\t",quote=F)
b_raw_filename = paste(opt$prefix,'.raw.bVal.txt',sep="")
write.table(raw_bVals,b_raw_filename,col.names = T,row.names = T,sep="\t",quote=F)
message("writing filtered table")
m_filename = paste(opt$prefix,'.ssnoob.mVal.txt',sep="")
write.table(mVals,m_filename,col.names = T,row.names = T,sep="\t",quote=F)
b_filename = paste(opt$prefix,'.ssnoob.bVal.txt',sep="")
write.table(bVals,b_filename,col.names = T,row.names = T,sep="\t",quote=F)