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

Uploaded
author testtool
date Mon, 12 Jun 2017 11:05:04 -0400
parents 24ac6f93cc3e
children
comparison
equal deleted inserted replaced
11:24ac6f93cc3e 12:083895dbc289
1 require("BiocGenerics", quietly = TRUE) 1 require("minfi", quietly = TRUE)
2 require("data.table", quietly = TRUE) 2
3 require("bumphunter", quietly = TRUE) 3 options(warn = -1)
4 options("download.file.method"="wget")
4 5
5 args <- commandArgs(trailingOnly = TRUE) 6 args <- commandArgs(trailingOnly = TRUE)
6 GSMTable = args[1] 7 input1 = args[1]
7 platform = args[2] 8 input2 = args[2]
8 Data_Table = args[3] 9 cutoff = as.numeric(args[3])
9 cutoff = as.numeric(args[4]) 10 B = as.numeric(args[4])
10 clusterSize = as.numeric(args[5]) 11 pickCutoffQ = as.numeric(args[5])
11 DMR = args[6] 12 output1 = args[6]
13 output2 = args[7]
12 14
13 TAB = fread(GSMTable) 15 GRset <- get(load(input1))
14 16
15 IlmnInfo = fread(platform) 17 pheno <- fread(input2)
16 18
17 gmSet = fread(Data_Table) 19 designMatrix <- model.matrix(~ pheno$Phenotype)
18 20
19 # bumphunter Run with processed data 21 dmrs <- bumphunter(GRset, design = designMatrix,
20 designMatrix <- model.matrix( ~ TAB$Phenotype) 22 cutoff =cutoff, B=B, type="Beta", pickCutoff=TRUE,
23 pickCutoffQ=pickCutoffQ)
21 24
22 bumps <- bumphunter( 25 DMRTable <- dmrs$table
23 as.matrix(gmSet),
24 design = designMatrix,
25 pos = IlmnInfo$BP,
26 cutoff = cutoff,
27 chr = IlmnInfo$CHR
28 )
29 26
30 # choose DMR's of a certain length threshold 27 write.table(DMRTable, output1)
31 DMRTable <- bumps$table[which(bumps$table$L >= clusterSize), ]
32 DMRInfo <- data.table(DMRTable$chr, DMRTable$start, DMRTable$end)
33 28
29 sign=sign(DMRTable$value)
34 30
35 write.table( 31 sign[sign==-1]="-"
36 DMRInfo, 32 sign[sign==1]="+"
37 DMR, 33
38 quote = F, 34 dmr_track<-cbind(as.character(DMRTable$chr),DMRTable$start,DMRTable$end,sign)
39 sep = "\t", 35
40 row.names = F, 36 write.table(dmr_track,output2,quote = FALSE, sep = "\t\t\t\t",row.names = FALSE, col.names = FALSE,append=TRUE)
41 col.names = F
42 )