|
0
|
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 ##########################################################################################################################################
|
