Mercurial > repos > pfrommolt > ngsrich

view 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
line wrap: on
line source

#!/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)

}
##########################################################################################################################################