htswanalysis/scripts/MACS_2_broadPeak.pm

   1 #!/usr/bin/perl -w
   2
   3 use strict;
   4 use warnings;
   5
   6
   7 ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   8 # Rami: May-19-2009
   9 # Transform MACS peakcalling output to broadPeak format (see doc at http://genomewiki.ucsc.edu/EncodeDCC/index.php/File_Formats)
  10 # Input can be BED file like: chr11   35886059        35886527        3228.29  (chr number | start pos  | end pos | score)
  11 ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  12
  13 #### Find Max score
  14 open (IFILE, "< $ARGV[0]") or die "Can't open INPUT file $ARGV[0]";
  15 ## vars
  16
  17 my $delimit = '[\t\s]+';
  18 my $name = '';
  19 my $signalValue = 0;
  20 my $end = 0;
  21 my $logef = 0;
  22 my $INTrec = '';
  23 my $max=-99999; # the max value
  24 my $i = 0;
  25 while(<IFILE>)
  26 {
  27   chomp;
  28   next if($_ =~ /^#/ || $_ =~ /^$/);
  29   $i++;
  30   # print "\nRead record $i: $_";
  31   if(!defined($_)) { $i++; print "\nRecord $i not defined."; next; }
  32   $INTrec='';
  33   my($chr,$start,$end,$signalValue) = split(/$delimit/,$_);
  34   if($chr =~ /^chr/)
  35   {
  36     if($max<$signalValue){$max=$signalValue;}
  37   }
  38 }
  39 ##print "\n ==== Max signalValue = $max ===\n";
  40 close IFILE;
  41 #################################################
  42
  43 ##### Generate 0-1000 score.
  44 my $BPFNAME = $ARGV[0].'.broadPeak';
  45 open (IFILE, "< $ARGV[0]") or die "Can't open INPUT file $ARGV[0]";
  46 open (BPFILE, "> $BPFNAME") or die "Can't open BROADPEAK output file $BPFNAME";
  47 my $rec = '';
  48 my $bigscore = 0; # 0-1000 score
  49 $i = 0;
  50 my $p = 0;
  51 while(<IFILE>)
  52 {
  53     $i++;
  54     if(!defined($_)) { $i++; print "\nRecord $i not defined."; next; }
  55     chomp;
  56     $rec = '';
  57
  58     my($chr,$start,$end,$signalValue) = split(/$delimit/,$_);
  59     if($chr =~ /^chr/)
  60     {
  61       $p++;
  62       $name = "peak$p";
  63       $bigscore = int( ($signalValue/$max)*1000 );
  64       $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');
  65       print BPFILE "$rec\n";
  66     }
  67     else
  68     {
  69       print "\n!! Invalid line ($i) skipped: $_\n";
  70     }
  71 }
  72
  73 close BPFILE;
  74 close IFILE;
  75
  76 print "\n==== F I N I S H E D. Created BROADPEAK file ($BPFNAME) from MACS BED file=====\n";
  77 print "\nMax signalValue = $max (corresponding to 1000)";
  78 print "\nTotal processed records = $i";
  79 print "\nNumber of valid peak records = $p";
  80 print "\n";
  81 exit;
  82
  83
  84 #################
  85 #Taken from http://encodewiki.ucsc.edu/EncodeDCC/index.php/File_Formats#broadPeak:_Broad_Peaks_.28or_Regions.29_Format   ON JUN/26/09
  86 #broadPeak: Broad Peaks (or Regions) Format
  87 #This format is used to provide called regions of signal enrichment based on pooled, normalized (interpreted) data. It is a BED 6+3 format.
  88 #    field        type      description
  89 #    ----         ---       ----------
  90 #    chrom        string    Name of the chromosome.
  91 #    chromStart   int       Starting position on the chromosome.
  92 #    chromEnd     int       Ending position on the chromosome.
  93 #    name         string    Name given to a region (preferably unique). Use `.` if no name is assigned.
  94 #    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.
  95 #    strand       char      +/- to denote strand or orientation (whenever applicable). Use '.' if no orientation is assigned.
  96 #    signalValue  float     Measurement of overall (usually, average) enrichment for the region.
  97 #    pValue       float     Measurement of statistical signficance (-log10). Use -1 if no pValue is assigned.
  98 #    qValue       float     Measurement of statistical significance using false discovery rate. Use -1 if no qValue is assigned.
  99