Mercurial > repos > pfrommolt > ngsrich

comparison NGSrich_0.5.5/R/eval_details.R @ 0:89ad0a9cca52 default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author pfrommolt
date Mon, 21 Nov 2011 08:12:19 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:89ad0a9cca52
1 #!/usr/bin/Rscript
2
3
4 ##Input parameters########################################################################################################################
5 xml=as.character(commandArgs()[6])
6 bed=as.character(commandArgs()[7])
7 cov=as.character(commandArgs()[8])
8 chrinfo=as.character(commandArgs()[9])
9 outdir=as.character(commandArgs()[10])
10 genome=as.character(commandArgs()[11])
11 ##########################################################################################################################################
12
13
14
15 ##Read data###############################################################################################################################
16 bedfile<-read.table(bed,header=FALSE,stringsAsFactors=FALSE)
17 wunknown<-which(bedfile$V4=="unknown")
18 bedfile$V4[wunknown]<-paste("unknown",seq(1,length(wunknown),1),sep="")
19 bedfile$V4<-gsub("/","-",bedfile$V4)
20 covdata<-readLines(cov,warn=FALSE)
21
22 xmldata<-readLines(xml,warn=FALSE)
23 nbrline<-xmldata[grep("<NumberReads>",xmldata)]
24 rlenline<-xmldata[grep("<ReadLength>",xmldata)]
25 nbreads<-as.numeric(strsplit(strsplit(nbrline,">")[[1]][2],"<")[[1]][1])
26 rlength<-as.numeric(strsplit(strsplit(rlenline,">")[[1]][2],"<")[[1]][1])
27
28 expect<-0
29 allchr<-levels(as.factor(bedfile$V1))
30 if(chrinfo!="none"){
31 cinfo<-read.table(chrinfo,header=FALSE,stringsAsFactors=FALSE)
32 gsize<-sum(as.numeric(cinfo[,2]))
33 expect=(rlength*nbreads)/gsize
34 }
35
36 dir.create(paste(outdir,"/chromosomes",sep=""),showWarnings=FALSE)
37 for(chr in allchr){dir.create(paste(outdir,"/chromosomes/",chr,sep=""),showWarnings=FALSE)}
38
39 tcov<-list()
40 tstarts<-numeric(0)
41 tcov0<-character(0)
42 firstt=1
43 for(i in 1:length(covdata)){
44 if(i%%10000==0){cat(i,"/",length(covdata),"\n")}
45 if(nchar(covdata[i])>10){
46 if(firstt==0){
47 tcov[[tname]]<-as.numeric(tcov0)
48 tcov0<-character(0)
49 }
50 firstt=0
51 nameel<-strsplit(covdata[i],"\t")
52 tname<-paste(nameel[[1]][1],":",nameel[[1]][2],"-",nameel[[1]][3],sep="")
53 }
54 else{tcov0<-c(tcov0,covdata[i])}
55 }
56 ##########################################################################################################################################
57
58
59 sortchr<-function(x){
60 if(sum(substr(x,1,3)=="chr")==length(x)){
61 x0<-strsplit(x,"chr")
62 xspl<-character(0)
63 for(i in 1:length(x0)){xspl<-c(xspl,x0[[i]][2])}
64 xnum<-xspl[!is.na(as.numeric(xspl))]
65 xchar<-setdiff(xspl,xnum)
66 xsorted<-paste("chr",c(sort(as.numeric(xnum)),sort(xchar)),sep="")
67 return(xsorted)
68 }
69 else{return(x)}
70 }
71
72
73 ##Create gene index#######################################################################################################################
74 allchr<-sortchr(levels(as.factor(bedfile$V1)))
75 outfile=paste(outdir,"/chromosomes/chromosomes.html",sep="")
76 cat(file=outfile,paste(
77 "<html>\n<head>\n<title>Index of Target Regions</title>\n",
78 "<style type=\"text/css\">\n body{font-family:sans-serif;}\n h2,h3{color: darkblue;}\n a{color:darkblue;}\n table.output td{ padding: 4px; background-color: lightskyblue; border: 1px solid #000; border-color: darkblue; }\n</style>\n",
79 "</head>\n\n",
80 "<script language=\"JavaScript\">\n var questionClass=\"ChrList\";\n function collapseAll(){\n var allSections=document.getElementsByTagName(\"div\");\n for(i=0;i<allSections.length;i++){\n if(allSections[i].className==questionClass){\n allSections[i].style.display = \"none\";\n }}}\n function expand(name){\n collapseAll();\n var newStyle=\"\";\n if(document.getElementById(name).style.display!=\"block\"){newStyle=\"block\";}\n else{newStyle=\"none\";}\n document.getElementById(name).style.display=newStyle;\n }\n</script>\n\n",
81 "<body onload=\"expand('",allchr[1],"')\">\n<h2>Show Details on Genewise Coverage</h2>",sep=""))
82
83 for(chr in allchr){
84 cat(file=outfile,paste("<a href=\"javascript:expand('",chr,"')\">",chr,"</a>\n",sep=""),append=TRUE)
85 }
86 cat(file=outfile,"<br/><br/>\n\n",append=TRUE)
87
88 for(chr in allchr){
89
90 thischr<-bedfile[bedfile$V1==chr,]
91 thischr<-thischr[order(thischr$V2),]
92 allgenes<-intersect(thischr$V4,thischr$V4)
93 cat(file=outfile,paste("<div class=\"ChrList\" id=\"",chr,"\"><table class=\"output\">\n<tr><td><b>Gene</b></td><td><b>Region</b></td><td><b>Link</b></td></tr>\n",sep=""),append=TRUE)
94 for(gene in allgenes){
95 genebed<-thischr[thischr$V4==gene,]
96 firststart=genebed$V2[1]
97 lastend<-genebed$V3[nrow(genebed)]
98 tname=paste(chr,":",firststart,"-",lastend,sep="")
99 if(substr(gene,1,7)=="unknown"){genelab="unknown"}
100 else{genelab=gene}
101 cat(file=outfile,paste("<tr><td>",genelab,"</td><td>",tname,"</td><td><a href=\"",chr,"/",gene,"/index.html\">Show Details</a></td></tr>\n",sep=""),append=TRUE)
102 }
103 cat(file=outfile,"</table>\n</div>\n",append=TRUE)
104
105 }
106 cat(file=outfile,"</body></html>\n",append=TRUE)
107 ##########################################################################################################################################
108
109
110
111 ##Create genewise reports#################################################################################################################
112 firstgene=1
113 partest<-cv<-goftest<-numeric(0)
114 genebed<-split(bedfile,bedfile$V4)
115
116 for(gene in names(genebed)){
117
118 thisbed<-genebed[[gene]]
119 chr=thisbed$V1[1]
120
121 if(substr(gene,1,7)=="unknown"){genelab="unknown"}
122 else{genelab=gene}
123
124 ##Write HTML header
125 outfile=paste(outdir,"/chromosomes/",chr,"/",gene,"/index.html",sep="")
126 dir.create(paste(outdir,"/chromosomes/",chr,"/",gene,sep=""),showWarnings=FALSE)
127 cat(file=outfile,paste("<html>\n<head>\n<title>Enrichment Details</title>\n<style type=\"text/css\">\n body{font-family:sans-serif;}\n h2,h3{color: darkblue;}\n a{color:darkblue;}\n table.output td{ padding: 4px; background-color: lightskyblue; border: 1px solid #000; border-color: darkblue; } table.noborder td{padding: 0px; border: 0px solid #000; border-color: lightskyblue}\n</style>\n</head>\n\n<body>\n<h2>Enrichment Details for ",genelab,"</h2>\n<table class=\"output\">\n<tr><td><b>Target Region</b></td><td><b>Coverage</b></td><td><b>Percentage</b></td><td><b>Enrichment<br/>Factor</b></td><td><b>Target Histogram</b></td><td><b>External Link</b></td></tr>\n",sep=""))
128
129 for(i in 1:nrow(thisbed)){
130
131 start=thisbed$V2[i]
132 end=thisbed$V3[i]
133 region=seq(start,end,1)
134 tname=paste(chr,":",start,"-",end,sep="")
135
136 if(is.element(tname,names(tcov))){
137 basecov=tcov[[tname]]
138 if(expect!=0){enrratio<-round(mean(basecov)/expect,1)}
139 else{enrratio<-"unknown"}
140 stddev=round(sd(basecov),2)
141
142 ##Create plots
143 png(file=paste(outdir,"/chromosomes/",chr,"/",gene,"/h_",tname,".png",sep=""))
144 plot(region,basecov,pch=".",xlab=paste("Physical Position [",chr,"]",sep=""),ylab="Coverage",col="tomato3")
145 edgesx<-c(region,end,start)
146 edgesy<-c(basecov,0,0)
147 polygon(edgesx,edgesy,col="tomato3",border=NA)
148 lines(c(start,start),c(-1000,1000000))
149 lines(c(end,end),c(-1000,1000000))
150 lines(c(start,end),c(mean(basecov),mean(basecov)),col="darkblue")
151 lines(c(start,end),c(mean(basecov)-stddev,mean(basecov)-stddev),lty="dashed",col="darkblue")
152 lines(c(start,end),c(mean(basecov)+stddev,mean(basecov)+stddev),lty="dashed",col="darkblue")
153 dev.off()
154
155 ##Write HTML table
156 ucsclink<-paste("http://www.genome.ucsc.edu/cgi-bin/hgTracks?&db=",genome,"&position=",chr,"%3A",start,"-",end,"&hgt.suggest=&pix=800&Submit=submit&hgsid=183341879",sep="")
157 cat(file=outfile,paste("<tr><td><b>",tname,"</b><br/>(",end-start+1," bp)</td><td><table class=\"noborder\"><tr><td><b>Mean:</b></td><td>",thisbed[i,5],"</td></tr><tr><td><b>Std Dev:</b></td><td>",stddev,"</td></tr></table></td><td><table class=\"noborder\"><tr><td><b>1x:</b></td><td>",thisbed[i,6],"%</td></tr><tr><td><b>5x:</b></td><td>",thisbed[i,7],"%</td></tr><tr><td><b>10x:</b></td><td>",thisbed[i,8],"%</td></tr><tr><td><b>20x:</b></td><td>",thisbed[i,9],"%</td></tr><tr><td><b>30x:</b></td><td>",thisbed[i,10],"%</td></tr></table></td><td>",enrratio,"</td><td><a href=\"h_",tname,".png\"><img src=\"h_",tname,".png\" width=160></img></a></td><td><a href=\"",ucsclink,"\">Show in Genome Browser</td></tr>\n",sep=""),append=TRUE)
158 }
159 }
160 cat(file=outfile,"</table>\n</body>\n</html>",append=TRUE)
161
162 }
163 ##########################################################################################################################################