annotate NGSrich_0.5.5/R/eval_enrichment.R @ 0:89ad0a9cca52 default tip

Uploaded
author pfrommolt
date Mon, 21 Nov 2011 08:12:19 -0500
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
1 #!/usr/bin/Rscript
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
2
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
3 xmlfile=as.character(commandArgs()[6])
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
4 bedfile=as.character(commandArgs()[7])
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
5 outdir=as.character(commandArgs()[8])
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
6 genome=as.character(commandArgs()[9])
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
7 poor=as.numeric(commandArgs()[10])
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
8 high=as.numeric(commandArgs()[11])
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
9 samplename=as.character(commandArgs()[12])
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
10 targetfile=as.character(commandArgs()[13])
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
11 details=as.numeric(commandArgs()[14])
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
12 cutoff=0
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
13 sdetails=""
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
14 #samplename0=strsplit(strsplit(xmlfile,"_")[[1]][1],"/")[[1]]
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
15 #samplename<-samplename0[length(samplename0)]
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
16 outfile=paste(outdir,"/",samplename,"_enrichment.html",sep="")
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
17
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
18 ##Read XML summary
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
19 cat("Reading XML file ... ")
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
20 xmltag<-function(line){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
21 return(strsplit(strsplit(line,">")[[1]][2],"</")[[1]][1])
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
22 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
23 xml<-readLines(xmlfile)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
24 numreads<-numeric(0)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
25 nextfold=0
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
26 for(line in xml){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
27 if(length(grep("ReadLength",line)>0)){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
28 readlength<-as.numeric(xmltag(line))
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
29 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
30 if(length(grep("NumberReads",line)>0)){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
31 numreads<-c(numreads,as.numeric(xmltag(line)))
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
32 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
33 if(length(grep("AvTargetCoverage",line)>0)){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
34 averagecov=as.numeric(xmltag(line))
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
35 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
36 if(details==1){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
37 sdetails="<a href=\"chromosomes/chromosomes.html\">[Show Details]</a><br/>"
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
38 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
39 if(length(grep("SDTargetCoverage",line)>0)){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
40 stddevcov=as.numeric(xmltag(line))
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
41 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
42 if(length(grep("TargetSize",line)>0)){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
43 targetsize=as.numeric(xmltag(line))
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
44 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
45 if(length(grep("<from1x>",line)>0)){nextfold=1}
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
46 else{
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
47 if((nextfold==1) && length(grep("PercBases",line))>0){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
48 sample1<-as.numeric(xmltag(line))
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
49 nextfold=0
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
50 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
51 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
52 if(length(grep("<from5x>",line)>0)){nextfold=5}
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
53 else{
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
54 if((nextfold==5) && length(grep("PercBases",line))>0){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
55 sample5<-as.numeric(xmltag(line))
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
56 nextfold=0
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
57 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
58 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
59 if(length(grep("<from10x>",line)>0)){nextfold=10}
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
60 else{
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
61 if((nextfold==10) && length(grep("PercBases",line))>0){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
62 sample10<-as.numeric(xmltag(line))
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
63 nextfold=0
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
64 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
65 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
66 if(length(grep("<from20x>",line)>0)){nextfold=20}
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
67 else{
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
68 if((nextfold==20) && length(grep("PercBases",line))>0){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
69 sample20<-as.numeric(xmltag(line))
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
70 nextfold=0
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
71 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
72 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
73 if(length(grep("<from30x>",line)>0)){nextfold=30}
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
74 else{
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
75 if((nextfold==30) && length(grep("PercBases",line))>0){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
76
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
77 sample30<-as.numeric(xmltag(line))
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
78 nextfold=0
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
79 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
80 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
81 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
82
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
83 numreads_total<-numreads[1]
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
84 numreads_target<-numreads[3]
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
85 tpkm=round(numreads_target/((targetsize/1000)*(numreads_total/1000000)),2)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
86
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
87 cat("ready.\n")
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
88
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
89
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
90 ##Read BED enrichment file and summarize output
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
91 cat("Reading BED file ... ")
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
92 bed<-read.table(bedfile,stringsAsFactors=FALSE)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
93 area0_2<-sum(bed$V5<2)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
94 area2_10<-sum((bed$V5>=2) & (bed$V5<10))
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
95 area10_20<-sum((bed$V5>=10) & (bed$V5<20))
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
96 area20_30<-sum((bed$V5>=20) & (bed$V5<30))
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
97 area30_50<-sum((bed$V5>=30) & (bed$V5<50))
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
98 area50_100<-sum((bed$V5>=50) & (bed$V5<100))
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
99 areagr100<-sum(bed$V5>100)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
100 cat("ready.\n")
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
101
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
102
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
103 ##Create pieplot
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
104 cat("Creating coverage pieplot ... ")
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
105 png(file=paste(outdir,"/plots/",samplename,"_pieplot.png",sep=""),width=580)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
106 par(mar=c(1,7,1,7))
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
107 pie(c(area0_2,area2_10,area10_20,area20_30,area30_50,area50_100,areagr100),
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
108 labels=c(paste("0x to 2x (",area0_2,")",sep=""),
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
109 paste("2x to 10x (",area2_10,")",sep=""),
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
110 paste("10x to 20x (",area10_20,")",sep=""),
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
111 paste("20x to 30x (",area20_30,")",sep=""),
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
112 paste("30x to 50x (",area30_50,")",sep=""),
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
113 paste("50x to 100x (",area50_100,")",sep=""),
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
114 paste("above 100x (",areagr100,")",sep="")),
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
115 col=c("gray30","gray40","gray50","gray60","gray70","gray80","gray90"),main="") ##Mean Coverage of Target Regions")
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
116 garbage<-dev.off()
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
117 cat("ready.\n")
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
118
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
119
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
120 cat("Preparing coverage barplots ... ")
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
121 maxgenemean=0
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
122 genes<-levels(as.factor(bed$V4))
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
123 genes<-genes[genes!="unknown"]
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
124 ngenes<-length(genes)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
125
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
126 for(i in 1:ngenes){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
127 genemean<-mean(bed[bed$V4==genes[i],5])
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
128 if(maxgenemean<genemean){maxgenemean=genemean}
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
129 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
130
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
131 maxgenemean=maxgenemean-(maxgenemean%%100)+100
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
132 if((cutoff>0) && (maxgenemean>cutoff)){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
133 maxgenemean=cutoff
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
134 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
135 cat("ready.\n")
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
136
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
137 fwidth<-function(genes){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
138 fwidth0=10*length(genes)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
139 if(length(genes)<=100){fwidth0=10*length(genes)}
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
140 if(length(genes)<=70){fwidth0=10*length(genes)}
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
141 if(length(genes)<=40){fwidth0=25*length(genes)}
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
142 if(length(genes)<=20){fwidth0=30*length(genes)}
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
143 if(length(genes)<=10){fwidth0=35*length(genes)}
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
144 if(length(genes)<=5){fwidth0=50*length(genes)}
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
145 return(fwidth0)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
146 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
147
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
148
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
149 chromosomes<-levels(as.factor(bed$V1))
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
150 if(length(genes)>=2000){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
151 cat("Creating coverage barplots ... ")
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
152 for(chr in chromosomes){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
153
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
154 chrbed<-bed[bed$V1==chr,]
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
155 genemean<-numeric(0)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
156 genes<-levels(as.factor(chrbed$V4))
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
157 genes<-genes[genes!="unknown"]
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
158
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
159 for(i in 1:length(genes)){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
160 if(cutoff>0){genemean[i]<-min(mean(chrbed[chrbed$V4==genes[i],5]),cutoff)}
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
161 else{genemean[i]<-mean(chrbed[chrbed$V4==genes[i],5])}
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
162 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
163
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
164 png(file=paste(outdir,"/plots/",samplename,"_target_coverage_",chr,".png",sep=""),width=fwidth(genes),height=450)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
165 par(mar=c(7,5,1,2))
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
166 barplot(as.numeric(genemean),names.arg=genes,las=2,
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
167 ylab=paste("Average Coverage at Gene Locus (",chr,")",sep=""),col="tomato3",xlim=c(0.043*length(genemean),1.2*length(genes)-0.043*length(genemean)),ylim=c(0,maxgenemean))
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
168 lineh=0
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
169 while(lineh<maxgenemean){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
170 lineh=lineh+100
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
171 lines(c(-100,1000+length(genes)),c(lineh,lineh),col="gray")
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
172 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
173 garbage<-dev.off()
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
174 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
175 cat("ready.\n")
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
176 } else{
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
177 cat("Creating coverage barplot ... ")
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
178 for(i in 1:length(genes)){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
179 if(cutoff>0){genemean[i]<-min(mean(bed[bed$V4==genes[i],5]),cutoff)}
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
180 else{genemean[i]<-mean(bed[bed$V4==genes[i],5])}
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
181 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
182
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
183
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
184 png(file=paste(outdir,"/plots/",samplename,"_target_coverage.png",sep=""),width=fwidth(genes),height=450)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
185 par(mar=c(7,5,1,2))
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
186 barplot(as.numeric(genemean),names.arg=genes,las=2,
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
187 ylab="Average Coverage at Gene Locus",col="tomato3",xlim=c(0,1.2*ngenes),ylim=c(0,maxgenemean))
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
188 lineh=0
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
189 while(lineh<maxgenemean){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
190 lineh=lineh+100
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
191 lines(c(-100,1000+length(genes)),c(lineh,lineh),col="gray")
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
192 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
193 garbage<-dev.off()
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
194 cat("ready.\n")
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
195 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
196
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
197
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
198 ##Searching for poorly/highly covered genes
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
199 cat("Searching for poorly (<",poor,"x) and highly (>",high,"x) covered genes ... ",sep="")
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
200 genes<-levels(as.factor(bed$V4))
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
201 chr<-start<-end<-gene1<-ave<-stddev<-gr1x<-gr5x<-gr10x<-gr20x<-gr30x<-numeric(0)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
202 for(gene in genes){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
203
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
204 bed0<-bed[bed$V4==gene,]
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
205 len<-bed0[,3]-bed0[,2]
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
206 fold1<-fold5<-fold10<-fold20<-fold30<-numeric(0)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
207 for(i in 1:nrow(bed0)){fold1<-c(fold1,as.numeric(bed0$V6[i]))}
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
208 for(i in 1:nrow(bed0)){fold5<-c(fold5,as.numeric(bed0$V7[i]))}
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
209 for(i in 1:nrow(bed0)){fold10<-c(fold10,as.numeric(bed0$V8[i]))}
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
210 for(i in 1:nrow(bed0)){fold20<-c(fold20,as.numeric(bed0$V9[i]))}
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
211
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
212 for(i in 1:nrow(bed0)){fold30<-c(fold30,as.numeric(bed0$V10[i]))}
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
213 fold1_gene<-paste(round(sum(fold1*len)/sum(len),2),"%",sep="")
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
214 fold5_gene<-paste(round(sum(fold5*len)/sum(len),2),"%",sep="")
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
215 fold10_gene<-paste(round(sum(fold10*len)/sum(len),2),"%",sep="")
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
216 fold20_gene<-paste(round(sum(fold20*len)/sum(len),2),"%",sep="")
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
217 fold30_gene<-paste(round(sum(fold30*len)/sum(len),2),"%",sep="")
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
218
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
219 chr<-c(chr,bed0[1,1])
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
220 start<-c(start,min(bed0[,2]))
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
221 end<-c(end,max(bed0[,3]))
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
222 gene1<-c(gene1,gene)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
223 ave<-c(ave,round(mean(bed0[,5]),2))
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
224 if(is.na(sd(bed0[,5]))){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
225 stddev<-c(stddev,0)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
226 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
227 else{
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
228 stddev<-c(stddev,round(sd(bed0[,5]),2))
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
229 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
230
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
231 gr1x<-c(gr1x,fold1_gene)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
232 gr5x<-c(gr5x,fold5_gene)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
233 gr10x<-c(gr10x,fold10_gene)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
234 gr20x<-c(gr20x,fold20_gene)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
235 gr30x<-c(gr30x,fold30_gene)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
236
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
237 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
238
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
239 bed_gene<-data.frame(chr,start,end,gene1,ave,stddev,gr1x,gr5x,gr10x,gr20x,gr30x)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
240
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
241 poorly_covered<-bed_gene[bed_gene$ave<poor,]
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
242 poorly_covered<-poorly_covered[order(poorly_covered[,5],decreasing=FALSE),]
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
243 highly_covered<-bed_gene[bed_gene$ave>high,]
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
244 highly_covered<-highly_covered[order(highly_covered[,5],decreasing=TRUE),]
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
245
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
246 linkpoor<-paste("http://www.genome.ucsc.edu/cgi-bin/hgTracks?&db=",genome,"&position=",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
247 poorly_covered$chr,"%3A",poorly_covered$start,"-",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
248 poorly_covered$end,"&hgt.suggest=&pix=800&Submit=submit&hgsid=183341879",sep="")
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
249 linkhigh<-paste("http://www.genome.ucsc.edu/cgi-bin/hgTracks?&db=",genome,"&position=",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
250 highly_covered$chr,"%3A",highly_covered$start,"-",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
251 highly_covered$end,"&hgt.suggest=&pix=800&Submit=submit&hgsid=183341879",sep="")
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
252 cat("ready.\n",sep="")
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
253
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
254
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
255 ##Output HTML document
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
256 cat("Writing HTML output ... ");
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
257
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
258 cat(file=outfile,paste(
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
259 "<html>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
260 "<head>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
261 "<title>Enrichment Performance</title>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
262 "<style type=\"text/css\">\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
263 " body{font-family:sans-serif;}\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
264 " h2,h3{color: darkblue;}\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
265 " a{color:darkblue;}\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
266 " table.output td{",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
267 " padding: 4px; background-color: lightskyblue;",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
268 " border: 1px solid #000; border-color: darkblue;",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
269 " }\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
270 "</style>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
271 "\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
272 "<script language=\"JavaScript\">\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
273 " var questionClass=\"chrView\";\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
274 " function collapseAll(){\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
275 " var allSections = document.getElementsByTagName(\"div\");\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
276 " for(i=0; i<allSections.length; i++){\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
277 " if(allSections[i].className==questionClass){\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
278 " allSections[i].style.display=\"none\";\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
279 " }\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
280 " }\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
281 " }\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
282 " function expand(name){\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
283 " collapseAll();\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
284 " var newStyle=\"\";\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
285 " if(document.getElementById(name).style.display!=\"block\"){\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
286 " newStyle=\"block\";\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
287 " }\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
288 " else{\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
289 " newStyle=\"none\";\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
290 " }\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
291 " document.getElementById(name).style.display=newStyle;\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
292 " }\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
293 "</script>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
294 "</head>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
295 "\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
296 "<body onload=\"expand('",chromosomes[1],"')\">",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
297 "<h2>Enrichment Performance of Sample ",samplename,"</h2>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
298 "<table>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
299 "<tr>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
300 "<td>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
301 "<h3>Summary Statistics</h3>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
302 "<table class=\"output\">\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
303 " <tr><td><b># Reads</b></td><td>",numreads_total,"</td></tr>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
304 " <tr><td><b># On Target &plusmn 100 bp</b></td><td>",numreads_target,"</td></tr>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
305 " <tr><td><b>Target Size (bp)</b><td>",targetsize,"</td></tr>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
306 " <tr><td><b># Target Regions</b><td>",nrow(bed),"</td></tr>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
307 " <tr><td><b>Coverage Mean</b></td><td>",averagecov,"</td></tr>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
308 " <tr><td><b>Coverage Std Dev</b></td><td>",stddevcov,"</td></tr>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
309 " <tr><td><b>Covered 1x</b></td><td>",sample1,"</td></tr>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
310 " <tr><td><b>Covered 5x</b></td><td>",sample5,"</td></tr>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
311 " <tr><td><b>Covered 10x</b></td><td>",sample10,"</td></tr>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
312 " <tr><td><b>Covered 20x</b></td><td>",sample20,"</td></tr>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
313 " <tr><td><b>Covered 30x</b></td><td>",sample30,"</td></tr>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
314 " <tr><td><b>TPKM</b></td><td>",tpkm,"</td></tr>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
315 "</table>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
316 "</td>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
317 "<td width=\"10%\"></td>",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
318 "<td>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
319 "<img src=\"plots/",samplename,"_pieplot.png\"></img>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
320 "</td>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
321 "</tr>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
322 "</table>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
323 "<br/>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
324 "<h2>Genewise Target Coverage</h2>",sdetails,"<br/>\n",sep=""))
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
325
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
326
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
327 if(ngenes>=2000){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
328 for(chromosome in chromosomes){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
329 cat(file=outfile,paste("<a href=\"javascript:expand('",chromosome,"')\">",chromosome,"</a>\n",sep=""),append=TRUE)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
330 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
331 for(chromosome in chromosomes){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
332 cat(file=outfile,paste("<div style=\"height:480px; overflow:auto;\" class=\"chrView\" id=\"",chromosome,"\"><img src=\"plots/",samplename,"_target_coverage_",chromosome,".png\"></img></div>\n",sep=""),append=TRUE)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
333 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
334 } else{
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
335 cat(file=outfile,paste("<div style=\"height:480px; overflow:auto;\"><img src=\"plots/",samplename,"_target_coverage.png\"></img></div>\n",sep=""),append=TRUE)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
336 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
337
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
338 cat(file=outfile,paste("<br/><br/>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
339 "<h2>Poorly Covered Genes (Cutoff: ",poor,"x)</h2>\n",sep=""),append=TRUE)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
340 if(nrow(poorly_covered)==0){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
341 cat(file=outfile,"<p>Nothing found for this cutoff.</p>",append=TRUE)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
342 } else{
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
343 cat(file=outfile,paste(
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
344 "<table class=\"output\">\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
345 " <tr>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
346 " <td><b>Region</b></td><td><b>Gene</b></td>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
347 " <td><b>Coverage Mean</b></td>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
348 " <td><b>Covered 1x</b></td><td><b>Covered 5x</b></td>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
349 " <td><b>Covered 10x</b></td><td><b>Covered 20x</b></td>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
350 " <td><b>Covered 30x</b></td><td><b>External Link</b></td>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
351 " </tr>\n",sep=""),append=TRUE)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
352
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
353 for(i in 1:nrow(poorly_covered)){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
354 cat(file=outfile,paste(
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
355 " <tr>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
356 " <td>",poorly_covered[i,1],":",poorly_covered[i,2],"-",poorly_covered[i,3],"</td>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
357 " <td>",poorly_covered[i,4],"</td>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
358 " <td>",poorly_covered[i,5],"</td>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
359 " <td>",poorly_covered[i,7],"</td>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
360 " <td>",poorly_covered[i,8],"</td>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
361 " <td>",poorly_covered[i,9],"</td>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
362 " <td>",poorly_covered[i,10],"</td>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
363 " <td>",poorly_covered[i,11],"</td>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
364 " <td><a href=\"",linkpoor[i],"\">Show in Genome Browser</a></td>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
365 " </tr>\n",sep=""),append=TRUE)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
366 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
367
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
368 cat(file=outfile,"</table><br/>\n",append=TRUE)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
369 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
370
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
371 cat(file=outfile,paste("<h2>Highly Covered Genes (Cutoff: ",high,"x)</h2>\n",sep=""),append=TRUE)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
372 if(nrow(highly_covered)==0){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
373 cat(file=outfile,"<p>Nothing found for this cutoff.</p>",append=TRUE)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
374 } else{
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
375 cat(file=outfile,paste(
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
376 "<table class=\"output\">\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
377 " <tr>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
378 " <td><b>Target</b></td><td><b>Gene</b></td><td><b>Coverage Mean</b></td>",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
379 " <td><b>Covered 1x</b></td><td><b>Covered 5x</b></td><td><b>Covered 10x</b></td>",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
380 " <td><b>Covered 20x</b></td><td><b>Covered 30x</b></td><td><b>External Link</b></td>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
381 " </tr>\n",sep=""),append=TRUE)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
382 for(i in 1:nrow(highly_covered)){
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
383 cat(file=outfile,paste(" <tr>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
384 " <td>",highly_covered[i,1],":",highly_covered[i,2],"-",highly_covered[i,3],"</td>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
385 " <td>",highly_covered[i,4],"</td>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
386 " <td>",highly_covered[i,5],"</td>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
387 " <td>",highly_covered[i,7],"</td>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
388 " <td>",highly_covered[i,8],"</td>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
389 " <td>",highly_covered[i,9],"</td>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
390 " <td>",highly_covered[i,10],"</td>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
391 " <td>",highly_covered[i,11],"</td>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
392 " <td><a href=\"",linkhigh[i],"\">Show in Genome Browser</a></td>\n",
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
393 " </tr>",sep=""),append=TRUE)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
394 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
395
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
396 cat(file=outfile,"</table>\n",append=TRUE)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
397 }
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
398 cat(file=outfile,"<br/>BED file used for specification of target regions:<br/>",targetfile,append=TRUE)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
399 cat(file=outfile,"</body>\n</html>\n",append=TRUE)
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
400
89ad0a9cca52 Uploaded
pfrommolt
parents:
diff changeset
401 cat("ready.\n")