#!/usr/bin/Rscript

xmlfile=as.character(commandArgs()[6])
bedfile=as.character(commandArgs()[7])
outdir=as.character(commandArgs()[8])
genome=as.character(commandArgs()[9])
poor=as.numeric(commandArgs()[10])
high=as.numeric(commandArgs()[11])
samplename=as.character(commandArgs()[12])
targetfile=as.character(commandArgs()[13])
details=as.numeric(commandArgs()[14])
cutoff=0
sdetails=""
#samplename0=strsplit(strsplit(xmlfile,"_")[[1]][1],"/")[[1]]
#samplename<-samplename0[length(samplename0)]
outfile=paste(outdir,"/",samplename,"_enrichment.html",sep="")

##Read XML summary
cat("Reading XML file ... ")
xmltag<-function(line){
	return(strsplit(strsplit(line,">")[[1]][2],"</")[[1]][1])
}
xml<-readLines(xmlfile)
numreads<-numeric(0)
nextfold=0
for(line in xml){
  if(length(grep("ReadLength",line)>0)){
    readlength<-as.numeric(xmltag(line))
  }
  if(length(grep("NumberReads",line)>0)){
    numreads<-c(numreads,as.numeric(xmltag(line)))
  }
  if(length(grep("AvTargetCoverage",line)>0)){
    averagecov=as.numeric(xmltag(line))
  }
  if(details==1){
	sdetails="<a href=\"chromosomes/chromosomes.html\">[Show Details]</a><br/>"
  }
  if(length(grep("SDTargetCoverage",line)>0)){
    stddevcov=as.numeric(xmltag(line))
  }
  if(length(grep("TargetSize",line)>0)){
    targetsize=as.numeric(xmltag(line))
  }
  if(length(grep("<from1x>",line)>0)){nextfold=1}
  else{
    if((nextfold==1) && length(grep("PercBases",line))>0){
      sample1<-as.numeric(xmltag(line))
      nextfold=0
    }
  }
  if(length(grep("<from5x>",line)>0)){nextfold=5}
  else{
    if((nextfold==5) && length(grep("PercBases",line))>0){
      sample5<-as.numeric(xmltag(line))
      nextfold=0
    }
  }
  if(length(grep("<from10x>",line)>0)){nextfold=10}	
  else{
    if((nextfold==10) && length(grep("PercBases",line))>0){
      sample10<-as.numeric(xmltag(line))
      nextfold=0
    }
  }
  if(length(grep("<from20x>",line)>0)){nextfold=20}
  else{
    if((nextfold==20) && length(grep("PercBases",line))>0){
      sample20<-as.numeric(xmltag(line))
      nextfold=0
    }
  }
  if(length(grep("<from30x>",line)>0)){nextfold=30}
  else{
    if((nextfold==30) && length(grep("PercBases",line))>0){
      
      sample30<-as.numeric(xmltag(line))
      nextfold=0
    }
  }
}

numreads_total<-numreads[1]
numreads_target<-numreads[3]
tpkm=round(numreads_target/((targetsize/1000)*(numreads_total/1000000)),2)

cat("ready.\n")


##Read BED enrichment file and summarize output
cat("Reading BED file ... ")
bed<-read.table(bedfile,stringsAsFactors=FALSE)
area0_2<-sum(bed$V5<2)
area2_10<-sum((bed$V5>=2) & (bed$V5<10))
area10_20<-sum((bed$V5>=10) & (bed$V5<20))
area20_30<-sum((bed$V5>=20) & (bed$V5<30))
area30_50<-sum((bed$V5>=30) & (bed$V5<50))
area50_100<-sum((bed$V5>=50) & (bed$V5<100))
areagr100<-sum(bed$V5>100)
cat("ready.\n")


##Create pieplot
cat("Creating coverage pieplot ... ")
png(file=paste(outdir,"/plots/",samplename,"_pieplot.png",sep=""),width=580)
par(mar=c(1,7,1,7))
pie(c(area0_2,area2_10,area10_20,area20_30,area30_50,area50_100,areagr100),
    labels=c(paste("0x to 2x (",area0_2,")",sep=""),
      paste("2x to 10x (",area2_10,")",sep=""),
      paste("10x to 20x (",area10_20,")",sep=""),
      paste("20x to 30x (",area20_30,")",sep=""),
      paste("30x to 50x (",area30_50,")",sep=""),
      paste("50x to 100x (",area50_100,")",sep=""),
      paste("above 100x (",areagr100,")",sep="")),
    col=c("gray30","gray40","gray50","gray60","gray70","gray80","gray90"),main="") ##Mean Coverage of Target Regions")
