|
0
|
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 figure = args[3]
|
|
|
135 upsetplot(anno, vennpie=TRUE)
|
|
|
136 dev.copy(png,'fig.png')
|
|
|
137 png(upsetplot(anno, vennpie=TRUE))
|
|
|
138 upsetplot(anno, vennpie=TRUE)
|
|
|
139 require("ChIPseeker",quietly = TRUE)
|
|
|
140 require("ChIPpeakAnno",quietly = TRUE)
|
|
|
141 upsetplot(anno, vennpie=TRUE)
|
|
|
142 figure<-('test-data/figure.png')
|
|
|
143 anno <- annotatePeak(peaks,level="gene", annoDb="org.Hs.eg.db" )
|
|
|
144 input=read.table(input,header = FALSE, sep="\t",stringsAsFactors=FALSE, quote="")
|
|
|
145 input<-("test-data/bump.bed")
|
|
|
146 input=read.table(input,header = FALSE, sep="\t",stringsAsFactors=FALSE, quote="")
|
|
|
147 peaks <- GRanges(seqnames=input[,1],ranges=IRanges
|
|
|
148 (start=input[,2],end= input[,3]))
|
|
|
149 anno <- annotatePeak(peaks,level="gene", annoDb="org.Hs.eg.db" )
|
|
|
150 peakAnno_genes <- as.data.frame(anno)
|
|
|
151 upsetplot(anno, vennpie=TRUE)
|
|
|
152 figure<-('test-data/figure.png')
|
|
|
153 png(figure)
|
|
|
154 dev.off()
|
|
|
155 png(figure)
|
|
|
156 upsetplot(anno, vennpie=TRUE)
|
|
1
|
157 org.Hs.eg.db
|
|
|
158 ??org.Hs.eg.db
|
|
|
159 require("ChIPseeker", quietly = TRUE)
|
|
|
160 require("ChIPpeakAnno", quietly = TRUE)
|
|
|
161 DMRInfo = read.table(
|
|
|
162 DMR,
|
|
|
163 header = FALSE,
|
|
|
164 sep = "\t",
|
|
|
165 stringsAsFactors = FALSE,
|
|
|
166 quote = ""
|
|
|
167 )
|
|
|
168 DMRPeaks <- GRanges(seqnames = DMRInfo[, 1],
|
|
|
169 ranges = IRanges
|
|
|
170 (start = DMRInfo[, 2], end = DMRInfo[, 3]))
|
|
|
171 annotatePeak <-
|
|
|
172 as.data.frame(annotatePeak(DMRPeaks, level = "gene", annoDb = "org.Hs.eg.db"))
|
|
|
173 ??org.Hs.eg.db
|
