comparison CoverageReport.pl @ 22:95062840f80f draft

Correction to png calls to use cairo instead of x11. thanks to Eric Enns for pointing this out.

21:06ccf2048b8f

my $spec = sprintf("%.1f",($ontarget / $totalmapped)*100);
# get min/max/median/average coverage => boxplot in R
open OUT, ">$wd/Rout/boxplot.R";
print OUT 'coverage <- read.table("../Targets.Global.Coverage",as.is=TRUE,sep="\t",header=FALSE)'."\n";
print OUT 'coverage <- coverage[,'.$covcol.']'."\n";
print OUT 'png(file="../Plots/CoverageBoxPlot.png", bg="white", width=240, height=480)'."\n";
print OUT 'boxplot(coverage,range=1.5,main="Target Region Coverage")'."\n";
print OUT 'graphics.off()'."\n";
close OUT;
system("cd $wd/Rout && Rscript boxplot.R");

print OUT '    i2 <- min(which(values$cov > mincov.x))'."\n";
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";
print OUT '  }'."\n";
print OUT '}'."\n";
# open output image and create plot
print OUT 'png(file="../Plots/CoverageNtPlot.png", bg="white", width=540, height=480)'."\n";
print OUT 'par(xaxs="i",yaxs="i")'."\n";
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";
print OUT 'lines(values$cov,values$count)'."\n";
print OUT 'if (mincov.x <= 1) {'."\n";
print OUT '  lines(c(mincov.x,mincov.x),c(0,mincov.y),lty=2,col="darkgreen")'."\n";

			else {
				$scale = 1;
			}
			my $width = 480 * $scale;
			my $height = 240 * $scale;
			print OUT 'png(file="../Plots/Coverage_'.$currgroup.'.png", bg="white", width='.$width.', height='.$height.')'."\n";
			print OUT 'ylim = c(0,max(max(log10(coverage),log10('.($thresh+20).'))))'."\n";
			print OUT 'mp <- barplot(log10(coverage),col=colors,main="Exon Coverage for '.$currgroup.'",ylab="Log10(Coverage)",ylim=ylim)'."\n";
			print OUT 'text(mp, log10(coverage) + '.(0.4/$scale).',format(coverage),xpd = TRUE,srt=90)'."\n";
			print OUT 'text(mp,par("usr")[3]-0.05,labels=entries,srt=45,adj=1,xpd=TRUE)'."\n";
			print OUT 'abline(h=log10('.$thresh.'),lwd=4,col=rgb(255,0,0,100,maxColorValue=255))'."\n";

		else {
			$scale = 1;
		}
		my $width = 480 * $scale;
		my $height = 240 * $scale;
		print OUT 'png(file="../Plots/Coverage_'.$currgroup.'.png", bg="white", width='.$width.', height='.$height.')'."\n";
		print OUT 'ylim = c(0,max(max(log10(coverage),log10('.($thresh+20).'))))'."\n";
		print OUT 'mp <- barplot(log10(coverage),col=colors,main="Exon Coverage for '.$currgroup.'",ylab="Log10(Coverage)", ylim=ylim)'."\n";
		print OUT 'text(mp, log10(coverage) + log10(2),format(coverage),xpd = TRUE,srt=90)'."\n";
		print OUT 'text(mp,par("usr")[3]-0.1,labels=entries,srt=45,adj=1,xpd=TRUE)'."\n";
		print OUT 'abline(h=log10('.$thresh.'),lwd=4,col=rgb(255,0,0,100,maxColorValue=255))'."\n";

		my $height = 240 ;
		my $exonstr = $exon;
		$exonstr =~ s/\s/_/g;
		$exon =~ s/_/ /g;
		$exon =~ s/\|/ /g;
		print OUT 'png(file="../Plots/Coverage_'.$exonstr.'.png", bg="white", width='.$width.', height='.$height.')'."\n";
		print OUT 'ylim = c(0,log10(max(max(coverage),'.($thresh+10).')))'."\n";
		if ($orient eq '-') {
			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";
			print OUT 'mtext("'.$subtitle.'")'."\n";
		}

		my $height = 240 ;
		my $exonstr = $exon;
		$exonstr =~ s/\s/_/g;
		$exon =~ s/_/ /g;
		$exon =~ s/\|/ /g;
		print OUT 'png(file="../Plots/Coverage_'.$exonstr.'.png", bg="white", width='.$width.', height='.$height.')'."\n";
		print OUT 'ylim = c(0,log10(max(max(coverage),'.($thresh+10).')))'."\n";
		if ($orient eq '-') {
			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";
			print OUT 'mtext("'.$subtitle.'")'."\n";
                }

