Mercurial > repos > testtool > geo_data
view getGRsetFromGEO/.Rhistory @ 11:983651783ce9 draft
Uploaded
author | testtool |
---|---|
date | Wed, 22 Mar 2017 09:01:06 -0400 |
parents | |
children |
line wrap: on
line source
TAB <- read.csv("input.csv") mysamples <- lapply(TAB$ID,function(x)getGEO(x)) input <- function(TAB) { if(is(TAB, "csv")){ TAB <- read.csv("input.csv")} else{ print("error in data file") }} input() TAB <- read.csv("input.csv")} TAB <- read.csv("input.csv") test_func <- function( clusterSize=2, cutoff=0.2, platform_id='HM450', genome_id='hg19') { args = commandArgs(trailingOnly=TRUE) methyl_file = args[1] ChiPseq_file = args[2] output_file = args[3] options(warn=-1) TAB=read.csv(methyl_file) ChiPseq=import(ChiPseq_file) if(is.null(TAB)){ stop("Must specify input files") }else{ mysamples <- lapply(TAB$ID,function(x) getGEO(x)) } wrappedFunction <- function(test_func) s0 <- lapply(mysamples,Table) id_ref<-lapply(s0,function(x)x$ID_REF) if(length(unique(id_ref)) != 1) { stop("Error different ID_REF for samples") } else if (is.null((unlist(unique(id_ref))))) { stop("NO GSM data avaliable") } else { values<-do.call("cbind",lapply(s0,function(x)x$VALUE)) colnames(values)=TAB$ID rownames(values)=id_ref[[1]] cg <- rownames(values) probe <- c(cg) hm450.hg19 <- getPlatform(platform=platform_id, genome=genome_id) probe.info <- hm450.hg19[probe] f <- data.table(probe=names(probe.info),CHR=as.data.frame(probe.info@seqnames)$value, BP=as.numeric(probe.info@elementMetadata$probeStart)) designMatrix <- model.matrix(~ TAB$Phenotype) DMR <- bumphunter(values, design = designMatrix, pos=f$BP,cutoff=cutoff,chr=f$CHR) MAT <- DMR$table[which(DMR$table$L>=clusterSize),] METH <- GRanges(seqnames=MAT$chr, ranges=IRanges (start=MAT$start, end=MAT$end), value_pos=MAT$value) peaks<-findOverlapsOfPeaks(probe.info,ChiPseq,maxgap=5000) p <- peaks$peaklist$`probe.info///ChiPseq` peakAnno <- annotatePeak(p, file=peakAnno) output_file <- annotatePeak(p, file=output_file) }} test_func <- function( clusterSize=2, cutoff=0.2, platform_id='HM450', genome_id='hg19') { args = commandArgs(trailingOnly=TRUE) methyl_file = args[1] ChiPseq_file = args[2] output_file = args[3] options(warn=-1) TAB=read.csv(methyl_file) ChiPseq=import(ChiPseq_file) if(is.null(TAB)){ stop("Must specify input files") }else{ mysamples <- lapply(TAB$ID,function(x) getGEO(x)) } wrappedFunction <- function(test_func) s0 <- lapply(mysamples,Table) id_ref<-lapply(s0,function(x)x$ID_REF) if(length(unique(id_ref)) != 1) { stop("Error different ID_REF for samples") } else if (is.null((unlist(unique(id_ref))))) { stop("NO GSM data avaliable") } else { values<-do.call("cbind",lapply(s0,function(x)x$VALUE)) colnames(values)=TAB$ID rownames(values)=id_ref[[1]] cg <- rownames(values) probe <- c(cg) hm450.hg19 <- getPlatform(platform=platform_id, genome=genome_id) probe.info <- hm450.hg19[probe] f <- data.table(probe=names(probe.info),CHR=as.data.frame(probe.info@seqnames)$value, BP=as.numeric(probe.info@elementMetadata$probeStart)) designMatrix <- model.matrix(~ TAB$Phenotype) DMR <- bumphunter(values, design = designMatrix, pos=f$BP,cutoff=cutoff,chr=f$CHR) MAT <- DMR$table[which(DMR$table$L>=clusterSize),] METH <- GRanges(seqnames=MAT$chr, ranges=IRanges (start=MAT$start, end=MAT$end), value_pos=MAT$value) peaks<-findOverlapsOfPeaks(probe.info,ChiPseq,maxgap=5000) p <- peaks$peaklist$`probe.info///ChiPseq` peakAnno <- annotatePeak(p, file=peakAnno) output_file <- annotatePeak(p, file=output_file) }} ) a<-c(3.2,3.7,3.6,3.9,3.7,3.5,3.8,4,3.5) a<-c(3.2,3.7,3.6,3.9,3.7,3.5,3.8,4,3.5) b<-c(2.9,2.7,2,6,2.7,2.5,2.5,3.2,3.1,3.3,2.8) a-b file<-read.csv("~/Documents/SS2.csv") head(file) e file t.test(file$WBT,file$WBA) t.test(file$WBT,file$WBA,paired=TRUE) mysamples <- lapply(TAB$ID,function(x) TablegetGEO(x)) require("GEOquery",quietly = TRUE) require("BiocGenerics",quietly = TRUE) args <- commandArgs(trailingOnly = TRUE) csv_file = args[1] #csv_file <- ("test-data/input.csv") TAB=read.csv(csv_file) if(is.null(TAB)){ stop("Must specify input files") }else{ mysamples <- lapply(TAB$ID,function(x) TablegetGEO(x)) } csv_file <- ("test-data/input.csv") TAB=read.csv(csv_file) TAB csv_file<-("test-data/input.csv") TAB = read.csv(csv_file) TAB ??bumphunter