garbage<-dev.off()
cat("ready.\n")


cat("Preparing coverage barplots ... ")
maxgenemean=0
genes<-levels(as.factor(bed$V4))
genes<-genes[genes!="unknown"]
ngenes<-length(genes)

for(i in 1:ngenes){
  genemean<-mean(bed[bed$V4==genes[i],5])
  if(maxgenemean<genemean){maxgenemean=genemean}
}

maxgenemean=maxgenemean-(maxgenemean%%100)+100
if((cutoff>0) && (maxgenemean>cutoff)){
  maxgenemean=cutoff
}
cat("ready.\n")

fwidth<-function(genes){
    fwidth0=10*length(genes)
    if(length(genes)<=100){fwidth0=10*length(genes)}
    if(length(genes)<=70){fwidth0=10*length(genes)}
    if(length(genes)<=40){fwidth0=25*length(genes)}
    if(length(genes)<=20){fwidth0=30*length(genes)}
    if(length(genes)<=10){fwidth0=35*length(genes)}
    if(length(genes)<=5){fwidth0=50*length(genes)}
    return(fwidth0)
}


chromosomes<-levels(as.factor(bed$V1))
if(length(genes)>=2000){
  cat("Creating coverage barplots ... ")
  for(chr in chromosomes){

    chrbed<-bed[bed$V1==chr,]
    genemean<-numeric(0)
    genes<-levels(as.factor(chrbed$V4))
    genes<-genes[genes!="unknown"]

    for(i in 1:length(genes)){
      if(cutoff>0){genemean[i]<-min(mean(chrbed[chrbed$V4==genes[i],5]),cutoff)}
      else{genemean[i]<-mean(chrbed[chrbed$V4==genes[i],5])}
    }

    png(file=paste(outdir,"/plots/",samplename,"_target_coverage_",chr,".png",sep=""),width=fwidth(genes),height=450)
    par(mar=c(7,5,1,2))
    barplot(as.numeric(genemean),names.arg=genes,las=2,
            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))
    lineh=0
    while(lineh<maxgenemean){
      lineh=lineh+100
      lines(c(-100,1000+length(genes)),c(lineh,lineh),col="gray")
    }
    garbage<-dev.off()
  }
  cat("ready.\n")
} else{
  cat("Creating coverage barplot ... ")
  for(i in 1:length(genes)){
    if(cutoff>0){genemean[i]<-min(mean(bed[bed$V4==genes[i],5]),cutoff)}
    else{genemean[i]<-mean(bed[bed$V4==genes[i],5])}
  }


  png(file=paste(outdir,"/plots/",samplename,"_target_coverage.png",sep=""),width=fwidth(genes),height=450)
  par(mar=c(7,5,1,2))
  barplot(as.numeric(genemean),names.arg=genes,las=2,
          ylab="Average Coverage at Gene Locus",col="tomato3",xlim=c(0,1.2*ngenes),ylim=c(0,maxgenemean))
  lineh=0
  while(lineh<maxgenemean){
    lineh=lineh+100
    lines(c(-100,1000+length(genes)),c(lineh,lineh),col="gray")
  }
  garbage<-dev.off()
  cat("ready.\n")
}


##Searching for poorly/highly covered genes
cat("Searching for poorly (<",poor,"x) and highly (>",high,"x) covered genes ... ",sep="")
genes<-levels(as.factor(bed$V4))
chr<-start<-end<-gene1<-ave<-stddev<-gr1x<-gr5x<-gr10x<-gr20x<-gr30x<-numeric(0)
for(gene in genes){

  bed0<-bed[bed$V4==gene,]
  len<-bed0[,3]-bed0[,2]
  fold1<-fold5<-fold10<-fold20<-fold30<-numeric(0)
 for(i in 1:nrow(bed0)){fold1<-c(fold1,as.numeric(bed0$V6[i]))}
 for(i in 1:nrow(bed0)){fold5<-c(fold5,as.numeric(bed0$V7[i]))}
 for(i in 1:nrow(bed0)){fold10<-c(fold10,as.numeric(bed0$V8[i]))}
 for(i in 1:nrow(bed0)){fold20<-c(fold20,as.numeric(bed0$V9[i]))}

 for(i in 1:nrow(bed0)){fold30<-c(fold30,as.numeric(bed0$V10[i]))}
  fold1_gene<-paste(round(sum(fold1*len)/sum(len),2),"%",sep="")
  fold5_gene<-paste(round(sum(fold5*len)/sum(len),2),"%",sep="")
  fold10_gene<-paste(round(sum(fold10*len)/sum(len),2),"%",sep="")
  fold20_gene<-paste(round(sum(fold20*len)/sum(len),2),"%",sep="")
  fold30_gene<-paste(round(sum(fold30*len)/sum(len),2),"%",sep="")
  
  chr<-c(chr,bed0[1,1])
  start<-c(start,min(bed0[,2]))
  end<-c(end,max(bed0[,3]))
  gene1<-c(gene1,gene)
  ave<-c(ave,round(mean(bed0[,5]),2))
  if(is.na(sd(bed0[,5]))){
    stddev<-c(stddev,0)
  }
  else{
    stddev<-c(stddev,round(sd(bed0[,5]),2))
  }
  
  gr1x<-c(gr1x,fold1_gene)
  gr5x<-c(gr5x,fold5_gene)
  gr10x<-c(gr10x,fold10_gene)
  gr20x<-c(gr20x,fold20_gene)
  gr30x<-c(gr30x,fold30_gene)
  
}

