|
11
|
1 TAB <- read.csv("input.csv")
|
|
|
2 mysamples <- lapply(TAB$ID,function(x)getGEO(x))
|
|
|
3 input <- function(TAB) { if(is(TAB, "csv")){
|
|
|
4 TAB <- read.csv("input.csv")}
|
|
|
5 else{
|
|
|
6 print("error in data file")
|
|
|
7 }}
|
|
|
8 input()
|
|
|
9 TAB <- read.csv("input.csv")}
|
|
|
10 TAB <- read.csv("input.csv")
|
|
|
11 test_func <- function(
|
|
|
12 clusterSize=2,
|
|
|
13 cutoff=0.2,
|
|
|
14 platform_id='HM450',
|
|
|
15 genome_id='hg19')
|
|
|
16 {
|
|
|
17 args = commandArgs(trailingOnly=TRUE)
|
|
|
18 methyl_file = args[1]
|
|
|
19 ChiPseq_file = args[2]
|
|
|
20 output_file = args[3]
|
|
|
21 options(warn=-1)
|
|
|
22 TAB=read.csv(methyl_file)
|
|
|
23 ChiPseq=import(ChiPseq_file)
|
|
|
24 if(is.null(TAB)){
|
|
|
25 stop("Must specify input files")
|
|
|
26 }else{
|
|
|
27 mysamples <- lapply(TAB$ID,function(x) getGEO(x))
|
|
|
28 }
|
|
|
29 wrappedFunction <- function(test_func)
|
|
|
30 s0 <- lapply(mysamples,Table)
|
|
|
31 id_ref<-lapply(s0,function(x)x$ID_REF)
|
|
|
32 if(length(unique(id_ref)) != 1) {
|
|
|
33 stop("Error different ID_REF for samples")
|
|
|
34 } else if (is.null((unlist(unique(id_ref))))) {
|
|
|
35 stop("NO GSM data avaliable")
|
|
|
36 } else {
|
|
|
37 values<-do.call("cbind",lapply(s0,function(x)x$VALUE))
|
|
|
38 colnames(values)=TAB$ID
|
|
|
39 rownames(values)=id_ref[[1]]
|
|
|
40 cg <- rownames(values)
|
|
|
41 probe <- c(cg)
|
|
|
42 hm450.hg19 <- getPlatform(platform=platform_id, genome=genome_id)
|
|
|
43 probe.info <- hm450.hg19[probe]
|
|
|
44 f <- data.table(probe=names(probe.info),CHR=as.data.frame(probe.info@seqnames)$value,
|
|
|
45 BP=as.numeric(probe.info@elementMetadata$probeStart))
|
|
|
46 designMatrix <- model.matrix(~ TAB$Phenotype)
|
|
|
47 DMR <- bumphunter(values, design = designMatrix,
|
|
|
48 pos=f$BP,cutoff=cutoff,chr=f$CHR)
|
|
|
49 MAT <- DMR$table[which(DMR$table$L>=clusterSize),]
|
|
|
50 METH <- GRanges(seqnames=MAT$chr,
|
|
|
51 ranges=IRanges
|
|
|
52 (start=MAT$start,
|
|
|
53 end=MAT$end),
|
|
|
54 value_pos=MAT$value)
|
|
|
55 peaks<-findOverlapsOfPeaks(probe.info,ChiPseq,maxgap=5000)
|
|
|
56 p <- peaks$peaklist$`probe.info///ChiPseq`
|
|
|
57 peakAnno <- annotatePeak(p, file=peakAnno)
|
|
|
58 output_file <- annotatePeak(p, file=output_file)
|
|
|
59 }}
|
|
|
60 test_func <- function(
|
|
|
61 clusterSize=2,
|
|
|
62 cutoff=0.2,
|
|
|
63 platform_id='HM450',
|
|
|
64 genome_id='hg19')
|
|
|
65 {
|
|
|
66 args = commandArgs(trailingOnly=TRUE)
|
|
|
67 methyl_file = args[1]
|
|
|
68 ChiPseq_file = args[2]
|
|
|
69 output_file = args[3]
|
|
|
70 options(warn=-1)
|
|
|
71 TAB=read.csv(methyl_file)
|
|
|
72 ChiPseq=import(ChiPseq_file)
|
|
|
73 if(is.null(TAB)){
|
|
|
74 stop("Must specify input files")
|
|
|
75 }else{
|
|
|
76 mysamples <- lapply(TAB$ID,function(x) getGEO(x))
|
|
|
77 }
|
|
|
78 wrappedFunction <- function(test_func)
|
|
|
79 s0 <- lapply(mysamples,Table)
|
|
|
80 id_ref<-lapply(s0,function(x)x$ID_REF)
|
|
|
81 if(length(unique(id_ref)) != 1) {
|
|
|
82 stop("Error different ID_REF for samples")
|
|
|
83 } else if (is.null((unlist(unique(id_ref))))) {
|
|
|
84 stop("NO GSM data avaliable")
|
|
|
85 } else {
|
|
|
86 values<-do.call("cbind",lapply(s0,function(x)x$VALUE))
|
|
|
87 colnames(values)=TAB$ID
|
|
|
88 rownames(values)=id_ref[[1]]
|
|
|
89 cg <- rownames(values)
|
|
|
90 probe <- c(cg)
|
|
|
91 hm450.hg19 <- getPlatform(platform=platform_id, genome=genome_id)
|
|
|
92 probe.info <- hm450.hg19[probe]
|
|
|
93 f <- data.table(probe=names(probe.info),CHR=as.data.frame(probe.info@seqnames)$value,
|
|
|
94 BP=as.numeric(probe.info@elementMetadata$probeStart))
|
|
|
95 designMatrix <- model.matrix(~ TAB$Phenotype)
|
|
|
96 DMR <- bumphunter(values, design = designMatrix,
|
|
|
97 pos=f$BP,cutoff=cutoff,chr=f$CHR)
|
|
|
98 MAT <- DMR$table[which(DMR$table$L>=clusterSize),]
|
|
|
99 METH <- GRanges(seqnames=MAT$chr,
|
|
|
100 ranges=IRanges
|
|
|
101 (start=MAT$start,
|
|
|
102 end=MAT$end),
|
|
|
103 value_pos=MAT$value)
|
|
|
104 peaks<-findOverlapsOfPeaks(probe.info,ChiPseq,maxgap=5000)
|
|
|
105 p <- peaks$peaklist$`probe.info///ChiPseq`
|
|
|
106 peakAnno <- annotatePeak(p, file=peakAnno)
|
|
|
107 output_file <- annotatePeak(p, file=output_file)
|
|
|
108 }}
|
|
|
109 )
|
|
|
110 a<-c(3.2,3.7,3.6,3.9,3.7,3.5,3.8,4,3.5)
|
|
|
111 a<-c(3.2,3.7,3.6,3.9,3.7,3.5,3.8,4,3.5)
|
|
|
112 b<-c(2.9,2.7,2,6,2.7,2.5,2.5,3.2,3.1,3.3,2.8)
|
|
|
113 a-b
|
|
|
114 file<-read.csv("~/Documents/SS2.csv")
|
|
|
115 head(file)
|
|
|
116 e
|
|
|
117 file
|
|
|
118 t.test(file$WBT,file$WBA)
|
|
|
119 t.test(file$WBT,file$WBA,paired=TRUE)
|
|
|
120 mysamples <- lapply(TAB$ID,function(x) TablegetGEO(x))
|
|
|
121 require("GEOquery",quietly = TRUE)
|
|
|
122 require("BiocGenerics",quietly = TRUE)
|
|
|
123 args <- commandArgs(trailingOnly = TRUE)
|
|
|
124 csv_file = args[1]
|
|
|
125 #csv_file <- ("test-data/input.csv")
|
|
|
126 TAB=read.csv(csv_file)
|
|
|
127 if(is.null(TAB)){
|
|
|
128 stop("Must specify input files")
|
|
|
129 }else{
|
|
|
130 mysamples <- lapply(TAB$ID,function(x) TablegetGEO(x))
|
|
|
131 }
|
|
|
132 csv_file <- ("test-data/input.csv")
|
|
|
133 TAB=read.csv(csv_file)
|
|
|
134 TAB
|
|
|
135 csv_file<-("test-data/input.csv")
|
|
|
136 TAB = read.csv(csv_file)
|
|
|
137 TAB
|
|
|
138 ??bumphunter
|
