Mercurial > repos > pjbriggs > ceas

view test-data/ceas_out1.log @ 0:f411ce97a351 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded initial version 1.0.2-2
author pjbriggs
date Tue, 30 Jun 2015 07:08:05 -0400
parents
children
line wrap: on
line source

ceas -- 0.9.9.7 (package version 1.0.2)
INFO  @ Tue, 23 Jun 2015 09:12:22: 
# ARGUMENTS: 
# name = ceas
# gene annotation table = galGal3.refGene
# BED file = ceas_in.bed
# WIG file = None
# extra BED file = None
# ChIP annotation = On
# gene-centered annotation =  On
# average profiling = Off
# dump profiles = Off
# re-annotation for genome background (ChIP region annotation) = False
# promoter sizes (ChIP region annotation) = 1000,2000,3000 bp
# downstream sizes (ChIP region annotation) = 1000,2000,3000 bp
# bidrectional promoter sizes (ChIP region annotation) = 2500,5000 bp
# span size (gene-centered annotation) = 3000 bp 
INFO  @ Tue, 23 Jun 2015 09:12:22: #1 read the gene table... 
INFO  @ Tue, 23 Jun 2015 09:12:22: #2 read the bed file of ChIP regions... 
INFO  @ Tue, 23 Jun 2015 09:12:22: #3 perform gene-centered annotation... 
INFO  @ Tue, 23 Jun 2015 09:12:22: #4 See ceas.xls for gene-centered annotation! 
INFO  @ Tue, 23 Jun 2015 09:12:22: #5 read the pre-computed genome bg annotation... 
INFO  @ Tue, 23 Jun 2015 09:12:22: #6 perform ChIP region annotation... 
INFO  @ Tue, 23 Jun 2015 09:12:22: #7 write a R script of ChIP region annotation... 