bed_gene<-data.frame(chr,start,end,gene1,ave,stddev,gr1x,gr5x,gr10x,gr20x,gr30x)

poorly_covered<-bed_gene[bed_gene$ave<poor,]
poorly_covered<-poorly_covered[order(poorly_covered[,5],decreasing=FALSE),]
highly_covered<-bed_gene[bed_gene$ave>high,]
highly_covered<-highly_covered[order(highly_covered[,5],decreasing=TRUE),]

linkpoor<-paste("http://www.genome.ucsc.edu/cgi-bin/hgTracks?&db=",genome,"&position=",
                poorly_covered$chr,"%3A",poorly_covered$start,"-",
                poorly_covered$end,"&hgt.suggest=&pix=800&Submit=submit&hgsid=183341879",sep="")
linkhigh<-paste("http://www.genome.ucsc.edu/cgi-bin/hgTracks?&db=",genome,"&position=",
                highly_covered$chr,"%3A",highly_covered$start,"-",
                highly_covered$end,"&hgt.suggest=&pix=800&Submit=submit&hgsid=183341879",sep="")
cat("ready.\n",sep="")


##Output HTML document
cat("Writing HTML output ... ");

cat(file=outfile,paste(
      "<html>\n",
      "<head>\n",
      "<title>Enrichment Performance</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",
      "\n",
      "<script language=\"JavaScript\">\n",
      "  var questionClass=\"chrView\";\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",
      "    }\n",
      "  }\n",
      "  function expand(name){\n",
      "    collapseAll();\n",
      "    var newStyle=\"\";\n",
      "    if(document.getElementById(name).style.display!=\"block\"){\n",
      "      newStyle=\"block\";\n",
      "    }\n",
      "    else{\n",
      "      newStyle=\"none\";\n",
      "    }\n",
      "    document.getElementById(name).style.display=newStyle;\n",
      "  }\n",
      "</script>\n",      
      "</head>\n",
      "\n",
      "<body onload=\"expand('",chromosomes[1],"')\">",
      "<h2>Enrichment Performance of Sample ",samplename,"</h2>\n",
      "<table>\n",
      "<tr>\n",
      "<td>\n",
      "<h3>Summary Statistics</h3>\n",
      "<table class=\"output\">\n",
      "  <tr><td><b># Reads</b></td><td>",numreads_total,"</td></tr>\n",
      "  <tr><td><b># On Target &plusmn 100 bp</b></td><td>",numreads_target,"</td></tr>\n",
      "  <tr><td><b>Target Size (bp)</b><td>",targetsize,"</td></tr>\n",
      "  <tr><td><b># Target Regions</b><td>",nrow(bed),"</td></tr>\n",
      "  <tr><td><b>Coverage Mean</b></td><td>",averagecov,"</td></tr>\n",
      "  <tr><td><b>Coverage Std Dev</b></td><td>",stddevcov,"</td></tr>\n",
      "  <tr><td><b>Covered 1x</b></td><td>",sample1,"</td></tr>\n",
      "  <tr><td><b>Covered 5x</b></td><td>",sample5,"</td></tr>\n",
      "  <tr><td><b>Covered 10x</b></td><td>",sample10,"</td></tr>\n",
      "  <tr><td><b>Covered 20x</b></td><td>",sample20,"</td></tr>\n",
      "  <tr><td><b>Covered 30x</b></td><td>",sample30,"</td></tr>\n",
      "  <tr><td><b>TPKM</b></td><td>",tpkm,"</td></tr>\n",
      "</table>\n",
      "</td>\n",
      "<td width=\"10%\"></td>",
      "<td>\n",
      "<img src=\"plots/",samplename,"_pieplot.png\"></img>\n",
      "</td>\n",
      "</tr>\n",
      "</table>\n",
      "<br/>\n",
      "<h2>Genewise Target Coverage</h2>",sdetails,"<br/>\n",sep=""))
      

