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 my $CHECK_ALL = 1;   # Global flag indicating that BED input and output is checked thoroughly.
  78
  79
  80 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> BED format <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  81
  82
  83 # Hash for converting BED keys to Biopieces keys.
  84
  85 my %BED2BIOPIECES = (
  86     chrom       => "CHR",
  87     chromStart  => "CHR_BEG",
  88     chromEnd    => "CHR_END",
  89     name        => "Q_ID",
  90     score       => "SCORE",
  91     strand      => "STRAND",
  92     thickStart  => "THICK_BEG",
  93     thickEnd    => "THICK_END",
  94     itemRgb     => "COLOR",
  95     blockCount  => "BLOCK_COUNT",
  96     blockSizes  => "BLOCK_LENS",
  97     blockStarts => "Q_BEGS",
  98 );
  99
 100
 101 # Hash for converting biopieces keys to BED keys.
 102
 103 my %BIOPIECES2BED = (
 104     CHR         => "chrom",
 105     CHR_BEG     => "chromStart",
 106     CHR_END     => "chromEnd",
 107     Q_ID        => "name",
 108     SCORE       => "score",
 109     STRAND      => "strand",
 110     THICK_BEG   => "thickStart",
 111     THICK_END   => "thickEnd",
 112     COLOR       => "itemRgb",
 113     BLOCK_COUNT => "blockCount",
 114     BLOCK_LENS  => "blockSizes",
 115     Q_BEGS      => "blockStarts",
 116 );
 117
 118
 119 sub bed_entry_get
 120 {
 121     # Martin A. Hansen, September 2008.
 122
 123     # Reads a BED entry given a filehandle.
 124
 125     my ( $fh,     # file handle
 126          $cols,   # columns to read - OPTIONAL (3,4,5,6,8,9 or 12)
 127        ) = @_;
 128
 129     # Returns a list.
 130
 131     my ( $line, @entry );
 132
 133     $line = <$fh>;
 134
 135     $line =~ tr/\n\r//d;    # some people have carriage returns in their BED files -> Grrrr
 136
 137     return if not defined $line;
 138
 139     if ( not defined $cols ) {
 140         $cols = 1 + $line =~ tr/\t//;
 141     }
 142
 143     @entry = split "\t", $line, $cols + 1;
 144
 145     pop @entry if scalar @entry > $cols;
 146
 147     bed_entry_check( \@entry ) if $CHECK_ALL;
 148
 149     return wantarray ? @entry : \@entry;
 150 }
 151
 152
 153 sub bed_entry_put
 154 {
 155     # Martin A. Hansen, September 2008.
 156
 157     # Writes a BED entry array to file.
 158
 159     my ( $entry,   # list
 160          $fh,      # file handle                   - OPTIONAL
 161          $cols,    # number of columns in BED file - OPTIONAL (3,4,5,6,8,9 or 12)
 162        ) = @_;
 163
 164     # Returns nothing.
 165
 166     if ( $cols and $cols < scalar @{ $entry } ) {
 167         @{ $entry } = @{ $entry }[ 0 .. $cols - 1 ];
 168     }
 169
 170     bed_entry_check( $entry ) if $CHECK_ALL;
 171
 172     $fh = \*STDOUT if not defined $fh;
 173
 174     print $fh join( "\t", @{ $entry } ), "\n";
 175 }
 176
 177
 178
 179 sub bed_entry_check
 180 {
 181     # Martin A. Hansen, November 2008.
 182
 183     # Checks a BED entry for integrity and
 184     # raises an error if there is a problem.
 185
 186     my ( $bed,   # array ref
 187        ) = @_;
 188
 189     # Returns nothing.
 190
 191     my ( $cols, @block_sizes, @block_starts );
 192
 193     $cols = scalar @{ $bed };
 194
 195     if ( $cols < 3 ) {
 196         Maasha::Common::error( qq(Bad BED entry - must contain at least 3 columns - not $cols) );
 197     }
 198
 199     if ( $cols > 12 ) {
 200         Maasha::Common::error( qq(Bad BED entry - must contains no more than 12 columns - not $cols) );
 201     }
 202
 203     if ( $bed->[ chrom ] =~ /\s/ ) {
 204         Maasha::Common::error( qq(Bad BED entry - no white space allowed in chrom field: "$bed->[ chrom ]") );
 205     }
 206
 207     if ( $bed->[ chromStart ] =~ /\D/ ) {
 208         Maasha::Common::error( qq(Bad BED entry - chromStart must be a whole number - not "$bed->[ chromStart ]") );
 209     }
 210
 211     if ( $bed->[ chromEnd ] =~ /\D/ ) {
 212         Maasha::Common::error( qq(Bad BED entry - chromEnd must be a whole number - not "$bed->[ chromEnd ]") );
 213     }
 214
 215     if ( $bed->[ chromEnd ] < $bed->[ chromStart ] ) {
 216         Maasha::Common::error( qq(Bad BED entry - chromEnd must be greater than chromStart - not "$bed->[ chromStart ] > $bed->[ chromEnd ]") );
 217     }
 218
 219     return if @{ $bed } == 3;
 220
 221     if ( $bed->[ name ] =~ /\s/ ) {
 222         Maasha::Common::error( qq(Bad BED entry - no white space allowed in name field: "$bed->[ name ]") );
 223     }
 224
 225     return if @{ $bed } == 4;
 226
 227     if ( $bed->[ score ] =~ /\D/ ) {
 228         Maasha::Common::error( qq(Bad BED entry - score must be a whole number - not "$bed->[ score ]") );
 229     }
 230
 231     if ( $bed->[ score ] < 0 or $bed->[ score ] > 1000 ) {
 232         Maasha::Common::error( qq(Bad BED entry - score must be between 0 and 1000 - not "$bed->[ score ]") );
 233     }
 234
 235     return if @{ $bed } == 5;
 236
 237     if ( $bed->[ strand ] ne '+' and $bed->[ strand ] ne '-' ) {
 238         Maasha::Common::error( qq(Bad BED entry - strand must be + or - not "$bed->[ strand ]") );
 239     }
 240
 241     return if @{ $bed } == 6;
 242
 243     if ( $bed->[ thickStart ] =~ /\D/ ) {
 244         Maasha::Common::error( qq(Bad BED entry - thickStart must be a whole number - not "$bed->[ thickStart ]") );
 245     }
 246
 247     if ( $bed->[ thickEnd ] =~ /\D/ ) {
 248         Maasha::Common::error( qq(Bad BED entry - thickEnd must be a whole number - not "$bed->[ thickEnd ]") );
 249     }
 250
 251     if ( $bed->[ thickEnd ] < $bed->[ thickStart ] ) {
 252         Maasha::Common::error( qq(Bad BED entry - thickEnd must be greater than thickStart - not "$bed->[ thickStart ] > $bed->[ thickEnd ]") );
 253     }
 254
 255     if ( $bed->[ thickStart ] < $bed->[ chromStart ] ) {
 256         Maasha::Common::error( qq(Bad BED entry - thickStart must be greater than chromStart - not "$bed->[ thickStart ] < $bed->[ chromStart ]") );
 257     }
 258
 259     if ( $bed->[ thickStart ] > $bed->[ chromEnd ] ) {
 260         Maasha::Common::error( qq(Bad BED entry - thickStart must be less than chromEnd - not "$bed->[ thickStart ] > $bed->[ chromEnd ]") );
 261     }
 262
 263     if ( $bed->[ thickEnd ] < $bed->[ chromStart ] ) {
 264         Maasha::Common::error( qq(Bad BED entry - thickEnd must be greater than chromStart - not "$bed->[ thickEnd ] < $bed->[ chromStart ]") );
 265     }
 266
 267     if ( $bed->[ thickEnd ] > $bed->[ chromEnd ] ) {
 268         Maasha::Common::error( qq(Bad BED entry - thickEnd must be less than chromEnd - not "$bed->[ thickEnd ] > $bed->[ chromEnd ]") );
 269     }
 270
 271     return if @{ $bed } == 8;
 272
 273     if ( $bed->[ itemRgb ] !~ /^(0|\d{1,3},\d{1,3},\d{1,3},?)$/ ) {
 274         Maasha::Common::error( qq(Bad BED entry - itemRgb must be 0 or in the form of 255,0,0 - not "$bed->[ itemRgb ]") );
 275     }
 276
 277     return if @{ $bed } == 9;
 278
 279     if ( $bed->[ blockCount ] =~ /\D/ ) {
 280         Maasha::Common::error( qq(Bad BED entry - blockCount must be a whole number - not "$bed->[ blockCount ]") );
 281     }
 282
 283     @block_sizes = split ",", $bed->[ blockSizes ];
 284
 285     if ( grep /\D/, @block_sizes ) {
 286         Maasha::Common::error( qq(Bad BED entry - blockSizes must be whole numbers - not "$bed->[ blockSizes ]") );
 287     }
 288
 289     if ( $bed->[ blockCount ] != scalar @block_sizes ) {
 290         Maasha::Common::error( qq(Bad BED entry - blockSizes "$bed->[ blockSizes ]" must match blockCount "$bed->[ blockCount ]") );
 291     }
 292
 293     @block_starts = split ",", $bed->[ blockStarts ];
 294
 295     if ( grep /\D/, @block_starts ) {
 296         Maasha::Common::error( qq(Bad BED entry - blockStarts must be whole numbers - not "$bed->[ blockStarts ]") );
 297     }
 298
 299     if ( $bed->[ blockCount ] != scalar @block_starts ) {
 300         Maasha::Common::error( qq(Bad BED entry - blockStarts "$bed->[ blockStarts ]" must match blockCount "$bed->[ blockCount ]") );
 301     }
 302
 303     if ( $bed->[ chromStart ] + $block_starts[ -1 ] + $block_sizes[ -1 ] != $bed->[ chromEnd ] ) {
 304         Maasha::Common::error( qq(Bad BED entry - chromStart + blockStarts[last] + blockSizes[last] must equal chromEnd: ) .
 305                                qq($bed->[ chromStart ] + $block_starts[ -1 ] + $block_sizes[ -1 ] != $bed->[ chromEnd ]) );
 306     }
 307 }
 308
 309
 310 sub bed_sort
 311 {
 312     # Martin A. hansen, November 2008.
 313
 314     # Sorts a given BED file according to a given
 315     # sorting scheme:
 316     # 1: chr AND chr_beg.
 317     # 2: chr AND strand AND chr_beg.
 318     # 3: chr_beg.
 319     # 4: strand AND chr_beg.
 320
 321     # Currently 'bed_sort' is used for sorting = a standalone c program
 322     # that uses qsort (unstable sort).
 323
 324     my ( $file_in,  # path to BED file.
 325          $scheme,   # sort scheme
 326          $cols,     # number of columns in BED file
 327        ) = @_;
 328
 329     # Returns nothing.
 330
 331     my ( $file_out );
 332
 333     $file_out = "$file_in.sort";
 334
 335     Maasha::Common::run( "bed_sort", "--sort $scheme --cols $cols $file_in > $file_out" );
 336
 337     if ( Maasha::Filesys::file_size( $file_in ) !=  Maasha::Filesys::file_size( $file_out ) ) {
 338         Maasha::Common::error( qq(bed_sort of file "$file_in" failed) );
 339     }
 340
 341     rename $file_out, $file_in;
 342 }
 343
 344
 345 sub bed_file_split_on_chr
 346 {
 347     # Martin A. Hansen, November 2008.
 348
 349     # Given a path to a BED file splits
 350     # this file into one file per chromosome
 351     # in a temporary directory. Returns a hash
 352     # with chromosome name as key and the corresponding
 353     # file as value.
 354
 355     my ( $file,   # path to BED file
 356          $dir,    # working directory
 357          $cols,   # number of BED columns to read - OPTIONAL
 358        ) = @_;
 359
 360     # Returns a hashref
 361
 362     my ( $fh_in, $fh_out, $entry, %fh_hash, %file_hash );
 363
 364     $fh_in = Maasha::Filesys::file_read_open( $file );
 365
 366     while ( $entry = bed_entry_get( $fh_in, $cols ) )
 367     {
 368         if ( not exists $file_hash{ $entry->[ chrom ] } )
 369         {
 370             $fh_hash{ $entry->[ chrom ] }   = Maasha::Filesys::file_write_open( "$dir/$entry->[ chrom ]" );
 371             $file_hash{ $entry->[ chrom ] } = "$dir/$entry->[ chrom ]";
 372         }
 373
 374         $fh_out = $fh_hash{ $entry->[ chrom ] };
 375
 376         Maasha::UCSC::BED::bed_entry_put( $entry, $fh_out, $cols );
 377     }
 378
 379     map { close $_ } keys %fh_hash;
 380
 381     return wantarray ? %file_hash : \%file_hash;
 382 }
 383
 384
 385 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> BIOPIECES <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 386
 387
 388 sub bed2biopiece
 389 {
 390     # Martin A. Hansen, November 2008.
 391
 392     # Converts a BED entry given as an arrayref
 393     # to a Biopiece record which is returned as
 394     # a hashref.
 395
 396     # IMPORTANT! The BED_END and THICK_END positions
 397     # will be the exact position in contrary to the
 398     # UCSC scheme.
 399
 400     my ( $bed_entry,   # BED entry as arrayref
 401        ) = @_;
 402
 403     # Returns a hashref
 404
 405     my ( $cols, %bp_record );
 406
 407     $cols = scalar @{ $bed_entry };
 408
 409     if ( not defined $bed_entry->[ chrom ] and
 410          not defined $bed_entry->[ chromStart ] and
 411          not defined $bed_entry->[ chromEnd ] )
 412     {
 413         return 0;
 414     }
 415
 416     $bp_record{ "REC_TYPE" } = "BED";
 417     $bp_record{ "BED_COLS" } = $cols;
 418     $bp_record{ "CHR" }      = $bed_entry->[ chrom ];
 419     $bp_record{ "CHR_BEG" }  = $bed_entry->[ chromStart ];
 420     $bp_record{ "CHR_END" }  = $bed_entry->[ chromEnd ] - 1;
 421     $bp_record{ "BED_LEN" }  = $bed_entry->[ chromEnd ] - $bed_entry->[ chromEnd ];
 422
 423     return wantarray ? %bp_record : \%bp_record if $cols == 3;
 424
 425     $bp_record{ "Q_ID" } = $bed_entry->[ name ];
 426
 427     return wantarray ? %bp_record : \%bp_record if $cols == 4;
 428
 429     $bp_record{ "SCORE" } = $bed_entry->[ score ];
 430
 431     return wantarray ? %bp_record : \%bp_record if $cols == 5;
 432
 433     $bp_record{ "STRAND" } = $bed_entry->[ strand ];
 434
 435     return wantarray ? %bp_record : \%bp_record if $cols == 6;
 436
 437     $bp_record{ "THICK_BEG" }   = $bed_entry->[ thickStart ];
 438     $bp_record{ "THICK_END" }   = $bed_entry->[ thickEnd ] - 1;
 439
 440     return wantarray ? %bp_record : \%bp_record if $cols == 8;
 441
 442     $bp_record{ "COLOR" }       = $bed_entry->[ itemRgb ];
 443
 444     return wantarray ? %bp_record : \%bp_record if $cols == 9;
 445
 446     $bp_record{ "BLOCK_COUNT" } = $bed_entry->[ blockCount ];
 447     $bp_record{ "BLOCK_LENS" }  = $bed_entry->[ blockSizes ];
 448     $bp_record{ "Q_BEGS" }      = $bed_entry->[ blockStarts ];
 449
 450     return wantarray ? %bp_record : \%bp_record;
 451 }
 452
 453
 454 sub biopiece2bed
 455 {
 456     # Martin A. Hansen, November 2008.
 457
 458     # Converts a Biopieces record given as a hashref
 459     # to a BED record which is returned as
 460     # an arrayref. As much as possible of the Biopiece
 461     # record is converted and undef is returned if
 462     # convertion failed.
 463
 464     # IMPORTANT! The BED_END and THICK_END positions
 465     # will be the inexact position used in the
 466     # UCSC scheme.
 467
 468     my ( $bp_record,   # hashref
 469          $cols,        # number of columns in BED entry - OPTIONAL (but faster)
 470        ) = @_;
 471
 472     # Returns arrayref.
 473
 474     my ( @bed_entry );
 475
 476     $cols ||= 12;   # max number of columns possible
 477
 478     $bed_entry[ chrom ]      = $bp_record->{ "CHR" }     ||
 479                                $bp_record->{ "S_ID" }    ||
 480                                return undef;
 481
 482     $bed_entry[ chromStart ] = $bp_record->{ "CHR_BEG" } ||
 483                                $bp_record->{ "S_BEG" }   ||
 484                                return undef;
 485
 486     $bed_entry[ chromEnd ]   = $bp_record->{ "CHR_END" } ||
 487                                $bp_record->{ "S_END" }   ||
 488                                return undef;
 489
 490     $bed_entry[ chromEnd ]++;
 491
 492     return wantarray ? @bed_entry : \@bed_entry if $cols == 3;
 493
 494     $bed_entry[ name ] = $bp_record->{ "Q_ID" } || return wantarray ? @bed_entry : \@bed_entry;
 495
 496     return wantarray ? @bed_entry : \@bed_entry if $cols == 4;
 497
 498     if ( exists $bp_record->{ "SCORE" } ) {
 499         $bed_entry[ score ] = $bp_record->{ "SCORE" };
 500     } else {
 501         return wantarray ? @bed_entry : \@bed_entry;
 502     }
 503
 504     return wantarray ? @bed_entry : \@bed_entry if $cols == 5;
 505
 506     if ( exists $bp_record->{ "STRAND" } ) {
 507         $bed_entry[ strand ] = $bp_record->{ "STRAND" };
 508     } else {
 509         return wantarray ? @bed_entry : \@bed_entry;
 510     }
 511
 512     return wantarray ? @bed_entry : \@bed_entry if $cols == 6;
 513
 514     if ( defined $bp_record->{ "THICK_BEG" }   and
 515          defined $bp_record->{ "THICK_END" } )
 516     {
 517         $bed_entry[ thickStart ] = $bp_record->{ "THICK_BEG" };
 518         $bed_entry[ thickEnd ]   = $bp_record->{ "THICK_END" };
 519     }
 520
 521     return wantarray ? @bed_entry : \@bed_entry if $cols == 8;
 522
 523     if ( exists $bp_record->{ "COLOR" } ) {
 524         $bed_entry[ itemRgb ] = $bp_record->{ "COLOR" };
 525     } else {
 526         return wantarray ? @bed_entry : \@bed_entry;
 527     }
 528
 529     return wantarray ? @bed_entry : \@bed_entry if $cols == 9;
 530
 531     if ( defined $bp_record->{ "BLOCK_COUNT" } and
 532          defined $bp_record->{ "BLOCK_LENS" }  and
 533          defined $bp_record->{ "Q_BEGS" } )
 534     {
 535         $bed_entry[ blockCount ]  = $bp_record->{ "BLOCK_COUNT" };
 536         $bed_entry[ blockSizes ]  = $bp_record->{ "BLOCK_LENS" };
 537         $bed_entry[ blockStarts ] = $bp_record->{ "Q_BEGS" };
 538         $bed_entry[ thickEnd ]++;
 539     }
 540
 541     return wantarray ? @bed_entry : \@bed_entry;
 542 }
 543
 544
 545 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> TAG CONTIGS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 546
 547
 548 sub tag_contigs_assemble
 549 {
 550     # Martin A. Hansen, November 2008.
 551
 552     # Returns a path to a BED file with tab contigs.
 553
 554     my ( $bed_file,   # path to BED file
 555          $chr,        # chromosome
 556          $strand,     # strand
 557          $dir,        # working directory
 558        ) = @_;
 559
 560     # Returns a string
 561
 562     my ( $fh_in, $fh_out, $array, $new_file, $pos, $id, $beg, $end, $score, $entry );
 563
 564     $fh_in = Maasha::Filesys::file_read_open( $bed_file );
 565
 566     $array = tag_contigs_assemble_array( $fh_in );
 567
 568     close $fh_in;
 569
 570     $new_file = "$bed_file.tag_contigs";
 571
 572     $fh_out = Maasha::Filesys::file_write_open( $new_file );
 573
 574     $pos = 0;
 575     $id  = 0;
 576
 577     while ( ( $beg, $end, $score ) = tag_contigs_scan( $array, $pos ) and $beg )
 578     {
 579         $entry->[ chrom ]      = $chr;
 580         $entry->[ chromStart ] = $beg;
 581         $entry->[ chromEnd ]   = $end;
 582         $entry->[ name ]       = sprintf( "TC%06d", $id );
 583         $entry->[ score ]      = $score;
 584         $entry->[ strand ]     = $strand;
 585
 586         bed_entry_put( $entry, $fh_out );
 587
 588         $pos = $end;
 589         $id++;
 590     }
 591
 592     close $fh_out;
 593
 594     return $new_file;
 595 }
 596
 597
 598 sub tag_contigs_assemble_array
 599 {
 600     # Martin A. Hansen, November 2008.
 601
 602     # Given a BED file with entries from only one
 603     # chromosome assembles tag contigs from these
 604     # ignoring strand information. Only tags with
 605     # a score higher than the clone count over
 606     # genomic loci (the score field) is included
 607     # in the tag contigs.
 608
 609     #       -----------              tags
 610     #          -------------
 611     #               ----------
 612     #                     ----------
 613     #       ========================  tag contig
 614
 615
 616     my ( $fh,   # file handle to BED file
 617        ) = @_;
 618
 619     # Returns an arrayref.
 620
 621     my ( $entry, $clones, $score, @array );
 622
 623     while ( $entry = bed_entry_get( $fh ) )
 624     {
 625         if ( $entry->[ name ] =~ /_(\d+)$/ )
 626         {
 627             $clones = $1;
 628
 629             if ( $entry->[ score ] )
 630             {
 631                 $score = int( $clones / $entry->[ score ] );
 632
 633                 map { $array[ $_ ] += $score } $entry->[ chromStart ] .. $entry->[ chromEnd ] - 1 if $score >= 1;
 634             }
 635         }
 636     }
 637
 638     return wantarray ? @array : \@array;
 639 }
 640
 641
 642 sub tag_contigs_scan
 643 {
 644     # Martin A. Hansen, November 2008.
 645
 646     # Scans an array with tag contigs and locates
 647     # the next contig from a given position. The
 648     # score of the tag contig is determined as the
 649     # maximum value of the tag contig. If a tag contig
 650     # is found a triple is returned with beg, end and score
 651     # otherwise an empty list is returned.
 652
 653     my ( $array,   # array to scan
 654          $beg,     # position to start scanning from
 655        ) = @_;
 656
 657     # Returns an arrayref.
 658
 659     my ( $end, $score );
 660
 661     $score = 0;
 662
 663     while ( $beg < scalar @{ $array } and not $array->[ $beg ] ) { $beg++ }
 664
 665     $end = $beg;
 666
 667     while ( $array->[ $end ] )
 668     {
 669         $score = Maasha::Calc::max( $score, $array->[ $end ] );
 670
 671         $end++
 672     }
 673
 674     if ( $score > 0 ) {
 675         return wantarray ? ( $beg, $end, $score ) : [ $beg, $end, $score ];
 676     } else {
 677         return wantarray ? () : [];
 678     }
 679 }
 680
 681
 682 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 683
 684
 685 1;
 686
 687
 688 __END__