R version 3.1.2 (2014-10-31) -- "Pumpkin Helmet"
Copyright (C) 2014 The R Foundation for Statistical Computing
Platform: x86_64-redhat-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> # ARGUMENTS: 
> # name = ceas
> # gene annotation table = galGal3.refGene
> # BED file = ceas_in.bed
> # WIG file = None
> # extra BED file = None
> # ChIP annotation = On
> # gene-centered annotation =  On
> # average profiling = Off
> # dump profiles = Off
> # re-annotation for genome background (ChIP region annotation) = False
> # promoter sizes (ChIP region annotation) = 1000,2000,3000 bp
> # downstream sizes (ChIP region annotation) = 1000,2000,3000 bp
> # bidrectional promoter sizes (ChIP region annotation) = 2500,5000 bp
> # span size (gene-centered annotation) = 3000 bp
> pdf("ceas.pdf",height=11.5,width=8.5)
> 
> # 09:12:22 Tue, 23 Jun 2015
> # 
> # ChIP annotation
> # 
> 
> 
> # 
> # Chromosomal Distribution
> # 
> 
> par(mar=c(4, 4, 5, 3.8),oma=c(4, 2, 4, 2))
> r0<-c(100.0)
> r1<-c(100.0)
> height<-rbind(r0,r1)
> names=c("26")
> mp<-barplot(height=height,names=names,beside=TRUE,horiz=TRUE,col=c("#5FA1C1","#EB9D86"),main="Chromosomal Distribution of ChIP Regions",xlab="Percentage %",ylab="Chromosome",border=FALSE,xlim=c(0.000000,183.333333),cex.names=1)
> text(x=c(100.0),y=mp[1,],label=c("100.0 %"),pos=4,offset=0.2,cex=0.9)
> text(x=c(100.0),y=mp[2,],label=c("100.0 % (<=4.9e-324)"),pos=4,offset=0.2,cex=0.9)
> legend("right",legend=c("Genome","ChIP (p-value)"),col=c("#5FA1C1","#EB9D86"),pch=15,bty="n")
> 
> # 
> # Promoter,Bipromoter,Downstream, Gene and Regions of interest
> # 
> 
> par(mfrow=c(4, 1),mar=c(4, 4, 5, 3.8),oma=c(4, 2, 4, 2))
> r0<-c(1.8532425688606797, 3.616851183410451, 5.322318854623416)
> r1<-c(0.0, 0.0, 0.0)
> height<-rbind(r0,r1)
> names=c("<=1000 bp","<=2000 bp","<=3000 bp")
> mp<-barplot(height=height,names=names,beside=TRUE,horiz=FALSE,col=c("#5FA1C1","#EB9D86"),main="Promoter",ylab="Percentage %",border=FALSE,ylim=c(0.000000,9.757585),cex.names=1)
> text(x=mp[1,],y=c(1.8532425688606797, 3.616851183410451, 5.322318854623416),label=c("1.9 %","3.6 %","5.3 %"),pos=3,offset=0.2)
> text(x=mp[2,],y=c(0.0, 0.0, 0.0),label=c("0.000 %
+ (0.981)","0.000 %
+ (0.964)","0.000 %
+ (0.947)"),pos=3,offset=0.2)
> legend("topleft",legend=c("Genome","ChIP (p-value)"),col=c("#5FA1C1","#EB9D86"),pch=15,bty="n")
> r0<-c(0.03876062889120376, 0.03876062889120376)
> r1<-c(0.0, 0.0)
> height<-rbind(r0,r1)
> names=c("<=2500 bp","<=5000 bp")
> mp<-barplot(height=height,names=names,beside=TRUE,horiz=FALSE,col=c("#5FA1C1","#EB9D86"),main="Bidirectional Promoter",ylab="Percentage %",border=FALSE,ylim=c(0.000000,0.071061),cex.names=1)
> text(x=mp[1,],y=c(0.03876062889120376, 0.03876062889120376),label=c("0.04 %","0.04 %"),pos=3,offset=0.2)
> text(x=mp[2,],y=c(0.0, 0.0),label=c("0.000 %
+ (1.000)","0.000 %
+ (1.000)"),pos=3,offset=0.2)
> legend("topleft",legend=c("Genome","ChIP (p-value)"),col=c("#5FA1C1","#EB9D86"),pch=15,bty="n")
> r0<-c(1.8290171758036773, 3.4690762857627364, 4.980740812519683)
> r1<-c(0.0, 0.0, 0.0)
> height<-rbind(r0,r1)
> names=c("<=1000 bp","<=2000 bp","<=3000 bp")
> mp<-barplot(height=height,names=names,beside=TRUE,horiz=FALSE,col=c("#5FA1C1","#EB9D86"),main="Downstream",ylab="Percentage %",border=FALSE,ylim=c(0.000000,9.131358),cex.names=1)
> text(x=mp[1,],y=c(1.8290171758036773, 3.4690762857627364, 4.980740812519683),label=c("1.8 %","3.5 %","5.0 %"),pos=3,offset=0.2)
> text(x=mp[2,],y=c(0.0, 0.0, 0.0),label=c("0.000 %
+ (0.982)","0.000 %
+ (0.965)","0.000 %
+ (0.950)"),pos=3,offset=0.2)
> legend("topleft",legend=c("Genome","ChIP (p-value)"),col=c("#5FA1C1","#EB9D86"),pch=15,bty="n")
> r0<-c(0.2034933016788197, 1.3978051793890356, 2.359553283752029, 19.734005184234114, 23.694856949054)
> r1<-c(0.0, 0.0, 0.0, 0.0, 0.0)
> height<-rbind(r0,r1)
> names=c("5'UTR","3'UTR","Coding Exon","Intron","All")
> mp<-barplot(height=height,names=names,beside=TRUE,horiz=FALSE,col=c("#5FA1C1","#EB9D86"),main="Gene",ylab="Percentage %",border=FALSE,ylim=c(0.000000,43.440571),cex.names=1)
> text(x=mp[1,],y=c(0.2034933016788197, 1.3978051793890356, 2.359553283752029, 19.734005184234114, 23.694856949054),label=c("0.2 %","1.4 %","2.4 %","19.7 %","23.7 %"),pos=3,offset=0.2)
> text(x=mp[2,],y=c(0.0, 0.0, 0.0, 0.0, 0.0),label=c("0.000 %
+ (0.998)","0.000 %
+ (0.986)","0.000 %
+ (0.976)","0.000 %
+ (0.803)","0.000 %
+ (0.763)"),pos=3,offset=0.2)
> legend("topleft",legend=c("Genome","ChIP (p-value)"),col=c("#5FA1C1","#EB9D86"),pch=15,bty="n")
> 
> # 
> # Distribution of Genome and ChIP regions over cis-regulatory element
> # Note that the x may be modified for better graphics in case a value is too small
> # Thus, look at the labels of the pie chart to get the real percentage values
> # 
> 
> par(mfcol=c(2, 2),mar=c(3, 3, 4, 2.8),oma=c(4, 2, 4, 2))
> x<-c(0.018532,0.017055,0.016037,0.017830,0.015092,0.014051,0.010000,0.013833,0.023014,0.192592,0.670292)
> pie(x=x,labels=c("1.9 %","1.7 %","1.6 %","1.8 %","1.5 %","1.4 %","0.2 %","1.4 %","2.3 %","19.3 %","67.0 %"),main="Genome",col=c("#445FA2","#EB9D86","#799F7A","#6C527F","#5FA1C1","#E8BB77","#A8C5EF","#FDCDB9","#C6E6B5","#F1D5EE","#B4E1F6"),clockwise=TRUE,border=FALSE,radius=0.9,cex=0.8,init.angle=90,density=100)
> x<-c(0.000000,1.000000)
> y<-c(0.000000,1.000000)
> plot(x, y,type="n",main="",xlab="",ylab="",frame=FALSE,axes=FALSE,xaxt="s",yaxt="s")
> legend("top",legend=c("Promoter (<=1000 bp): 1.9 %","Promoter (1000-2000 bp): 1.7 %","Promoter (2000-3000 bp): 1.6 %","Downstream (<=1000 bp): 1.8 %","Downstream (1000-2000 bp): 1.5 %","Downstream (2000-3000 bp): 1.4 %","5'UTR: 0.2 %","3'UTR: 1.4 %","Coding exon: 2.3 %","Intron: 19.3 %","Distal intergenic: 67.0 %"),col=c("#445FA2","#EB9D86","#799F7A","#6C527F","#5FA1C1","#E8BB77","#A8C5EF","#FDCDB9","#C6E6B5","#F1D5EE","#B4E1F6"),pch=15,bty="n")
> x<-c(0.010000,0.010000,0.010000,0.010000,0.010000,0.010000,0.010000,0.010000,0.010000,0.010000,1.000000)
> pie(x=x,labels=c("0.000 %","0.000 %","0.000 %","0.000 %","0.000 %","0.000 %","0.000 %","0.000 %","0.000 %","0.000 %","100.0 %"),main="ChIP",col=c("#445FA2","#EB9D86","#799F7A","#6C527F","#5FA1C1","#E8BB77","#A8C5EF","#FDCDB9","#C6E6B5","#F1D5EE","#B4E1F6"),clockwise=TRUE,border=FALSE,radius=0.9,cex=0.8,init.angle=90,density=100)
> x<-c(0.000000,1.000000)
> y<-c(0.000000,1.000000)
> plot(x, y,type="n",main="",xlab="",ylab="",frame=FALSE,axes=FALSE,xaxt="s",yaxt="s")
> legend("top",legend=c("Promoter (<=1000 bp): 0.000 %","Promoter (1000-2000 bp): 0.000 %","Promoter (2000-3000 bp): 0.000 %","Downstream (<=1000 bp): 0.000 %","Downstream (1000-2000 bp): 0.000 %","Downstream (2000-3000 bp): 0.000 %","5'UTR: 0.000 %","3'UTR: 0.000 %","Coding exon: 0.000 %","Intron: 0.000 %","Distal intergenic: 100.0 %"),col=c("#445FA2","#EB9D86","#799F7A","#6C527F","#5FA1C1","#E8BB77","#A8C5EF","#FDCDB9","#C6E6B5","#F1D5EE","#B4E1F6"),pch=15,bty="n")
> 
> # 
> # ChIP regions over the genome
> # 
> 
> par(mar=c(4, 4, 5, 3.8),oma=c(4, 2, 4, 2))
> layout(matrix(c(1, 0, 2, 2), 2, 2, byrow = TRUE),widths=c(1, 1),heights=c(1, 5))
> x<-c(0.000000,2.515610)
> y<-c(0.000000,1.000000)
> plot(x, y,type="n",main="Distribution of Peak Heights",xlab="",ylab="",xlim=c(0.000000,2.515610),ylim=c(0.000000,1.000000),frame=FALSE,xaxt="s",yaxt="n",cex=0.9)
> x<-c(0.000000,2.515610,2.515610,0.000000)
> y<-c(0.000000,0.000000,1.000000,1.000000)
> polygon(x,y,col=c("black"))
> x <- c(0.000000,0.169726,0.339451,0.509177,0.678903,0.848628,1.018354,1.188079,1.357805,1.527531,1.697256,1.866982,2.036708,2.206433,2.376159)
> y<-c(0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.800000)
> lines(x, y,xlim=c(0, 2.51561),ylim=c(0, 1),type="l",col=c("cyan"),lwd=2)
> x<-c(4119129.000000,4119130.000000)
> y<-c(0.855556,1.144444)
> plot(x, y,type="n",main="ChIP Regions (Peaks) over Chromosomes",xlab="Chromosome Size (bp)",ylab="Chromosome",xlim=c(4119129.000000,4119130.000000),ylim=c(0.855556,1.144444),frame=FALSE,xaxt="s",yaxt="n")
> start <- c(4119129)
> end <- c(4119130)
> vals <- c(2.51561)
> vals[vals > 2.51561] <- 2.51561
> vals[vals < 0] <- 0
> heights <- 0.288889 * ((vals - 0)/(2.51561 - 0)) + 0.855555555556
> for (i in 1:length(heights)) {
+ 	polygon(x=c(start[i], end[i], end[i], start[i]), y=c(0.855555555556, 0.855555555556, heights[i], heights[i]), col=c("#CC0000"), border=c("#CC0000"))
+ }
> mtext("26",side=2,line=0,outer=FALSE,at=1.0)
> dev.off()
null device 
          1 
> 
INFO  @ Tue, 23 Jun 2015 09:12:22: #... cong! See ceas.pdf for the graphical results of CEAS!