view annoPeakFigure/.Rhistory @ 5:1de0d1429a3c draft default tip

Uploaded

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
require("ChIPseeker", quietly = TRUE)
require("ChIPpeakAnno", quietly = TRUE)
options(warn = -1)
#args <- commandArgs(trailingOnly = TRUE)
#DMR = args[1]
#annoPeakFigure = args[2]
setwd('/Users/katarzynamurat/Desktop/galaxy/test-data')
DMR <- ("DMR.bed")
DMRInfo = read.table(
DMR,
header = FALSE,
sep = "\t",
stringsAsFactors = FALSE,
quote = ""
)
DMRPeaks <- GRanges(seqnames = DMRInfo[, 1],
ranges = IRanges
(start = DMRInfo[, 2], end = DMRInfo[, 3]))
annotatePeak <- (annotatePeak(DMRPeaks, level = "gene", annoDb = "org.Hs.eg.db"))
annoPeakFigure <- ('annoPeakFigure.png')
png(file = annoPeakFigure, width = 1500,height = 700, units = "px", bg = "white")
genes= lapply(annotatePeak, function(i) as.data.frame(i)$geneId)
vennplot(genes)
dev.off()
genes= lapply(as.list(annotatePeak), function(i) as.data.frame(i)$geneId)
vennplot(genes)
peakAnnoList <- lapply(DMRPeaks, annotatePeak)
peakAnnoList
annotatePeak <- (annotatePeak(DMRPeaks, level = "gene", annoDb = "org.Hs.eg.db"))
peakAnnoList <- lapply(DMRPeaks, annotatePeak)
peakAnnoList
annoPeakFigure <- ('annoPeakFigure.png')
png(file = annoPeakFigure, width = 1500,height = 700, units = "px", bg = "white")
plotAnnoBar(peakAnno)
dev.off()
annotatePeak <- (annotatePeak(DMRPeaks, level = "gene", annoDb = "org.Hs.eg.db"))
annoPeakFigure <- ('annoPeakFigure.png')
png(file = annoPeakFigure, width = 1500,height = 700, units = "px", bg = "white")
plotAnnoBar(annotatePeak)
dev.off()
png(file = annoPeakFigure, width = 800,height = 400, units = "px", bg = "white")
plotAnnoBar(annotatePeak)
dev.off()
png(file = annoPeakFigure, width = 800,height = 400, units = "px", bg = "white", pointsize = 12)
plotAnnoBar(annotatePeak)
dev.off()
annoPeakFigure <- ('annoPeakFigure.png')
png(file = annoPeakFigure, width = 800,height = 400, units = "px", bg = "white", pointsize = 12)
plotAnnoBar(annotatePeak)
dev.off()
options(warn = -1)
args <- commandArgs(trailingOnly = TRUE)
DMR = args[1]
annoPeakFigure = args[2]
DMR <- ("test-data/DMR.bed")
DMRInfo = read.table(
DMR,
header = FALSE,
sep = "\t",
stringsAsFactors = FALSE,
quote = ""
)
DMRPeaks <- GRanges(seqnames = DMRInfo[, 1],
ranges = IRanges
(start = DMRInfo[, 2], end = DMRInfo[, 3]))
annotatePeak <- (annotatePeak(DMRPeaks, level = "gene", annoDb = "org.Hs.eg.db"))
annoPeakFigure <- ('test-data/annoPeakFigure.png')
png(file = annoPeakFigure,
width = 1200,
height = 600)
plotAnnoBar(annotatePeak)
upsetplot(annotatePeak, vennpie = TRUE)
dev.off()
annoPeakFigure1 <- ('test-data/annoPeakUpset.png')
annoPeakFigure2 <- ('test-data/annoPeakAnnoBar.png')
png(file = annoPeakFigure1,
width = 1200,
height = 600)
plotAnnoBar(annotatePeak)
dev.off()
png(file = annoPeakFigure2,
width = 1200,
height = 600)
upsetplot(annotatePeak, vennpie = TRUE)
dev.off()
getwd()
setwd("/Users/katarzynamurat/Desktop/annoPeakFigure")
annoPeakFigure1 <- ('test-data/annoPeakUpset.png')
annoPeakFigure2 <- ('test-data/annoPeakAnnoBar.png')
png(file = annoPeakFigure1,
width = 1200,
height = 600)
plotAnnoBar(annotatePeak)
dev.off()
png(file = annoPeakFigure2,
width = 1200,
height = 600)
upsetplot(annotatePeak, vennpie = TRUE)
dev.off()
require("ChIPseeker", quietly = TRUE)
require("ChIPpeakAnno", quietly = TRUE)
options(warn = -1)
DMR <- ("test-data/DMR.bed")
DMRInfo = read.table(
DMR,
header = FALSE,
sep = "\t",
stringsAsFactors = FALSE,
quote = ""
)
DMRPeaks <- GRanges(seqnames = DMRInfo[, 1],
ranges = IRanges
(start = DMRInfo[, 2], end = DMRInfo[, 3]))
annotatePeak <- (annotatePeak(DMRPeaks, level = "gene", annoDb = "org.Hs.eg.db"))
annoPeakFigure1 <- ('test-data/annoPeakUpset.png')
annoPeakFigure2 <- ('test-data/annoPeakAnnoBar.png')
png(file = annoPeakFigure1,
width = 1200,
height = 600)
plotAnnoBar(annotatePeak)
dev.off()
png(file = annoPeakFigure2,
width = 1200,
height = 600)
upsetplot(annotatePeak, vennpie = TRUE)
dev.off()
getwd()
??GRanges