Mercurial > repos > geert-vandeweyer > coverage_report

comparison CoverageReport.pl @ 22:95062840f80f draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Correction to png calls to use cairo instead of x11. thanks to Eric Enns for pointing this out.
author geert-vandeweyer
date Mon, 17 Nov 2014 07:10:11 -0500
parents 86df3f847a72
children fd788f9db899
comparison
equal deleted inserted replaced
21:06ccf2048b8f 22:95062840f80f
206 my $spec = sprintf("%.1f",($ontarget / $totalmapped)*100); 206 my $spec = sprintf("%.1f",($ontarget / $totalmapped)*100);
207 # get min/max/median/average coverage => boxplot in R 207 # get min/max/median/average coverage => boxplot in R
208 open OUT, ">$wd/Rout/boxplot.R"; 208 open OUT, ">$wd/Rout/boxplot.R";
209 print OUT 'coverage <- read.table("../Targets.Global.Coverage",as.is=TRUE,sep="\t",header=FALSE)'."\n"; 209 print OUT 'coverage <- read.table("../Targets.Global.Coverage",as.is=TRUE,sep="\t",header=FALSE)'."\n";
210 print OUT 'coverage <- coverage[,'.$covcol.']'."\n"; 210 print OUT 'coverage <- coverage[,'.$covcol.']'."\n";
211 print OUT 'png(file="../Plots/CoverageBoxPlot.png", bg="white", width=240, height=480)'."\n"; 211 print OUT 'png(file="../Plots/CoverageBoxPlot.png", bg="white", width=240, height=480,type=c("cairo"))'."\n";
212 print OUT 'boxplot(coverage,range=1.5,main="Target Region Coverage")'."\n"; 212 print OUT 'boxplot(coverage,range=1.5,main="Target Region Coverage")'."\n";
213 print OUT 'graphics.off()'."\n"; 213 print OUT 'graphics.off()'."\n";
214 close OUT; 214 close OUT;
215 system("cd $wd/Rout && Rscript boxplot.R"); 215 system("cd $wd/Rout && Rscript boxplot.R");
216 216
272 print OUT ' i2 <- min(which(values$cov > mincov.x))'."\n"; 272 print OUT ' i2 <- min(which(values$cov > mincov.x))'."\n";
273 print OUT ' mincov.y <- ((values[i2,"count"] - values[i1,"count"])/(values[i2,"cov"] - values[i1,"cov"]))*(mincov.x - values[i1,"cov"]) + values[i1,"count"]'."\n"; 273 print OUT ' mincov.y <- ((values[i2,"count"] - values[i1,"count"])/(values[i2,"cov"] - values[i1,"cov"]))*(mincov.x - values[i1,"cov"]) + values[i1,"count"]'."\n";
274 print OUT ' }'."\n"; 274 print OUT ' }'."\n";
275 print OUT '}'."\n"; 275 print OUT '}'."\n";
276 # open output image and create plot 276 # open output image and create plot
277 print OUT 'png(file="../Plots/CoverageNtPlot.png", bg="white", width=540, height=480)'."\n"; 277 print OUT 'png(file="../Plots/CoverageNtPlot.png", bg="white", width=540, height=480,type=c("cairo"))'."\n";
278 print OUT 'par(xaxs="i",yaxs="i")'."\n"; 278 print OUT 'par(xaxs="i",yaxs="i")'."\n";
279 print OUT 'plot(values$cov,values$count,ylim=c(0,100),pch=".",main="Cumulative Normalised Base-Coverage Plot",xlab="Normalizalised Coverage",ylab="Cumulative Nr. Of Bases")'."\n"; 279 print OUT 'plot(values$cov,values$count,ylim=c(0,100),pch=".",main="Cumulative Normalised Base-Coverage Plot",xlab="Normalizalised Coverage",ylab="Cumulative Nr. Of Bases")'."\n";
280 print OUT 'lines(values$cov,values$count)'."\n"; 280 print OUT 'lines(values$cov,values$count)'."\n";
281 print OUT 'if (mincov.x <= 1) {'."\n"; 281 print OUT 'if (mincov.x <= 1) {'."\n";
282 print OUT ' lines(c(mincov.x,mincov.x),c(0,mincov.y),lty=2,col="darkgreen")'."\n"; 282 print OUT ' lines(c(mincov.x,mincov.x),c(0,mincov.y),lty=2,col="darkgreen")'."\n";
421 else { 421 else {
422 $scale = 1; 422 $scale = 1;
423 } 423 }
424 my $width = 480 * $scale; 424 my $width = 480 * $scale;
425 my $height = 240 * $scale; 425 my $height = 240 * $scale;
426 print OUT 'png(file="../Plots/Coverage_'.$currgroup.'.png", bg="white", width='.$width.', height='.$height.')'."\n"; 426 print OUT 'png(file="../Plots/Coverage_'.$currgroup.'.png", bg="white", width='.$width.', height='.$height.',type=c("cairo"))'."\n";
427 print OUT 'ylim = c(0,max(max(log10(coverage),log10('.($thresh+20).'))))'."\n"; 427 print OUT 'ylim = c(0,max(max(log10(coverage),log10('.($thresh+20).'))))'."\n";
428 print OUT 'mp <- barplot(log10(coverage),col=colors,main="Exon Coverage for '.$currgroup.'",ylab="Log10(Coverage)",ylim=ylim)'."\n"; 428 print OUT 'mp <- barplot(log10(coverage),col=colors,main="Exon Coverage for '.$currgroup.'",ylab="Log10(Coverage)",ylim=ylim)'."\n";
429 print OUT 'text(mp, log10(coverage) + '.(0.4/$scale).',format(coverage),xpd = TRUE,srt=90)'."\n"; 429 print OUT 'text(mp, log10(coverage) + '.(0.4/$scale).',format(coverage),xpd = TRUE,srt=90)'."\n";
430 print OUT 'text(mp,par("usr")[3]-0.05,labels=entries,srt=45,adj=1,xpd=TRUE)'."\n"; 430 print OUT 'text(mp,par("usr")[3]-0.05,labels=entries,srt=45,adj=1,xpd=TRUE)'."\n";
431 print OUT 'abline(h=log10('.$thresh.'),lwd=4,col=rgb(255,0,0,100,maxColorValue=255))'."\n"; 431 print OUT 'abline(h=log10('.$thresh.'),lwd=4,col=rgb(255,0,0,100,maxColorValue=255))'."\n";
463 else { 463 else {
464 $scale = 1; 464 $scale = 1;
465 } 465 }
466 my $width = 480 * $scale; 466 my $width = 480 * $scale;
467 my $height = 240 * $scale; 467 my $height = 240 * $scale;
468 print OUT 'png(file="../Plots/Coverage_'.$currgroup.'.png", bg="white", width='.$width.', height='.$height.')'."\n"; 468 print OUT 'png(file="../Plots/Coverage_'.$currgroup.'.png", bg="white", width='.$width.', height='.$height.',type=c("cairo"))'."\n";
469 print OUT 'ylim = c(0,max(max(log10(coverage),log10('.($thresh+20).'))))'."\n"; 469 print OUT 'ylim = c(0,max(max(log10(coverage),log10('.($thresh+20).'))))'."\n";
470 print OUT 'mp <- barplot(log10(coverage),col=colors,main="Exon Coverage for '.$currgroup.'",ylab="Log10(Coverage)", ylim=ylim)'."\n"; 470 print OUT 'mp <- barplot(log10(coverage),col=colors,main="Exon Coverage for '.$currgroup.'",ylab="Log10(Coverage)", ylim=ylim)'."\n";
471 print OUT 'text(mp, log10(coverage) + log10(2),format(coverage),xpd = TRUE,srt=90)'."\n"; 471 print OUT 'text(mp, log10(coverage) + log10(2),format(coverage),xpd = TRUE,srt=90)'."\n";
472 print OUT 'text(mp,par("usr")[3]-0.1,labels=entries,srt=45,adj=1,xpd=TRUE)'."\n"; 472 print OUT 'text(mp,par("usr")[3]-0.1,labels=entries,srt=45,adj=1,xpd=TRUE)'."\n";
473 print OUT 'abline(h=log10('.$thresh.'),lwd=4,col=rgb(255,0,0,100,maxColorValue=255))'."\n"; 473 print OUT 'abline(h=log10('.$thresh.'),lwd=4,col=rgb(255,0,0,100,maxColorValue=255))'."\n";
559 my $height = 240 ; 559 my $height = 240 ;
560 my $exonstr = $exon; 560 my $exonstr = $exon;
561 $exonstr =~ s/\s/_/g; 561 $exonstr =~ s/\s/_/g;
562 $exon =~ s/_/ /g; 562 $exon =~ s/_/ /g;
563 $exon =~ s/\|/ /g; 563 $exon =~ s/\|/ /g;
564 print OUT 'png(file="../Plots/Coverage_'.$exonstr.'.png", bg="white", width='.$width.', height='.$height.')'."\n"; 564 print OUT 'png(file="../Plots/Coverage_'.$exonstr.'.png", bg="white", width='.$width.', height='.$height.',type=c("cairo"))'."\n";
565 print OUT 'ylim = c(0,log10(max(max(coverage),'.($thresh+10).')))'."\n"; 565 print OUT 'ylim = c(0,log10(max(max(coverage),'.($thresh+10).')))'."\n";
566 if ($orient eq '-') { 566 if ($orient eq '-') {
567 print OUT 'plot(positions,log10(coverage),type="n",main="Coverage for '.$exon.'",ylab="log10(Coverage)",ylim=ylim,xlab="Position",xlim=rev(range(positions)),sub="(Transcribed from minus strand)")'."\n"; 567 print OUT 'plot(positions,log10(coverage),type="n",main="Coverage for '.$exon.'",ylab="log10(Coverage)",ylim=ylim,xlab="Position",xlim=rev(range(positions)),sub="(Transcribed from minus strand)")'."\n";
568 print OUT 'mtext("'.$subtitle.'")'."\n"; 568 print OUT 'mtext("'.$subtitle.'")'."\n";
569 } 569 }
656 my $height = 240 ; 656 my $height = 240 ;
657 my $exonstr = $exon; 657 my $exonstr = $exon;
658 $exonstr =~ s/\s/_/g; 658 $exonstr =~ s/\s/_/g;
659 $exon =~ s/_/ /g; 659 $exon =~ s/_/ /g;
660 $exon =~ s/\|/ /g; 660 $exon =~ s/\|/ /g;
661 print OUT 'png(file="../Plots/Coverage_'.$exonstr.'.png", bg="white", width='.$width.', height='.$height.')'."\n"; 661 print OUT 'png(file="../Plots/Coverage_'.$exonstr.'.png", bg="white", width='.$width.', height='.$height.',type=c("cairo"))'."\n";
662 print OUT 'ylim = c(0,log10(max(max(coverage),'.($thresh+10).')))'."\n"; 662 print OUT 'ylim = c(0,log10(max(max(coverage),'.($thresh+10).')))'."\n";
663 if ($orient eq '-') { 663 if ($orient eq '-') {
664 print OUT 'plot(positions,log10(coverage),type="n",main="Coverage for '.$exon.'",ylab="log10(Coverage)",ylim=ylim,xlab="Position",xlim=rev(range(positions)),sub="(Transcribed from minus strand)")'."\n"; 664 print OUT 'plot(positions,log10(coverage),type="n",main="Coverage for '.$exon.'",ylab="log10(Coverage)",ylim=ylim,xlab="Position",xlim=rev(range(positions)),sub="(Transcribed from minus strand)")'."\n";
665 print OUT 'mtext("'.$subtitle.'")'."\n"; 665 print OUT 'mtext("'.$subtitle.'")'."\n";
666 } 666 }