22:95062840f80f

my $spec = sprintf("%.1f",($ontarget / $totalmapped)*100);
# get min/max/median/average coverage => boxplot in R
open OUT, ">$wd/Rout/boxplot.R";
print OUT 'coverage <- read.table("../Targets.Global.Coverage",as.is=TRUE,sep="\t",header=FALSE)'."\n";
print OUT 'coverage <- coverage[,'.$covcol.']'."\n";
print OUT 'png(file="../Plots/CoverageBoxPlot.png", bg="white", width=240, height=480,type=c("cairo"))'."\n";
print OUT 'boxplot(coverage,range=1.5,main="Target Region Coverage")'."\n";
print OUT 'graphics.off()'."\n";
close OUT;
system("cd $wd/Rout && Rscript boxplot.R");

print OUT '    i2 <- min(which(values$cov > mincov.x))'."\n";
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";
print OUT '  }'."\n";
print OUT '}'."\n";
# open output image and create plot
print OUT 'png(file="../Plots/CoverageNtPlot.png", bg="white", width=540, height=480,type=c("cairo"))'."\n";
print OUT 'par(xaxs="i",yaxs="i")'."\n";
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";
print OUT 'lines(values$cov,values$count)'."\n";
print OUT 'if (mincov.x <= 1) {'."\n";
print OUT '  lines(c(mincov.x,mincov.x),c(0,mincov.y),lty=2,col="darkgreen")'."\n";

			else {
				$scale = 1;
			}
			my $width = 480 * $scale;
			my $height = 240 * $scale;
			print OUT 'png(file="../Plots/Coverage_'.$currgroup.'.png", bg="white", width='.$width.', height='.$height.',type=c("cairo"))'."\n";
			print OUT 'ylim = c(0,max(max(log10(coverage),log10('.($thresh+20).'))))'."\n";
			print OUT 'mp <- barplot(log10(coverage),col=colors,main="Exon Coverage for '.$currgroup.'",ylab="Log10(Coverage)",ylim=ylim)'."\n";
			print OUT 'text(mp, log10(coverage) + '.(0.4/$scale).',format(coverage),xpd = TRUE,srt=90)'."\n";
			print OUT 'text(mp,par("usr")[3]-0.05,labels=entries,srt=45,adj=1,xpd=TRUE)'."\n";
			print OUT 'abline(h=log10('.$thresh.'),lwd=4,col=rgb(255,0,0,100,maxColorValue=255))'."\n";

		else {
			$scale = 1;
		}
		my $width = 480 * $scale;
		my $height = 240 * $scale;
		print OUT 'png(file="../Plots/Coverage_'.$currgroup.'.png", bg="white", width='.$width.', height='.$height.',type=c("cairo"))'."\n";
		print OUT 'ylim = c(0,max(max(log10(coverage),log10('.($thresh+20).'))))'."\n";
		print OUT 'mp <- barplot(log10(coverage),col=colors,main="Exon Coverage for '.$currgroup.'",ylab="Log10(Coverage)", ylim=ylim)'."\n";
		print OUT 'text(mp, log10(coverage) + log10(2),format(coverage),xpd = TRUE,srt=90)'."\n";
		print OUT 'text(mp,par("usr")[3]-0.1,labels=entries,srt=45,adj=1,xpd=TRUE)'."\n";
		print OUT 'abline(h=log10('.$thresh.'),lwd=4,col=rgb(255,0,0,100,maxColorValue=255))'."\n";

		my $height = 240 ;
		my $exonstr = $exon;
		$exonstr =~ s/\s/_/g;
		$exon =~ s/_/ /g;
		$exon =~ s/\|/ /g;
		print OUT 'png(file="../Plots/Coverage_'.$exonstr.'.png", bg="white", width='.$width.', height='.$height.',type=c("cairo"))'."\n";
		print OUT 'ylim = c(0,log10(max(max(coverage),'.($thresh+10).')))'."\n";
		if ($orient eq '-') {
			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";
			print OUT 'mtext("'.$subtitle.'")'."\n";
		}

		my $height = 240 ;
		my $exonstr = $exon;
		$exonstr =~ s/\s/_/g;
		$exon =~ s/_/ /g;
		$exon =~ s/\|/ /g;
		print OUT 'png(file="../Plots/Coverage_'.$exonstr.'.png", bg="white", width='.$width.', height='.$height.',type=c("cairo"))'."\n";
		print OUT 'ylim = c(0,log10(max(max(coverage),'.($thresh+10).')))'."\n";
		if ($orient eq '-') {
			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";
			print OUT 'mtext("'.$subtitle.'")'."\n";
                }