|
|
|
|
|
|
|
|
|
|
|
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") |
|
|
|
|
|
|
|
|
|
|
|
# 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) |