Mercurial > repos > pfrommolt > ngsrich
diff NGSrich_0.5.5/R/eval_details.R @ 0:89ad0a9cca52 default tip
Uploaded
| author | pfrommolt |
|---|---|
| date | Mon, 21 Nov 2011 08:12:19 -0500 |
| parents | |
| children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/NGSrich_0.5.5/R/eval_details.R Mon Nov 21 08:12:19 2011 -0500 @@ -0,0 +1,163 @@ +#!/usr/bin/Rscript + + +##Input parameters######################################################################################################################## +xml=as.character(commandArgs()[6]) +bed=as.character(commandArgs()[7]) +cov=as.character(commandArgs()[8]) +chrinfo=as.character(commandArgs()[9]) +outdir=as.character(commandArgs()[10]) +genome=as.character(commandArgs()[11]) +########################################################################################################################################## + + + +##Read data############################################################################################################################### +bedfile<-read.table(bed,header=FALSE,stringsAsFactors=FALSE) +wunknown<-which(bedfile$V4=="unknown") +bedfile$V4[wunknown]<-paste("unknown",seq(1,length(wunknown),1),sep="") +bedfile$V4<-gsub("/","-",bedfile$V4) +covdata<-readLines(cov,warn=FALSE) + +xmldata<-readLines(xml,warn=FALSE) +nbrline<-xmldata[grep("<NumberReads>",xmldata)] +rlenline<-xmldata[grep("<ReadLength>",xmldata)] +nbreads<-as.numeric(strsplit(strsplit(nbrline,">")[[1]][2],"<")[[1]][1]) +rlength<-as.numeric(strsplit(strsplit(rlenline,">")[[1]][2],"<")[[1]][1]) + +expect<-0 +allchr<-levels(as.factor(bedfile$V1)) +if(chrinfo!="none"){ + cinfo<-read.table(chrinfo,header=FALSE,stringsAsFactors=FALSE) + gsize<-sum(as.numeric(cinfo[,2])) + expect=(rlength*nbreads)/gsize +} + +dir.create(paste(outdir,"/chromosomes",sep=""),showWarnings=FALSE) +for(chr in allchr){dir.create(paste(outdir,"/chromosomes/",chr,sep=""),showWarnings=FALSE)} + +tcov<-list() +tstarts<-numeric(0) +tcov0<-character(0) +firstt=1 +for(i in 1:length(covdata)){ + if(i%%10000==0){cat(i,"/",length(covdata),"\n")} + if(nchar(covdata[i])>10){ + if(firstt==0){ + tcov[[tname]]<-as.numeric(tcov0) + tcov0<-character(0) + } + firstt=0 + nameel<-strsplit(covdata[i],"\t") + tname<-paste(nameel[[1]][1],":",nameel[[1]][2],"-",nameel[[1]][3],sep="") + } + else{tcov0<-c(tcov0,covdata[i])} +} +########################################################################################################################################## + + +sortchr<-function(x){ + if(sum(substr(x,1,3)=="chr")==length(x)){ + x0<-strsplit(x,"chr") + xspl<-character(0) + for(i in 1:length(x0)){xspl<-c(xspl,x0[[i]][2])} + xnum<-xspl[!is.na(as.numeric(xspl))] + xchar<-setdiff(xspl,xnum) + xsorted<-paste("chr",c(sort(as.numeric(xnum)),sort(xchar)),sep="") + return(xsorted) + } + else{return(x)} +} + + +##Create gene index####################################################################################################################### +allchr<-sortchr(levels(as.factor(bedfile$V1))) +outfile=paste(outdir,"/chromosomes/chromosomes.html",sep="") +cat(file=outfile,paste( + "<html>\n<head>\n<title>Index of Target Regions</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; }\n</style>\n", + "</head>\n\n", + "<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", + "<body onload=\"expand('",allchr[1],"')\">\n<h2>Show Details on Genewise Coverage</h2>",sep="")) + +for(chr in allchr){ + cat(file=outfile,paste("<a href=\"javascript:expand('",chr,"')\">",chr,"</a>\n",sep=""),append=TRUE) +} +cat(file=outfile,"<br/><br/>\n\n",append=TRUE) + +for(chr in allchr){ + + thischr<-bedfile[bedfile$V1==chr,] + thischr<-thischr[order(thischr$V2),] + allgenes<-intersect(thischr$V4,thischr$V4) + 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) + for(gene in allgenes){ + genebed<-thischr[thischr$V4==gene,] + firststart=genebed$V2[1] + lastend<-genebed$V3[nrow(genebed)] + tname=paste(chr,":",firststart,"-",lastend,sep="") + if(substr(gene,1,7)=="unknown"){genelab="unknown"} + else{genelab=gene} + 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) + } + cat(file=outfile,"</table>\n</div>\n",append=TRUE) + +} +cat(file=outfile,"</body></html>\n",append=TRUE) +########################################################################################################################################## + + + +##Create genewise reports################################################################################################################# +firstgene=1 +partest<-cv<-goftest<-numeric(0) +genebed<-split(bedfile,bedfile$V4) + +for(gene in names(genebed)){ + + thisbed<-genebed[[gene]] + chr=thisbed$V1[1] + + if(substr(gene,1,7)=="unknown"){genelab="unknown"} + else{genelab=gene} + + ##Write HTML header + outfile=paste(outdir,"/chromosomes/",chr,"/",gene,"/index.html",sep="") + dir.create(paste(outdir,"/chromosomes/",chr,"/",gene,sep=""),showWarnings=FALSE) + 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="")) + + for(i in 1:nrow(thisbed)){ + + start=thisbed$V2[i] + end=thisbed$V3[i] + region=seq(start,end,1) + tname=paste(chr,":",start,"-",end,sep="") + + if(is.element(tname,names(tcov))){ + basecov=tcov[[tname]] + if(expect!=0){enrratio<-round(mean(basecov)/expect,1)} + else{enrratio<-"unknown"} + stddev=round(sd(basecov),2) + + ##Create plots + png(file=paste(outdir,"/chromosomes/",chr,"/",gene,"/h_",tname,".png",sep="")) + plot(region,basecov,pch=".",xlab=paste("Physical Position [",chr,"]",sep=""),ylab="Coverage",col="tomato3") + edgesx<-c(region,end,start) + edgesy<-c(basecov,0,0) + polygon(edgesx,edgesy,col="tomato3",border=NA) + lines(c(start,start),c(-1000,1000000)) + lines(c(end,end),c(-1000,1000000)) + lines(c(start,end),c(mean(basecov),mean(basecov)),col="darkblue") + lines(c(start,end),c(mean(basecov)-stddev,mean(basecov)-stddev),lty="dashed",col="darkblue") + lines(c(start,end),c(mean(basecov)+stddev,mean(basecov)+stddev),lty="dashed",col="darkblue") + dev.off() + + ##Write HTML table + 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="") + 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) + } + } + cat(file=outfile,"</table>\n</body>\n</html>",append=TRUE) + +} +##########################################################################################################################################