Mercurial > repos > testtool > annotate_peak

comparison annotatePeak/.Rhistory @ 0:5d70366e7c6f draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author testtool
date Mon, 06 Feb 2017 10:53:57 -0500
parents
children 92387cb81962
comparison
equal deleted inserted replaced
-1:000000000000 0:5d70366e7c6f
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)