From: Rami Rauch <rrauch@stanford.edu>
Date: Fri, 26 Jun 2009 21:43:02 +0000 (+0000)
Subject: (no commit message)
X-Git-Url: http://woldlab.caltech.edu/gitweb/?p=htsworkflow.git;a=commitdiff_plain;h=7f4f00a4c6fa8551bff347b12964dad85c1828b6

---

diff --git a/htswanalysis/scripts/MACS_2_broadPeak.pm b/htswanalysis/scripts/MACS_2_broadPeak.pm
new file mode 100755
index 0000000..e4578a2
--- /dev/null
+++ b/htswanalysis/scripts/MACS_2_broadPeak.pm
@@ -0,0 +1,99 @@
+#!/usr/bin/perl -w
+
+use strict;
+use warnings;
+
+ 
+##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
+# Rami: May-19-2009 
+# Transform MACS peakcalling output to broadPeak format (see doc at http://genomewiki.ucsc.edu/EncodeDCC/index.php/File_Formats)
+# Input can be BED file like: chr11   35886059        35886527        3228.29  (chr number | start pos  | end pos | score)
+##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+#### Find Max score
+open (IFILE, "< $ARGV[0]") or die "Can't open INPUT file $ARGV[0]";
+## vars
+
+my $delimit = '[\t\s]+';
+my $name = '';
+my $signalValue = 0;
+my $end = 0;
+my $logef = 0;
+my $INTrec = '';
+my $max=-99999; # the max value
+my $i = 0;
+while(<IFILE>)
+{
+  chomp;
+  next if($_ =~ /^#/ || $_ =~ /^$/);
+  $i++;
+  # print "\nRead record $i: $_";
+  if(!defined($_)) { $i++; print "\nRecord $i not defined."; next; }
+  $INTrec='';
+  my($chr,$start,$end,$signalValue) = split(/$delimit/,$_);  
+  if($chr =~ /^chr/)
+  {  
+    if($max<$signalValue){$max=$signalValue;}
+  }
+}
+##print "\n ==== Max signalValue = $max ===\n";
+close IFILE;
+#################################################
+
+##### Generate 0-1000 score. 
+my $BPFNAME = $ARGV[0].'.broadPeak';
+open (IFILE, "< $ARGV[0]") or die "Can't open INPUT file $ARGV[0]";
+open (BPFILE, "> $BPFNAME") or die "Can't open BROADPEAK output file $BPFNAME";
+my $rec = '';
+my $bigscore = 0; # 0-1000 score
+$i = 0;
+my $p = 0; 
+while(<IFILE>)
+{                                                                                                                                                      
+    $i++;
+    if(!defined($_)) { $i++; print "\nRecord $i not defined."; next; }
+    chomp;
+    $rec = '';
+    
+    my($chr,$start,$end,$signalValue) = split(/$delimit/,$_);
+    if($chr =~ /^chr/)
+    {
+      $p++;
+      $name = "peak$p";
+      $bigscore = int( ($signalValue/$max)*1000 );
+      $rec = sprintf("%s\t%d\t%d\t%s\t%d\t%s\t%0.3f\t%s\t%s",$chr,$start,$end,$name,$bigscore,'.',$signalValue,'-1','-1');
+      print BPFILE "$rec\n";
+    }
+    else
+    {
+      print "\n!! Invalid line ($i) skipped: $_\n";
+    }
+}
+
+close BPFILE;
+close IFILE;
+
+print "\n==== F I N I S H E D. Created BROADPEAK file ($BPFNAME) from MACS BED file=====\n";
+print "\nMax signalValue = $max (corresponding to 1000)";
+print "\nTotal processed records = $i"; 
+print "\nNumber of valid peak records = $p";
+print "\n";
+exit;
+
+
+#################
+#Taken from http://encodewiki.ucsc.edu/EncodeDCC/index.php/File_Formats#broadPeak:_Broad_Peaks_.28or_Regions.29_Format   ON JUN/26/09
+#broadPeak: Broad Peaks (or Regions) Format
+#This format is used to provide called regions of signal enrichment based on pooled, normalized (interpreted) data. It is a BED 6+3 format.
+#    field        type      description
+#    ----         ---       ----------
+#    chrom        string    Name of the chromosome.
+#    chromStart   int       Starting position on the chromosome.
+#    chromEnd     int       Ending position on the chromosome.
+#    name         string    Name given to a region (preferably unique). Use `.` if no name is assigned.
+#    score        int       Indicates how dark the peak will be displayed in the browser (1-1000). If `0`, the DCC will assign this based on signal value. Ideally average signalValue per base spread between 100-1000.
+#    strand       char      +/- to denote strand or orientation (whenever applicable). Use '.' if no orientation is assigned.
+#    signalValue  float     Measurement of overall (usually, average) enrichment for the region.
+#    pValue       float     Measurement of statistical signficance (-log10). Use -1 if no pValue is assigned.
+#    qValue       float     Measurement of statistical significance using false discovery rate. Use -1 if no qValue is assigned. 
+