code_perl/Maasha/UCSC/BED.pm

   1 package Maasha::UCSC::BED;
   2
   3 # Copyright (C) 2007-2008 Martin A. Hansen.
   4
   5 # This program is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU General Public License
   7 # as published by the Free Software Foundation; either version 2
   8 # of the License, or (at your option) any later version.
   9
  10 # This program is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 # GNU General Public License for more details.
  14
  15 # You should have received a copy of the GNU General Public License
  16 # along with this program; if not, write to the Free Software
  17 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18
  19 # http://www.gnu.org/copyleft/gpl.html
  20
  21
  22 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  23
  24
  25 # Routines for interaction with Browser Extensible DATA (BED) entries and files.
  26
  27 # Read about the BED format here: http://genome.ucsc.edu/FAQ/FAQformat#format1
  28
  29
  30 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  31
  32
  33 use strict;
  34
  35 use Data::Dumper;
  36 use Maasha::Common;
  37 use Maasha::Filesys;
  38 use Maasha::Calc;
  39
  40 use vars qw( @ISA @EXPORT_OK );
  41
  42 require Exporter;
  43
  44 @ISA = qw( Exporter );
  45
  46 @EXPORT_OK = qw(
  47 );
  48
  49 use constant {
  50     chrom       => 0,   # BED field names
  51     chromStart  => 1,
  52     chromEnd    => 2,
  53     name        => 3,
  54     score       => 4,
  55     strand      => 5,
  56     thickStart  => 6,
  57     thickEnd    => 7,
  58     itemRgb     => 8,
  59     blockCount  => 9,
  60     blockSizes  => 10,
  61     blockStarts => 11,
  62     CHR         => 0,    # Biopieces field names
  63     CHR_BEG     => 1,
  64     CHR_END     => 2,
  65     Q_ID        => 3,
  66     SCORE       => 4,
  67     STRAND      => 5,
  68     THICK_BEG   => 6,
  69     THICK_END   => 7,
  70     COLOR       => 8,
  71     BLOCK_COUNT => 9,
  72     BLOCK_LENS  => 10,
  73     Q_BEGS      => 11,
  74 };
  75
  76
  77 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> BED format <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  78
  79
  80 # Hash for converting BED keys to Biopieces keys.
  81
  82 my %BED2BIOPIECES = (
  83     chrom       => "CHR",
  84     chromStart  => "CHR_BEG",
  85     chromEnd    => "CHR_END",
  86     name        => "Q_ID",
  87     score       => "SCORE",
  88     strand      => "STRAND",
  89     thickStart  => "THICK_BEG",
  90     thickEnd    => "THICK_END",
  91     itemRgb     => "COLOR",
  92     blockCount  => "BLOCK_COUNT",
  93     blockSizes  => "BLOCK_LENS",
  94     blockStarts => "Q_BEGS",
  95 );
  96
  97
  98 # Hash for converting biopieces keys to BED keys.
  99
 100 my %BIOPIECES2BED = (
 101     CHR         => "chrom",
 102     CHR_BEG     => "chromStart",
 103     CHR_END     => "chromEnd",
 104     Q_ID        => "name",
 105     SCORE       => "score",
 106     STRAND      => "strand",
 107     THICK_BEG   => "thickStart",
 108     THICK_END   => "thickEnd",
 109     COLOR       => "itemRgb",
 110     BLOCK_COUNT => "blockCount",
 111     BLOCK_LENS  => "blockSizes",
 112     Q_BEGS      => "blockStarts",
 113 );
 114
 115
 116 sub bed_entry_get
 117 {
 118     # Martin A. Hansen, September 2008.
 119
 120     # Reads a BED entry given a filehandle.
 121
 122     my ( $fh,     # file handle
 123          $cols,   # columns to read               - OPTIONAL (3,4,5,6,8,9 or 12)
 124          $check,  # check integrity of BED values - OPTIONAL
 125        ) = @_;
 126
 127     # Returns a hash.
 128
 129     my ( $line, %entry );
 130
 131     $line = <$fh>;
 132
 133     $line =~ tr/\n\r//d;    # some people have carriage returns in their BED files -> Grrrr
 134
 135     return if not defined $line;
 136
 137     if ( not defined $cols ) {
 138         $cols = 1 + $line =~ tr/\t//;
 139     }
 140
 141     if ( $cols == 3 )
 142     {
 143         (
 144             $entry{ 'chrom' },
 145             $entry{ 'chromStart' },
 146             $entry{ 'chromEnd' }
 147         ) = split "\t", $line, $cols + 1;
 148     }
 149     elsif ( $cols == 4 )
 150     {
 151         (
 152             $entry{ 'chrom' },
 153             $entry{ 'chromStart' },
 154             $entry{ 'chromEnd' },
 155             $entry{ 'name' }
 156         ) = split "\t", $line, $cols + 1;
 157     }
 158     elsif ( $cols == 5 )
 159     {
 160         (
 161             $entry{ 'chrom' },
 162             $entry{ 'chromStart' },
 163             $entry{ 'chromEnd' },
 164             $entry{ 'name' },
 165             $entry{ 'score' }
 166         ) = split "\t", $line, $cols + 1;
 167     }
 168     elsif ( $cols == 6 )
 169     {
 170         (
 171             $entry{ 'chrom' },
 172             $entry{ 'chromStart' },
 173             $entry{ 'chromEnd' },
 174             $entry{ 'name' },
 175             $entry{ 'score' },
 176             $entry{ 'strand' }
 177         ) = split "\t", $line, $cols + 1;
 178     }
 179     elsif ( $cols == 8 )
 180     {
 181         (
 182             $entry{ 'chrom' },
 183             $entry{ 'chromStart' },
 184             $entry{ 'chromEnd' },
 185             $entry{ 'name' },
 186             $entry{ 'score' },
 187             $entry{ 'strand' },
 188             $entry{ 'thickStart' },
 189             $entry{ 'thickEnd' },
 190         ) = split "\t", $line, $cols + 1;
 191     }
 192     elsif ( $cols == 9 )
 193     {
 194         (
 195             $entry{ 'chrom' },
 196             $entry{ 'chromStart' },
 197             $entry{ 'chromEnd' },
 198             $entry{ 'name' },
 199             $entry{ 'score' },
 200             $entry{ 'strand' },
 201             $entry{ 'thickStart' },
 202             $entry{ 'thickEnd' },
 203             $entry{ 'itemRgb' },
 204         ) = split "\t", $line, $cols + 1;
 205     }
 206     elsif ( $cols == 12 )
 207     {
 208         (
 209             $entry{ 'chrom' },
 210             $entry{ 'chromStart' },
 211             $entry{ 'chromEnd' },
 212             $entry{ 'name' },
 213             $entry{ 'score' },
 214             $entry{ 'strand' },
 215             $entry{ 'thickStart' },
 216             $entry{ 'thickEnd' },
 217             $entry{ 'itemRgb' },
 218             $entry{ 'blockCount' },
 219             $entry{ 'blockSizes' },
 220             $entry{ 'blockStarts' }
 221         ) = split "\t", $line, $cols + 1;
 222     }
 223     else
 224     {
 225         Maasha::Common::error( qq(Bad BED entry column count: $cols) );
 226     }
 227
 228     bed_entry_check( \%entry ) if $check;
 229
 230     return wantarray ? %entry : \%entry;
 231 }
 232
 233
 234 sub bed_entry_get_list
 235 {
 236     # Martin A. Hansen, September 2008.
 237
 238     # Reads a BED entry given a filehandle.
 239
 240     my ( $fh,     # file handle
 241          $cols,   # columns to read               - OPTIONAL (3,4,5,6,8,9 or 12)
 242          $check,  # check integrity of BED values - OPTIONAL
 243        ) = @_;
 244
 245     # Returns a list.
 246
 247     my ( $line, @entry );
 248
 249     $line = <$fh>;
 250
 251     $line =~ tr/\n\r//d;    # some people have carriage returns in their BED files -> Grrrr
 252
 253     return if not defined $line;
 254
 255     if ( not defined $cols ) {
 256         $cols = 1 + $line =~ tr/\t//;
 257     }
 258
 259     @entry = split "\t", $line, $cols + 1;
 260
 261     pop @entry if scalar @entry > $cols;
 262
 263     bed_entry_check( \@entry ) if $check;
 264
 265     return wantarray ? @entry : \@entry;
 266 }
 267
 268
 269 sub bed_entry_put
 270 {
 271     # Martin A. Hansen, September 2008.
 272
 273     # Writes a BED entry array to file.
 274
 275     my ( $entry,   # hash
 276          $fh,      # file handle                   - OPTIONAL
 277          $cols,    # number of columns in BED file - OPTIONAL (3,4,5,6,8,9 or 12)
 278          $check,   # check integrity of BED values - OPTIONAL
 279        ) = @_;
 280
 281     # Returns nothing.
 282
 283     $cols = scalar keys %{ $entry } if not $cols;
 284
 285     bed_entry_check( $entry ) if $check;
 286
 287     $fh = \*STDOUT if not defined $fh;
 288
 289     if ( $cols == 3 )
 290     {
 291         print $fh join( "\t",
 292             $entry->{ 'chrom' },
 293             $entry->{ 'chromStart' },
 294             $entry->{ 'chromEnd' },
 295         ), "\n";
 296     }
 297     elsif ( $cols == 4 )
 298     {
 299         print $fh join( "\t",
 300             $entry->{ 'chrom' },
 301             $entry->{ 'chromStart' },
 302             $entry->{ 'chromEnd' },
 303             $entry->{ 'name' },
 304         ), "\n";
 305     }
 306     elsif ( $cols == 5 )
 307     {
 308         print $fh join( "\t",
 309             $entry->{ 'chrom' },
 310             $entry->{ 'chromStart' },
 311             $entry->{ 'chromEnd' },
 312             $entry->{ 'name' },
 313             $entry->{ 'score' },
 314         ), "\n";
 315     }
 316     elsif ( $cols == 6 )
 317     {
 318         print $fh join( "\t",
 319             $entry->{ 'chrom' },
 320             $entry->{ 'chromStart' },
 321             $entry->{ 'chromEnd' },
 322             $entry->{ 'name' },
 323             $entry->{ 'score' },
 324             $entry->{ 'strand' },
 325         ), "\n";
 326     }
 327     elsif ( $cols == 8 )
 328     {
 329         print $fh join( "\t",
 330             $entry->{ 'chrom' },
 331             $entry->{ 'chromStart' },
 332             $entry->{ 'chromEnd' },
 333             $entry->{ 'name' },
 334             $entry->{ 'score' },
 335             $entry->{ 'strand' },
 336             $entry->{ 'thickStart' },
 337             $entry->{ 'thickEnd' },
 338         ), "\n";
 339     }
 340     elsif ( $cols == 9 )
 341     {
 342         print $fh join( "\t",
 343             $entry->{ 'chrom' },
 344             $entry->{ 'chromStart' },
 345             $entry->{ 'chromEnd' },
 346             $entry->{ 'name' },
 347             $entry->{ 'score' },
 348             $entry->{ 'strand' },
 349             $entry->{ 'thickStart' },
 350             $entry->{ 'thickEnd' },
 351             $entry->{ 'itemRgb' },
 352         ), "\n";
 353     }
 354     elsif ( $cols == 12 )
 355     {
 356         print $fh join( "\t",
 357             $entry->{ 'chrom' },
 358             $entry->{ 'chromStart' },
 359             $entry->{ 'chromEnd' },
 360             $entry->{ 'name' },
 361             $entry->{ 'score' },
 362             $entry->{ 'strand' },
 363             $entry->{ 'thickStart' },
 364             $entry->{ 'thickEnd' },
 365             $entry->{ 'itemRgb' },
 366             $entry->{ 'blockCount' },
 367             $entry->{ 'blockSizes' },
 368             $entry->{ 'blockStarts' },
 369         ), "\n";
 370     }
 371     else
 372     {
 373         Maasha::Common::error( qq(Bad BED entry column count: $cols) );
 374     }
 375 }
 376
 377
 378 sub bed_entry_check
 379 {
 380     # Martin A. Hansen, November 2008.
 381
 382     # Checks a BED entry for integrity and
 383     # raises an error if there is a problem.
 384
 385     my ( $bed,   # hash ref
 386        ) = @_;
 387
 388     # Returns nothing.
 389
 390     my ( $cols, @block_sizes, @block_starts );
 391
 392     $cols = scalar keys %{ $bed };
 393
 394     if ( $cols < 3 ) {
 395         Maasha::Common::error( qq(Bad BED entry - must contain at least 3 columns - not $cols) );
 396     }
 397
 398     if ( $cols > 12 ) {
 399         Maasha::Common::error( qq(Bad BED entry - must contains no more than 12 columns - not $cols) );
 400     }
 401
 402     if ( $bed->{ 'chrom' } =~ /\s/ ) {
 403         Maasha::Common::error( qq(Bad BED entry - no white space allowed in chrom field: "$bed->{ 'chrom' }") );
 404     }
 405
 406     if ( $bed->{ 'chromStart' } =~ /\D/ ) {
 407         Maasha::Common::error( qq(Bad BED entry - chromStart must be a whole number - not "$bed->{ 'chromStart' }") );
 408     }
 409
 410     if ( $bed->{ 'chromEnd' } =~ /\D/ ) {
 411         Maasha::Common::error( qq(Bad BED entry - chromEnd must be a whole number - not "$bed->{ 'chromEnd' }") );
 412     }
 413
 414     if ( $bed->{ 'chromEnd' } < $bed->{ 'chromStart' } ) {
 415         Maasha::Common::error( qq(Bad BED entry - chromEnd must be greater than chromStart - not "$bed->{ 'chromStart' } > $bed->{ 'chromEnd' }") );
 416     }
 417
 418     return if $cols == 3;
 419
 420     if ( $bed->{ 'name' } =~ /\s/ ) {
 421         Maasha::Common::error( qq(Bad BED entry - no white space allowed in name field: "$bed->{ 'name' }") );
 422     }
 423
 424     return if $cols == 4;
 425
 426     if ( $bed->{ 'score' } =~ /\D/ ) {
 427         Maasha::Common::error( qq(Bad BED entry - score must be a whole number - not "$bed->{ 'score' }") );
 428     }
 429
 430     # if ( $bed->{ 'score' } < 0 or $bed->{ 'score' } > 1000 ) { # disabled - too restrictive !
 431     if ( $bed->{ 'score' } < 0 ) {
 432         Maasha::Common::error( qq(Bad BED entry - score must be between 0 and 1000 - not "$bed->{ 'score' }") );
 433     }
 434
 435     return if $cols == 5;
 436
 437     if ( $bed->{ 'strand' } ne '+' and $bed->{ 'strand' } ne '-' ) {
 438         Maasha::Common::error( qq(Bad BED entry - strand must be + or - not "$bed->{ 'strand' }") );
 439     }
 440
 441     return if $cols == 6;
 442
 443     if ( $bed->{ 'thickStart' } =~ /\D/ ) {
 444         Maasha::Common::error( qq(Bad BED entry - thickStart must be a whole number - not "$bed->{ 'thickStart' }") );
 445     }
 446
 447     if ( $bed->{ 'thickEnd' } =~ /\D/ ) {
 448         Maasha::Common::error( qq(Bad BED entry - thickEnd must be a whole number - not "$bed->{ 'thickEnd' }") );
 449     }
 450
 451     if ( $bed->{ 'thickEnd' } < $bed->{ 'thickStart' } ) {
 452         Maasha::Common::error( qq(Bad BED entry - thickEnd must be greater than thickStart - not "$bed->{ 'thickStart' } > $bed->{ 'thickEnd' }") );
 453     }
 454
 455     if ( $bed->{ 'thickStart' } < $bed->{ 'chromStart' } ) {
 456         Maasha::Common::error( qq(Bad BED entry - thickStart must be greater than chromStart - not "$bed->{ 'thickStart' } < $bed->{ 'chromStart' }") );
 457     }
 458
 459     if ( $bed->{ 'thickStart' } > $bed->{ 'chromEnd' } ) {
 460         Maasha::Common::error( qq(Bad BED entry - thickStart must be less than chromEnd - not "$bed->{ 'thickStart' } > $bed->{ 'chromEnd' }") );
 461     }
 462
 463     if ( $bed->{ 'thickEnd' } < $bed->{ 'chromStart' } ) {
 464         Maasha::Common::error( qq(Bad BED entry - thickEnd must be greater than chromStart - not "$bed->{ 'thickEnd' } < $bed->{ 'chromStart' }") );
 465     }
 466
 467     if ( $bed->{ 'thickEnd' } > $bed->{ 'chromEnd' } ) {
 468         Maasha::Common::error( qq(Bad BED entry - thickEnd must be less than chromEnd - not "$bed->{ 'thickEnd' } > $bed->{ 'chromEnd' }") );
 469     }
 470
 471     return if $cols == 8;
 472
 473     if ( $bed->{ 'itemRgb' } !~ /^(0|\d{1,3},\d{1,3},\d{1,3},?)$/ ) {
 474         Maasha::Common::error( qq(Bad BED entry - itemRgb must be 0 or in the form of 255,0,0 - not "$bed->{ 'itemRgb' }") );
 475     }
 476
 477     return if $cols == 9;
 478
 479     if ( $bed->{ 'blockCount' } =~ /\D/ ) {
 480         Maasha::Common::error( qq(Bad BED entry - blockCount must be a whole number - not "$bed->{ 'blockCount' }") );
 481     }
 482
 483     @block_sizes = split ",", $bed->{ 'blockSizes' };
 484
 485     if ( grep /\D/, @block_sizes ) {
 486         Maasha::Common::error( qq(Bad BED entry - blockSizes must be whole numbers - not "$bed->{ 'blockSizes' }") );
 487     }
 488
 489     if ( $bed->{ 'blockCount' } != scalar @block_sizes ) {
 490         Maasha::Common::error( qq(Bad BED entry - blockSizes "$bed->{ 'blockSizes' }" must match blockCount "$bed->{ 'blockCount' }") );
 491     }
 492
 493     @block_starts = split ",", $bed->{ 'blockStarts' };
 494
 495     if ( grep /\D/, @block_starts ) {
 496         Maasha::Common::error( qq(Bad BED entry - blockStarts must be whole numbers - not "$bed->{ 'blockStarts' }") );
 497     }
 498
 499     if ( $bed->{ 'blockCount' } != scalar @block_starts ) {
 500         Maasha::Common::error( qq(Bad BED entry - blockStarts "$bed->{ 'blockStarts' }" must match blockCount "$bed->{ 'blockCount' }") );
 501     }
 502
 503     if ( $bed->{ 'chromStart' } + $block_starts[ -1 ] + $block_sizes[ -1 ] != $bed->{ 'chromEnd' } ) {
 504         Maasha::Common::error( qq(Bad BED entry - chromStart + blockStarts[last] + blockSizes[last] must equal chromEnd: ) .
 505                                qq($bed->{ 'chromStart' } + $block_starts[ -1 ] + $block_sizes[ -1 ] != $bed->{ 'chromEnd' }) );
 506     }
 507 }
 508
 509
 510 sub bed_sort
 511 {
 512     # Martin A. hansen, November 2008.
 513
 514     # Sorts a given BED file according to a given
 515     # sorting scheme:
 516     # 1: chr AND chr_beg.
 517     # 2: chr AND strand AND chr_beg.
 518     # 3: chr_beg.
 519     # 4: strand AND chr_beg.
 520
 521     # Currently 'bed_sort' is used for sorting = a standalone c program
 522     # that uses qsort (unstable sort).
 523
 524     my ( $file_in,  # path to BED file.
 525          $scheme,   # sort scheme
 526          $cols,     # number of columns in BED file
 527        ) = @_;
 528
 529     # Returns nothing.
 530
 531     my ( $file_out );
 532
 533     $file_out = "$file_in.sort";
 534
 535     Maasha::Common::run( "bed_sort", "--sort $scheme --cols $cols $file_in > $file_out" );
 536
 537     if ( Maasha::Filesys::file_size( $file_in ) !=  Maasha::Filesys::file_size( $file_out ) ) {
 538         Maasha::Common::error( qq(bed_sort of file "$file_in" failed) );
 539     }
 540
 541     rename $file_out, $file_in;
 542 }
 543
 544
 545 sub bed_file_split_on_chr
 546 {
 547     # Martin A. Hansen, November 2008.
 548
 549     # Given a path to a BED file splits
 550     # this file into one file per chromosome
 551     # in a temporary directory. Returns a hash
 552     # with chromosome name as key and the corresponding
 553     # file as value.
 554
 555     my ( $file,   # path to BED file
 556          $dir,    # working directory
 557          $cols,   # number of BED columns to read - OPTIONAL
 558        ) = @_;
 559
 560     # Returns a hashref
 561
 562     my ( $fh_in, $fh_out, $entry, %fh_hash, %file_hash );
 563
 564     $fh_in = Maasha::Filesys::file_read_open( $file );
 565
 566     while ( $entry = bed_entry_get( $fh_in, $cols ) )
 567     {
 568         if ( not exists $file_hash{ $entry->[ chrom ] } )
 569         {
 570             $fh_hash{ $entry->[ chrom ] }   = Maasha::Filesys::file_write_open( "$dir/$entry->[ chrom ]" );
 571             $file_hash{ $entry->[ chrom ] } = "$dir/$entry->[ chrom ]";
 572         }
 573
 574         $fh_out = $fh_hash{ $entry->[ chrom ] };
 575
 576         Maasha::UCSC::BED::bed_entry_put( $entry, $fh_out, $cols );
 577     }
 578
 579     map { close $_ } keys %fh_hash;
 580
 581     return wantarray ? %file_hash : \%file_hash;
 582 }
 583
 584
 585 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> BIOPIECES <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 586
 587
 588 sub bed2biopiece
 589 {
 590     # Martin A. Hansen, November 2008.
 591
 592     # Converts a BED entry given as a hashref
 593     # to a Biopiece record which is returned as
 594     # a hashref.
 595
 596     # IMPORTANT! The BED_END and THICK_END positions
 597     # will be the exact position in contrary to the
 598     # UCSC scheme.
 599
 600     my ( $bed_entry,   # BED entry as hashref
 601        ) = @_;
 602
 603     # Returns a hashref
 604
 605     my ( $cols, %bp_record );
 606
 607     $cols = scalar keys %{ $bed_entry };
 608
 609     if ( not exists $bed_entry->{ 'chrom' }      and
 610          not exists $bed_entry->{ 'chromStart' } and
 611          not exists $bed_entry->{ 'chromEnd' } )
 612     {
 613         return 0;
 614     }
 615
 616     $bp_record{ 'REC_TYPE' } = 'BED';
 617     $bp_record{ 'BED_COLS' } = $cols;
 618     $bp_record{ 'CHR' }      = $bed_entry->{ 'chrom' };
 619     $bp_record{ 'CHR_BEG' }  = $bed_entry->{ 'chromStart' };
 620     $bp_record{ 'CHR_END' }  = $bed_entry->{ 'chromEnd' } - 1;
 621     $bp_record{ 'BED_LEN' }  = $bed_entry->{ 'chromEnd' } - $bed_entry->{ 'chromStart' };
 622
 623     return wantarray ? %bp_record : \%bp_record if $cols == 3;
 624
 625     $bp_record{ 'Q_ID' } = $bed_entry->{ 'name' };
 626
 627     return wantarray ? %bp_record : \%bp_record if $cols == 4;
 628
 629     $bp_record{ 'SCORE' } = $bed_entry->{ 'score' };
 630
 631     return wantarray ? %bp_record : \%bp_record if $cols == 5;
 632
 633     $bp_record{ 'STRAND' } = $bed_entry->{ 'strand' };
 634
 635     return wantarray ? %bp_record : \%bp_record if $cols == 6;
 636
 637     $bp_record{ 'THICK_BEG' }   = $bed_entry->{ 'thickStart' };
 638     $bp_record{ 'THICK_END' }   = $bed_entry->{ 'thickEnd' } - 1;
 639
 640     return wantarray ? %bp_record : \%bp_record if $cols == 8;
 641
 642     $bp_record{ 'COLOR' }       = $bed_entry->{ 'itemRgb' };
 643
 644     return wantarray ? %bp_record : \%bp_record if $cols == 9;
 645
 646     $bp_record{ 'BLOCK_COUNT' } = $bed_entry->{ 'blockCount' };
 647     $bp_record{ 'BLOCK_LENS' }  = $bed_entry->{ 'blockSizes' };
 648     $bp_record{ 'Q_BEGS' }      = $bed_entry->{ 'blockStarts' };
 649
 650     return wantarray ? %bp_record : \%bp_record;
 651 }
 652
 653
 654 sub biopiece2bed
 655 {
 656     # Martin A. Hansen, November 2008.
 657
 658     # Converts a Biopieces record given as a hashref
 659     # to a BED record which is returned as
 660     # an arrayref. As much as possible of the Biopiece
 661     # record is converted and undef is returned if
 662     # convertion failed.
 663
 664     # IMPORTANT! The BED_END and THICK_END positions
 665     # will be the inexact position used in the
 666     # UCSC scheme.
 667
 668     my ( $bp_record,   # hashref
 669          $cols,        # number of columns in BED entry - OPTIONAL (but faster)
 670        ) = @_;
 671
 672     # Returns arrayref.
 673
 674     my ( %bed_entry );
 675
 676     $cols ||= 12;   # max number of columns possible
 677
 678     $bed_entry{ 'chrom' } = $bp_record->{ 'CHR' }  ||
 679                             $bp_record->{ 'S_ID' } ||
 680                             return undef;
 681
 682     if ( exists $bp_record->{ 'CHR_BEG' } ) {
 683         $bed_entry{ 'chromStart' } = $bp_record->{ 'CHR_BEG' };
 684     } elsif ( exists $bp_record->{ 'S_BEG' } ) {
 685         $bed_entry{ 'chromStart' } = $bp_record->{ 'S_BEG' };
 686     } else {
 687         return undef;
 688     }
 689
 690     if ( exists $bp_record->{ 'CHR_END' } ) {
 691         $bed_entry{ 'chromEnd' } = $bp_record->{ 'CHR_END' };
 692     } elsif ( exists $bp_record->{ 'S_END' }) {
 693         $bed_entry{ 'chromEnd' } = $bp_record->{ 'S_END' };
 694     } else {
 695         return undef;
 696     }
 697
 698     $bed_entry{ 'chromEnd' }++;
 699
 700     return wantarray ? %bed_entry : \%bed_entry if $cols == 3;
 701
 702     $bed_entry{ 'name' } = $bp_record->{ 'Q_ID' } || return wantarray ? %bed_entry : \%bed_entry;
 703
 704     return wantarray ? %bed_entry : \%bed_entry if $cols == 4;
 705
 706     if ( exists $bp_record->{ 'SCORE' } ) {
 707         $bed_entry{ 'score' } = $bp_record->{ 'SCORE' };
 708     } else {
 709         return wantarray ? %bed_entry : \%bed_entry;
 710     }
 711
 712     return wantarray ? %bed_entry : \%bed_entry if $cols == 5;
 713
 714     if ( exists $bp_record->{ 'STRAND' } ) {
 715         $bed_entry{ 'strand' } = $bp_record->{ 'STRAND' };
 716     } else {
 717         return wantarray ? %bed_entry : \%bed_entry;
 718     }
 719
 720     return wantarray ? %bed_entry : \%bed_entry if $cols == 6;
 721
 722     if ( defined $bp_record->{ 'THICK_BEG' }   and
 723          defined $bp_record->{ 'THICK_END' } )
 724     {
 725         $bed_entry{ 'thickStart' } = $bp_record->{ 'THICK_BEG' };
 726         $bed_entry{ 'thickEnd' }   = $bp_record->{ 'THICK_END' };
 727     }
 728
 729     return wantarray ? %bed_entry : \%bed_entry if $cols == 8;
 730
 731     if ( exists $bp_record->{ 'COLOR' } ) {
 732         $bed_entry{ 'itemRgb' } = $bp_record->{ 'COLOR' };
 733     } else {
 734         return wantarray ? %bed_entry : \%bed_entry;
 735     }
 736
 737     return wantarray ? %bed_entry : \%bed_entry if $cols == 9;
 738
 739     if ( defined $bp_record->{ 'BLOCK_COUNT' } and
 740          defined $bp_record->{ 'BLOCK_LENS' }  and
 741          defined $bp_record->{ 'Q_BEGS' } )
 742     {
 743         $bed_entry{ 'blockCount' }  = $bp_record->{ 'BLOCK_COUNT' };
 744         $bed_entry{ 'blockSizes' }  = $bp_record->{ 'BLOCK_LENS' };
 745         $bed_entry{ 'blockStarts' } = $bp_record->{ 'Q_BEGS' };
 746         $bed_entry{ 'thickEnd' }++;
 747     }
 748
 749     return wantarray ? %bed_entry : \%bed_entry;
 750 }
 751
 752
 753 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> TAG CONTIGS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 754
 755
 756 sub tag_contigs_assemble
 757 {
 758     # Martin A. Hansen, November 2008.
 759
 760     # Returns a path to a BED file with tab contigs.
 761
 762     my ( $bed_file,   # path to BED file
 763          $chr,        # chromosome
 764          $strand,     # strand
 765          $dir,        # working directory
 766        ) = @_;
 767
 768     # Returns a string
 769
 770     my ( $fh_in, $fh_out, $array, $new_file, $pos, $id, $beg, $end, $score, $entry );
 771
 772     $fh_in = Maasha::Filesys::file_read_open( $bed_file );
 773
 774     $array = tag_contigs_assemble_array( $fh_in );
 775
 776     close $fh_in;
 777
 778     $new_file = "$bed_file.tag_contigs";
 779
 780     $fh_out = Maasha::Filesys::file_write_open( $new_file );
 781
 782     $pos = 0;
 783     $id  = 0;
 784
 785     while ( ( $beg, $end, $score ) = tag_contigs_scan( $array, $pos ) and $beg )
 786     {
 787         $entry->[ chrom ]      = $chr;
 788         $entry->[ chromStart ] = $beg;
 789         $entry->[ chromEnd ]   = $end;
 790         $entry->[ name ]       = sprintf( "TC%06d", $id );
 791         $entry->[ score ]      = $score;
 792         $entry->[ strand ]     = $strand;
 793
 794         bed_entry_put( $entry, $fh_out );
 795
 796         $pos = $end;
 797         $id++;
 798     }
 799
 800     close $fh_out;
 801
 802     return $new_file;
 803 }
 804
 805
 806 sub tag_contigs_assemble_array
 807 {
 808     # Martin A. Hansen, November 2008.
 809
 810     # Given a BED file with entries from only one
 811     # chromosome assembles tag contigs from these
 812     # ignoring strand information. Only tags with
 813     # a score higher than the clone count over
 814     # genomic loci (the score field) is included
 815     # in the tag contigs.
 816
 817     #       -----------              tags
 818     #          -------------
 819     #               ----------
 820     #                     ----------
 821     #       ========================  tag contig
 822
 823
 824     my ( $fh,   # file handle to BED file
 825        ) = @_;
 826
 827     # Returns an arrayref.
 828
 829     my ( $entry, $clones, $score, @array );
 830
 831     while ( $entry = bed_entry_get( $fh ) )
 832     {
 833         if ( $entry->[ name ] =~ /_(\d+)$/ )
 834         {
 835             $clones = $1;
 836
 837             if ( $entry->[ score ] )
 838             {
 839                 $score = int( $clones / $entry->[ score ] );
 840
 841                 map { $array[ $_ ] += $score } $entry->[ chromStart ] .. $entry->[ chromEnd ] - 1 if $score >= 1;
 842             }
 843         }
 844     }
 845
 846     return wantarray ? @array : \@array;
 847 }
 848
 849
 850 sub tag_contigs_scan
 851 {
 852     # Martin A. Hansen, November 2008.
 853
 854     # Scans an array with tag contigs and locates
 855     # the next contig from a given position. The
 856     # score of the tag contig is determined as the
 857     # maximum value of the tag contig. If a tag contig
 858     # is found a triple is returned with beg, end and score
 859     # otherwise an empty list is returned.
 860
 861     my ( $array,   # array to scan
 862          $beg,     # position to start scanning from
 863        ) = @_;
 864
 865     # Returns an arrayref.
 866
 867     my ( $end, $score );
 868
 869     $score = 0;
 870
 871     while ( $beg < scalar @{ $array } and not $array->[ $beg ] ) { $beg++ }
 872
 873     $end = $beg;
 874
 875     while ( $array->[ $end ] )
 876     {
 877         $score = Maasha::Calc::max( $score, $array->[ $end ] );
 878
 879         $end++
 880     }
 881
 882     if ( $score > 0 ) {
 883         return wantarray ? ( $beg, $end, $score ) : [ $beg, $end, $score ];
 884     } else {
 885         return wantarray ? () : [];
 886     }
 887 }
 888
 889
 890 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 891
 892
 893 1;
 894
 895
 896 __END__