annotate annotatePeak/.Rhistory @ 0:5d70366e7c6f draft

Uploaded
author testtool
date Mon, 06 Feb 2017 10:53:57 -0500
parents
children 92387cb81962
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
1 TAB <- read.csv("input.csv")
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
2 mysamples <- lapply(TAB$ID,function(x)getGEO(x))
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
3 input <- function(TAB) { if(is(TAB, "csv")){
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
4 TAB <- read.csv("input.csv")}
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
5 else{
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
6 print("error in data file")
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
7 }}
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
8 input()
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
9 TAB <- read.csv("input.csv")}
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
10 TAB <- read.csv("input.csv")
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
11 test_func <- function(
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
12 clusterSize=2,
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
13 cutoff=0.2,
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
14 platform_id='HM450',
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
15 genome_id='hg19')
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
16 {
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
17 args = commandArgs(trailingOnly=TRUE)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
18 methyl_file = args[1]
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
19 ChiPseq_file = args[2]
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
20 output_file = args[3]
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
21 options(warn=-1)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
22 TAB=read.csv(methyl_file)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
23 ChiPseq=import(ChiPseq_file)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
24 if(is.null(TAB)){
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
25 stop("Must specify input files")
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
26 }else{
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
27 mysamples <- lapply(TAB$ID,function(x) getGEO(x))
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
28 }
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
29 wrappedFunction <- function(test_func)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
30 s0 <- lapply(mysamples,Table)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
31 id_ref<-lapply(s0,function(x)x$ID_REF)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
32 if(length(unique(id_ref)) != 1) {
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
33 stop("Error different ID_REF for samples")
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
34 } else if (is.null((unlist(unique(id_ref))))) {
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
35 stop("NO GSM data avaliable")
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
36 } else {
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
37 values<-do.call("cbind",lapply(s0,function(x)x$VALUE))
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
38 colnames(values)=TAB$ID
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
39 rownames(values)=id_ref[[1]]
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
40 cg <- rownames(values)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
41 probe <- c(cg)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
42 hm450.hg19 <- getPlatform(platform=platform_id, genome=genome_id)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
43 probe.info <- hm450.hg19[probe]
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
44 f <- data.table(probe=names(probe.info),CHR=as.data.frame(probe.info@seqnames)$value,
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
45 BP=as.numeric(probe.info@elementMetadata$probeStart))
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
46 designMatrix <- model.matrix(~ TAB$Phenotype)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
47 DMR <- bumphunter(values, design = designMatrix,
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
48 pos=f$BP,cutoff=cutoff,chr=f$CHR)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
49 MAT <- DMR$table[which(DMR$table$L>=clusterSize),]
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
50 METH <- GRanges(seqnames=MAT$chr,
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
51 ranges=IRanges
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
52 (start=MAT$start,
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
53 end=MAT$end),
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
54 value_pos=MAT$value)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
55 peaks<-findOverlapsOfPeaks(probe.info,ChiPseq,maxgap=5000)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
56 p <- peaks$peaklist$`probe.info///ChiPseq`
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
57 peakAnno <- annotatePeak(p, file=peakAnno)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
58 output_file <- annotatePeak(p, file=output_file)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
59 }}
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
60 test_func <- function(
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
61 clusterSize=2,
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
62 cutoff=0.2,
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
63 platform_id='HM450',
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
64 genome_id='hg19')
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
65 {
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
66 args = commandArgs(trailingOnly=TRUE)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
67 methyl_file = args[1]
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
68 ChiPseq_file = args[2]
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
69 output_file = args[3]
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
70 options(warn=-1)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
71 TAB=read.csv(methyl_file)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
72 ChiPseq=import(ChiPseq_file)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
73 if(is.null(TAB)){
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
74 stop("Must specify input files")
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
75 }else{
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
76 mysamples <- lapply(TAB$ID,function(x) getGEO(x))
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
77 }
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
78 wrappedFunction <- function(test_func)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
79 s0 <- lapply(mysamples,Table)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
80 id_ref<-lapply(s0,function(x)x$ID_REF)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
81 if(length(unique(id_ref)) != 1) {
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
82 stop("Error different ID_REF for samples")
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
83 } else if (is.null((unlist(unique(id_ref))))) {
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
84 stop("NO GSM data avaliable")
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
85 } else {
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
86 values<-do.call("cbind",lapply(s0,function(x)x$VALUE))
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
87 colnames(values)=TAB$ID
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
88 rownames(values)=id_ref[[1]]
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
89 cg <- rownames(values)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
90 probe <- c(cg)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
91 hm450.hg19 <- getPlatform(platform=platform_id, genome=genome_id)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
92 probe.info <- hm450.hg19[probe]
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
93 f <- data.table(probe=names(probe.info),CHR=as.data.frame(probe.info@seqnames)$value,
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
94 BP=as.numeric(probe.info@elementMetadata$probeStart))
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
95 designMatrix <- model.matrix(~ TAB$Phenotype)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
96 DMR <- bumphunter(values, design = designMatrix,
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
97 pos=f$BP,cutoff=cutoff,chr=f$CHR)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
98 MAT <- DMR$table[which(DMR$table$L>=clusterSize),]
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
99 METH <- GRanges(seqnames=MAT$chr,
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
100 ranges=IRanges
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
101 (start=MAT$start,
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
102 end=MAT$end),
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
103 value_pos=MAT$value)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
104 peaks<-findOverlapsOfPeaks(probe.info,ChiPseq,maxgap=5000)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
105 p <- peaks$peaklist$`probe.info///ChiPseq`
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
106 peakAnno <- annotatePeak(p, file=peakAnno)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
107 output_file <- annotatePeak(p, file=output_file)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
108 }}
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
109 )
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
110 a<-c(3.2,3.7,3.6,3.9,3.7,3.5,3.8,4,3.5)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
111 a<-c(3.2,3.7,3.6,3.9,3.7,3.5,3.8,4,3.5)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
112 b<-c(2.9,2.7,2,6,2.7,2.5,2.5,3.2,3.1,3.3,2.8)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
113 a-b
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
114 file<-read.csv("~/Documents/SS2.csv")
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
115 head(file)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
116 e
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
117 file
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
118 t.test(file$WBT,file$WBA)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
119 t.test(file$WBT,file$WBA,paired=TRUE)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
120 mysamples <- lapply(TAB$ID,function(x) TablegetGEO(x))
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
121 require("GEOquery",quietly = TRUE)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
122 require("BiocGenerics",quietly = TRUE)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
123 args <- commandArgs(trailingOnly = TRUE)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
124 csv_file = args[1]
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
125 #csv_file <- ("test-data/input.csv")
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
126 TAB=read.csv(csv_file)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
127 if(is.null(TAB)){
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
128 stop("Must specify input files")
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
129 }else{
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
130 mysamples <- lapply(TAB$ID,function(x) TablegetGEO(x))
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
131 }
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
132 csv_file <- ("test-data/input.csv")
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
133 TAB=read.csv(csv_file)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
134 figure = args[3]
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
135 upsetplot(anno, vennpie=TRUE)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
136 dev.copy(png,'fig.png')
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
137 png(upsetplot(anno, vennpie=TRUE))
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
138 upsetplot(anno, vennpie=TRUE)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
139 require("ChIPseeker",quietly = TRUE)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
140 require("ChIPpeakAnno",quietly = TRUE)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
141 upsetplot(anno, vennpie=TRUE)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
142 figure<-('test-data/figure.png')
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
143 anno <- annotatePeak(peaks,level="gene", annoDb="org.Hs.eg.db" )
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
144 input=read.table(input,header = FALSE, sep="\t",stringsAsFactors=FALSE, quote="")
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
145 input<-("test-data/bump.bed")
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
146 input=read.table(input,header = FALSE, sep="\t",stringsAsFactors=FALSE, quote="")
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
147 peaks <- GRanges(seqnames=input[,1],ranges=IRanges
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
148 (start=input[,2],end= input[,3]))
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
149 anno <- annotatePeak(peaks,level="gene", annoDb="org.Hs.eg.db" )
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
150 peakAnno_genes <- as.data.frame(anno)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
151 upsetplot(anno, vennpie=TRUE)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
152 figure<-('test-data/figure.png')
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
153 png(figure)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
154 dev.off()
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
155 png(figure)
5d70366e7c6f Uploaded
testtool
parents:
diff changeset
156 upsetplot(anno, vennpie=TRUE)