if(ngenes>=2000){
for(chromosome in chromosomes){
  cat(file=outfile,paste("<a href=\"javascript:expand('",chromosome,"')\">",chromosome,"</a>\n",sep=""),append=TRUE)
}
for(chromosome in chromosomes){
  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)
}
} else{
  cat(file=outfile,paste("<div style=\"height:480px; overflow:auto;\"><img src=\"plots/",samplename,"_target_coverage.png\"></img></div>\n",sep=""),append=TRUE)
}

cat(file=outfile,paste("<br/><br/>\n",
      "<h2>Poorly Covered Genes (Cutoff: ",poor,"x)</h2>\n",sep=""),append=TRUE)
if(nrow(poorly_covered)==0){
	cat(file=outfile,"<p>Nothing found for this cutoff.</p>",append=TRUE)
} else{
cat(file=outfile,paste(
      "<table class=\"output\">\n",
      "  <tr>\n",
      "    <td><b>Region</b></td><td><b>Gene</b></td>\n",
      "    <td><b>Coverage Mean</b></td>\n",
      "    <td><b>Covered 1x</b></td><td><b>Covered 5x</b></td>\n",
      "    <td><b>Covered 10x</b></td><td><b>Covered 20x</b></td>\n",
      "    <td><b>Covered 30x</b></td><td><b>External Link</b></td>\n",
      "  </tr>\n",sep=""),append=TRUE)

for(i in 1:nrow(poorly_covered)){
cat(file=outfile,paste(
      "  <tr>\n",
      "    <td>",poorly_covered[i,1],":",poorly_covered[i,2],"-",poorly_covered[i,3],"</td>\n",
      "    <td>",poorly_covered[i,4],"</td>\n",
      "    <td>",poorly_covered[i,5],"</td>\n",
      "    <td>",poorly_covered[i,7],"</td>\n",
      "    <td>",poorly_covered[i,8],"</td>\n",
      "    <td>",poorly_covered[i,9],"</td>\n",
      "    <td>",poorly_covered[i,10],"</td>\n",
      "    <td>",poorly_covered[i,11],"</td>\n",
      "    <td><a href=\"",linkpoor[i],"\">Show in Genome Browser</a></td>\n",
      "  </tr>\n",sep=""),append=TRUE)
}

cat(file=outfile,"</table><br/>\n",append=TRUE)
}

cat(file=outfile,paste("<h2>Highly Covered Genes (Cutoff: ",high,"x)</h2>\n",sep=""),append=TRUE)
if(nrow(highly_covered)==0){
 cat(file=outfile,"<p>Nothing found for this cutoff.</p>",append=TRUE)
} else{
cat(file=outfile,paste(
      "<table class=\"output\">\n",
      "  <tr>\n",
      "    <td><b>Target</b></td><td><b>Gene</b></td><td><b>Coverage Mean</b></td>",
      "    <td><b>Covered 1x</b></td><td><b>Covered 5x</b></td><td><b>Covered 10x</b></td>",
      "    <td><b>Covered 20x</b></td><td><b>Covered 30x</b></td><td><b>External Link</b></td>\n",
      "  </tr>\n",sep=""),append=TRUE)
for(i in 1:nrow(highly_covered)){
cat(file=outfile,paste("  <tr>\n",
      "    <td>",highly_covered[i,1],":",highly_covered[i,2],"-",highly_covered[i,3],"</td>\n",
      "    <td>",highly_covered[i,4],"</td>\n",
      "    <td>",highly_covered[i,5],"</td>\n",
      "    <td>",highly_covered[i,7],"</td>\n",
      "    <td>",highly_covered[i,8],"</td>\n",
      "    <td>",highly_covered[i,9],"</td>\n",
      "    <td>",highly_covered[i,10],"</td>\n",
      "    <td>",highly_covered[i,11],"</td>\n",
      "    <td><a href=\"",linkhigh[i],"\">Show in Genome Browser</a></td>\n",
      "  </tr>",sep=""),append=TRUE)
}

cat(file=outfile,"</table>\n",append=TRUE)
}
cat(file=outfile,"<br/>BED file used for specification of target regions:<br/>",targetfile,append=TRUE)
cat(file=outfile,"</body>\n</html>\n",append=TRUE)

cat("ready.\n")