Quest_2_bed conversaion added.

diff --git a/htswanalysis/scripts/QuEST_2_BED.pm b/htswanalysis/scripts/QuEST_2_BED.pm
new file mode 100755 (executable)
index 0000000..fc8d1d8
--- /dev/null
+++ b/htswanalysis/scripts/QuEST_2_BED.pm
@@ -0,0 +1,132 @@
+#!/usr/bin/perl -w
+
+use strict;
+use warnings;
+
+ 
+##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
+# Transform QuEST peakcalling .out file to BED (graph)
+# Input rec example: chr8 w: 433 locmax: 5.68905 at: 126511693 nulls: 0.145256 ef: 27.6411204454239
+# w: we don't need.
+# locmax:is the absolute height of the peak. it does not correspond to the actual tag reads, because not all tags contribute equally to the peak determination.
+# at: + 1 is the start position (to adjust for 0 base open, which means start at 1)
+# at: + 2is the end position.
+# nulls: is the control value
+# ef: is the foldof enrichment at that peak, one base
+##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+sub lnTolog2($) 
+{
+  return log(shift)/log(2); 
+}
+
+# 1. Format into Tab delimited and remove field names inside records
+
+my $track_name = $ARGV[1] || "QuEST Peak Calling";
+
+##############################################################################
+# 2. Find EF Max and nulls MIN
+my $INTFNAME = $ARGV[0].'.INT';
+## my $BEDFNAME = $ARGV[0].'.bed';
+open (IFILE, "< $ARGV[0]") or die "Can't open file $ARGV[0]";
+open (INTFILE, "> $INTFNAME") or die "Can't open INTERMIDIATE output file";
+
+## vars
+my $i = 0;
+my $delimit = '[\t\s]+';
+my $name = '';
+my $score = 0;
+my $end = 0;
+my $logef = 0;
+my $INTrec = '';
+my $max=-999; # the max ef value
+my $min_ctrl=999; # the minimum background value
+my $seqwindow = 50;
+## print INTFILE "chr\tstart\tend\tname\tef\tlog(ef)\tscore(=rand(log(ef))\n";
+my $r=0; 
+
+while(<IFILE>)
+{
+  chomp;
+  next if($_ =~ /^#/ || $_ =~ /^$/);
+ ## chr8 w: 433 locmax: 5.68905 at: 126511693 nulls: 0.145256 ef: 27.6411204454239
+  $r++;
+  $i++;
+  # print "\nRead record $i: $_";
+  if(!defined($_)) { $i++; print "\nRecord $i not defined."; next; }
+  $INTrec='';
+  my($chr,$x1,$w,$x2,$locmax,$x3,$start,$x4,$nulls,$x5,$ef) = split(/ /,$_);  
+  $start++; # to adjust for 0 base open
+  # apply arbitrary range of 50 bases around the peak
+
+  $end = $start + $seqwindow;
+  $start -= $seqwindow;  
+
+  $name = "peak$i";
+  if($ef eq "NA") {
+    $score="NA";
+  } else {
+    $score = log($ef)/log(2);
+    if($max<$score){$max=$score;}
+    if($min_ctrl>$nulls){$min_ctrl=$nulls;} 
+  }
+  
+  $INTrec = "$chr\t$start\t$end\t$name\t$locmax\t$ef\t$logef\t$score";
+  print INTFILE "$INTrec\n";
+}
+close IFILE;
+close INTFILE;
+print "\n ==== Max EF = $max ===\n";
+
+#################################################
+
+##### Generate 0-1000 score and the BED graph score 
+my $BEDFNAME = $ARGV[0].'.bed';
+my $BEDGRAPHFNAME = $ARGV[0].'.bedgraph';
+open (IFILE, "< $INTFNAME") or die "Can't open file $INTFNAME";
+open (BEDFILE, "> $BEDFNAME") or die "Can't open BED output file $BEDFNAME";
+open (BEDGRAPHFILE, "> $BEDGRAPHFNAME") or die "Can't open BED output file $BEDGRAPHFNAME";
+my $b=0;
+my $BEDrec = '';
+my $BEDGRAPHrec = '';
+my $bigscore = 0; # 0-1000 score
+my $graphscore = 0;
+print BEDGRAPHFILE "track type=bedGraph name=\"$track_name\" description=\"$track_name\" viewLimits=0:20 autoscale=off visibility=full color=200,100,0 altColor=0,100,200 priority=20\n";
+
+while(<IFILE>)
+{
+    $b++;                                                                                                                                                      
+    if(!defined($_)) { $b++; print "\nRecord $b not defined."; next; }
+    chomp;
+    $BEDrec = '';
+    my($chr,$start,$end,$name,$locmax,$ef,$logef,$score) = split(/$delimit/,$_);
+    $bigscore = $score;
+    $graphscore = $score;
+    
+    if($graphscore =~ /^NA$/) { $graphscore = log($locmax/$min_ctrl)/log(2);}
+
+    if($graphscore > $max){$max = $graphscore;}
+    $bigscore = int( ($graphscore/$max)*1000 );
+
+    $BEDrec = "$chr\t$start\t$end\t$name\t$bigscore";
+    $BEDGRAPHrec = sprintf("%s\t%d\t%d\t%0.3f",$chr,$start,$end,$graphscore); ## ,$bigscore,$score,$locmax,$ef);
+    print BEDFILE "$BEDrec\n";
+    print BEDGRAPHFILE "$BEDGRAPHrec\n";
+}
+
+
+close BEDFILE;
+close BEDGRAPHFILE;
+close IFILE;
+
+print "\n==== F I N I S H E D. Created BED files processing $i records ===============\n";
+print "\n ==== Max EF(After clac NA values) = $max ===\n";
+print "\n ==== Min Background = $min_ctrl ===\n";
+print "\n( Temp file is $INTFNAME )";
+print "\n\nOUTPUT:  BED file is $BEDFNAME";
+print "\nOUTPUT:    BEDgraph file is $BEDGRAPHFNAME\n";
+
+if($r != $b){ print "\nSOMETHING'S WRONG: different number of records processed into BED ($b), compared to INPUT ($r).\n";}
+
+exit;
+