diff findDMR/findDMR.R @ 12:083895dbc289 draft

Uploaded

--- a/findDMR/findDMR.R	Thu May 11 11:26:02 2017 -0400
+++ b/findDMR/findDMR.R	Mon Jun 12 11:05:04 2017 -0400
@@ -1,42 +1,36 @@
-require("BiocGenerics", quietly = TRUE)
-require("data.table", quietly = TRUE)
-require("bumphunter", quietly = TRUE)
+require("minfi", quietly = TRUE)
+
+options(warn = -1)
+options("download.file.method"="wget")
 
 args <- commandArgs(trailingOnly = TRUE)
-GSMTable = args[1]
-platform = args[2]
-Data_Table = args[3]
-cutoff = as.numeric(args[4])
-clusterSize = as.numeric(args[5])
-DMR = args[6]
+input1 = args[1]
+input2 = args[2]
+cutoff = as.numeric(args[3])
+B = as.numeric(args[4]) 
+pickCutoffQ = as.numeric(args[5])
+output1 = args[6]
+output2 = args[7]
 
-TAB = fread(GSMTable)
-
-IlmnInfo = fread(platform)
+GRset <- get(load(input1))
 
-gmSet = fread(Data_Table)
+pheno <- fread(input2)
 
-# bumphunter Run with processed data
-designMatrix <- model.matrix( ~ TAB$Phenotype)
+designMatrix <- model.matrix(~ pheno$Phenotype)
 
-bumps <- bumphunter(
-  as.matrix(gmSet),
-  design = designMatrix,
-  pos = IlmnInfo$BP,
-  cutoff = cutoff,
-  chr = IlmnInfo$CHR
-)
+dmrs <- bumphunter(GRset, design = designMatrix,
+                   cutoff =cutoff, B=B, type="Beta",  pickCutoff=TRUE,
+                   pickCutoffQ=pickCutoffQ)
+
+DMRTable <- dmrs$table
+
+write.table(DMRTable, output1)
 
-# choose DMR's of a certain length threshold
-DMRTable <- bumps$table[which(bumps$table$L >= clusterSize), ]
-DMRInfo <- data.table(DMRTable$chr, DMRTable$start, DMRTable$end)
-
+sign=sign(DMRTable$value)
 
-write.table(
-  DMRInfo,
-  DMR,
-  quote = F,
-  sep = "\t",
-  row.names = F,
-  col.names = F
-)
+sign[sign==-1]="-"
+sign[sign==1]="+"
+
+dmr_track<-cbind(as.character(DMRTable$chr),DMRTable$start,DMRTable$end,sign)
+
+write.table(dmr_track,output2,quote = FALSE, sep = "\t\t\t\t",row.names = FALSE, col.names = FALSE,append=TRUE)
\ No newline at end of file