Mercurial > repos > geert-vandeweyer > dc_genotyper

comparison DC_Genotyper.pl @ 19:8938f339ed37 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added output of pdf report with distributions
author geert-vandeweyer
date Mon, 29 Sep 2014 05:19:48 -0400
parents 3ba482a2dd0e
children 8262299f8f3c
comparison
equal deleted inserted replaced
18:93f4d7524823 19:8938f339ed37
296 print OUT "dists <- read.table(file='".$opts{'a'}."', header=TRUE, as.is=TRUE)\n"; 296 print OUT "dists <- read.table(file='".$opts{'a'}."', header=TRUE, as.is=TRUE)\n";
297 print OUT "pdf(file='".$opts{'d'}."',paper='a4',onefile=TRUE)\n"; 297 print OUT "pdf(file='".$opts{'d'}."',paper='a4',onefile=TRUE)\n";
298 print OUT "par(mfrow=c(2, 2))\n"; 298 print OUT "par(mfrow=c(2, 2))\n";
299 print OUT "for (i in 1:nrow(dists)) {\n"; 299 print OUT "for (i in 1:nrow(dists)) {\n";
300 print OUT " if (dists[i,'avg'] > 0.5) {\n"; 300 print OUT " if (dists[i,'avg'] > 0.5) {\n";
301 print OUT " x <- seq(0.85,1,length=1000))\n"; 301 print OUT " x <- seq(0.85,1,length=1000)\n";
302 print OUT " y <- dnorm(x,mean=dists[i,'avg'],sd=dists[i,'sd']))\n"; 302 print OUT " y <- dnorm(x,mean=dists[i,'avg'],sd=dists[i,'sd'])\n";
303 print OUT " plot(x,y,main=paste('Distribution for allele \"',dists[i,'allele'],'\"',sep=''),xlab='Allelic Ratio',type='l',lwd=1))\n"; 303 print OUT " plot(x,y,main=paste('Distribution for allele \"',dists[i,'allele'],'\"',sep=''),xlab='Allelic Ratio',type='l',lwd=1)\n";
304 print OUT " abline(v=(dists[i,'avg']-3*dists[i,'sd']),col='red'))\n"; 304 print OUT " abline(v=(dists[i,'avg']-3*dists[i,'sd']),col='red')\n";
305 print OUT " text(0.855,max(y-0.5),paste(c('avg: ',round(dists[i,'avg'],digits=5),'\\nsd: ',round(dists[i,'sd'],digits=5),'\\nN: ',dists[i,'N']),sep=' ',collapse=''),adj=c(0,1)))\n"; 305 print OUT " text(0.855,max(y-0.5),paste(c('avg: ',round(dists[i,'avg'],digits=5),'\\nsd: ',round(dists[i,'sd'],digits=5),'\\nN: ',dists[i,'N']),sep=' ',collapse=''),adj=c(0,1))\n";
306 print OUT " } else {)\n"; 306 print OUT " } else {\n";
307 print OUT " x <- seq(0,0.15,length=1000))\n"; 307 print OUT " x <- seq(0,0.15,length=1000)\n";
308 print OUT " y <- dnorm(x,dists[i,'avg'],sd=dists[i,'sd']))\n"; 308 print OUT " y <- dnorm(x,dists[i,'avg'],sd=dists[i,'sd'])\n";
309 print OUT " plot(x,y,main=paste('Distribution for \"',dists[i,'allele'],'\" variation',sep=''),xlab='Allelic Ratio',type='l',lwd=1))\n"; 309 print OUT " plot(x,y,main=paste('Distribution for \"',dists[i,'allele'],'\" variation',sep=''),xlab='Allelic Ratio',type='l',lwd=1)\n";
310 print OUT " abline(v=(dists[i,'avg']+3*dists[i,'sd']),col='red'))\n"; 310 print OUT " abline(v=(dists[i,'avg']+3*dists[i,'sd']),col='red')\n";
311 print OUT " text(0.1,max(y-0.5),paste(c('avg: ',round(dists[i,'avg'],digits=5),'\\nsd: ',round(dists[i,'sd'],digits=5),'\\nN: ',dists[i,'N']),sep=' ',collapse=''),adj=c(0,1)))\n"; 311 print OUT " text(0.1,max(y-0.5),paste(c('avg: ',round(dists[i,'avg'],digits=5),'\\nsd: ',round(dists[i,'sd'],digits=5),'\\nN: ',dists[i,'N']),sep=' ',collapse=''),adj=c(0,1))\n";
312 print OUT " })\n"; 312 print OUT " }\n";
313 print OUT "})\n"; 313 print OUT "}\n";
314 print OUT "dev.off())\n"; 314 print OUT "dev.off()\n";
315 close OUT; 315 close OUT;
316 system("cd $wd && Rscript MakePlots.R > /dev/null 2>&1"); 316 system("cd $wd && Rscript MakePlots.R > /dev/null 2>&1");
317 317
318 ## CALL SNPs 318 ## CALL SNPs
319 ############ 319 ############
472 for (my $pos = 0; $pos < length($seq); $pos++) { 472 for (my $pos = 0; $pos < length($seq); $pos++) {
473 $ref_alleles{($pos+$start)} = substr($seq,$pos,1); 473 $ref_alleles{($pos+$start)} = substr($seq,$pos,1);
474 } 474 }
475 ## get counts. 475 ## get counts.
476 my $target = "$chr:$start-$stop"; 476 my $target = "$chr:$start-$stop";
477 my $command = "igvtools count -w 1 --bases --query '$target' '$bam' '$wd/WIGS/$chr-$start-$stop.wig' '$igvgenome' > /dev/null 2>&1"; 477 my $command = "igvtools count -w 1 --bases --query '$target' '$wd/input.bam' '$wd/WIGS/$chr-$start-$stop.wig' '$igvgenome' > /dev/null 2>&1";
478 system($command); 478 system($command);
479 open WIG, "$wd/WIGS/$chr-$start-$stop.wig"; 479 open WIG, "$wd/WIGS/$chr-$start-$stop.wig";
480 my $h = <WIG>; 480 my $h = <WIG>;
481 $h = <WIG>; 481 $h = <WIG>;
482 $h = <WIG>; 482 $h = <WIG>;