You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

5 yıl önce
5 yıl önce
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869
  1. suppressPackageStartupMessages(library("optparse"))
  2. suppressPackageStartupMessages(library("stats"))
  3. # specify our desired options in a list
  4. # by default OptionParser will add an help option equivalent to
  5. # make_option(c("-h", "--help"), action="store_true", default=FALSE,
  6. # help="Show this help message and exit")
  7. option_list <- list(
  8. make_option(c("-p", "--prefix"), type="character",default="./",
  9. help="The output files prefix [default ./]"),
  10. make_option(c("-i", "--input"),type="character", default=NULL,
  11. help="The directory of input EPIC files. required!")
  12. )
  13. # get command line options, if help option encountered print help and exit,
  14. # otherwise if options not found on command line then set defaults,
  15. opt <- parse_args(OptionParser(option_list=option_list))
  16. # load libraries
  17. library("minfi")
  18. library("IlluminaHumanMethylationEPICmanifest")
  19. library("IlluminaHumanMethylationEPICanno.ilm10b4.hg19")
  20. ## 1. data import
  21. targets <- read.metharray.sheet(opt$input)
  22. rgSet <- read.metharray.exp(targets = targets)
  23. targets$ID <- paste(targets$Sample_Group,targets$Sample_Name,sep=".")
  24. sampleNames(rgSet) <- targets$ID
  25. #phenoData <- pData(rgSet)
  26. message("data import:finished")
  27. ## 2. Quality control
  28. # 2.1 qc report by minfi
  29. # qcReport(rgSet, sampNames=targets$ID, sampGroups=targets$Sample_Group,pdf="qcReport.pdf")
  30. # 2.2 data filtering
  31. # get detected p value
  32. # a. remove samples with average p value less than 0.05
  33. detP <- detectionP(rgSet)
  34. message("calculate p value:finished")
  35. # raw
  36. mSetRaw <- preprocessRaw(rgSet)
  37. raw_mVals <- getM(mSetRaw)
  38. raw_bVals <- getBeta(mSetRaw)
  39. message("m value and beta value output:finished")
  40. # normalization
  41. MSet.noob <- preprocessNoob(rgSet)
  42. RSet <- ratioConvert(MSet.noob, what = "both", keepCN = TRUE)
  43. bVals <- getBeta(RSet)
  44. mVals <- getM(RSet)
  45. # write output
  46. message("saving R data")
  47. rdata_filename = paste(opt$prefix, '.RData',sep="")
  48. save(rgSet, targets, detP, file = rdata_filename)
  49. message("writing raw table")
  50. m_raw_filename = paste(opt$prefix,'.raw.mVal.txt',sep="")
  51. write.table(raw_mVals,m_raw_filename,col.names = T,row.names = T,sep="\t",quote=F)
  52. b_raw_filename = paste(opt$prefix,'.raw.bVal.txt',sep="")
  53. write.table(raw_bVals,b_raw_filename,col.names = T,row.names = T,sep="\t",quote=F)
  54. message("writing filtered table")
  55. m_filename = paste(opt$prefix,'.ssnoob.mVal.txt',sep="")
  56. write.table(mVals,m_filename,col.names = T,row.names = T,sep="\t",quote=F)
  57. b_filename = paste(opt$prefix,'.ssnoob.bVal.txt',sep="")
  58. write.table(bVals,b_filename,col.names = T,row.names = T,sep="\t",quote=F)