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 list.
 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     @entry = split "\t", $line, $cols + 1;
 142
 143     pop @entry if scalar @entry > $cols;
 144
 145     bed_entry_check( \@entry ) if $check;
 146
 147     return wantarray ? @entry : \@entry;
 148 }
 149
 150
 151 sub bed_entry_put
 152 {
 153     # Martin A. Hansen, September 2008.
 154
 155     # Writes a BED entry array to file.
 156
 157     my ( $entry,   # list
 158          $fh,      # file handle                   - OPTIONAL
 159          $cols,    # number of columns in BED file - OPTIONAL (3,4,5,6,8,9 or 12)
 160          $check,   # check integrity of BED values - OPTIONAL
 161        ) = @_;
 162
 163     # Returns nothing.
 164
 165     if ( $cols and $cols < scalar @{ $entry } ) {
 166         @{ $entry } = @{ $entry }[ 0 .. $cols - 1 ];
 167     }
 168
 169     bed_entry_check( $entry ) if $check;
 170
 171     $fh = \*STDOUT if not defined $fh;
 172
 173     print $fh join( "\t", @{ $entry } ), "\n";
 174 }
 175
 176
 177
 178 sub bed_entry_check
 179 {
 180     # Martin A. Hansen, November 2008.
 181
 182     # Checks a BED entry for integrity and
 183     # raises an error if there is a problem.
 184
 185     my ( $bed,   # array ref
 186        ) = @_;
 187
 188     # Returns nothing.
 189
 190     my ( $cols, @block_sizes, @block_starts );
 191
 192     $cols = scalar @{ $bed };
 193
 194     if ( $cols < 3 ) {
 195         Maasha::Common::error( qq(Bad BED entry - must contain at least 3 columns - not $cols) );
 196     }
 197
 198     if ( $cols > 12 ) {
 199         Maasha::Common::error( qq(Bad BED entry - must contains no more than 12 columns - not $cols) );
 200     }
 201
 202     if ( $bed->[ chrom ] =~ /\s/ ) {
 203         Maasha::Common::error( qq(Bad BED entry - no white space allowed in chrom field: "$bed->[ chrom ]") );
 204     }
 205
 206     if ( $bed->[ chromStart ] =~ /\D/ ) {
 207         Maasha::Common::error( qq(Bad BED entry - chromStart must be a whole number - not "$bed->[ chromStart ]") );
 208     }
 209
 210     if ( $bed->[ chromEnd ] =~ /\D/ ) {
 211         Maasha::Common::error( qq(Bad BED entry - chromEnd must be a whole number - not "$bed->[ chromEnd ]") );
 212     }
 213
 214     if ( $bed->[ chromEnd ] < $bed->[ chromStart ] ) {
 215         Maasha::Common::error( qq(Bad BED entry - chromEnd must be greater than chromStart - not "$bed->[ chromStart ] > $bed->[ chromEnd ]") );
 216     }
 217
 218     return if $cols == 3;
 219
 220     if ( $bed->[ name ] =~ /\s/ ) {
 221         Maasha::Common::error( qq(Bad BED entry - no white space allowed in name field: "$bed->[ name ]") );
 222     }
 223
 224     return if $cols == 4;
 225
 226     if ( $bed->[ score ] =~ /\D/ ) {
 227         Maasha::Common::error( qq(Bad BED entry - score must be a whole number - not "$bed->[ score ]") );
 228     }
 229
 230     # if ( $bed->[ score ] < 0 or $bed->[ score ] > 1000 ) { # disabled - too restrictive !
 231     if ( $bed->[ score ] < 0 ) {
 232         Maasha::Common::error( qq(Bad BED entry - score must be between 0 and 1000 - not "$bed->[ score ]") );
 233     }
 234
 235     return if $cols == 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 $cols == 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 $cols == 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 $cols == 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->[ chromStart ];
 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     if ( defined $bp_record->{ "CHR_BEG" } ) {
 483         $bed_entry[ chromStart ] = $bp_record->{ "CHR_BEG" };
 484     } elsif ( defined $bp_record->{ "S_BEG" } ) {
 485         $bed_entry[ chromStart ] = $bp_record->{ "S_BEG" };
 486     } else {
 487         return undef;
 488     }
 489
 490     if ( defined $bp_record->{ "CHR_END" } ) {
 491         $bed_entry[ chromEnd ] = $bp_record->{ "CHR_END" };
 492     } elsif ( defined $bp_record->{ "S_END" }) {
 493         $bed_entry[ chromEnd ] = $bp_record->{ "S_END" };
 494     } else {
 495         return undef;
 496     }
 497
 498     $bed_entry[ chromEnd ]++;
 499
 500     return wantarray ? @bed_entry : \@bed_entry if $cols == 3;
 501
 502     $bed_entry[ name ] = $bp_record->{ "Q_ID" } || return wantarray ? @bed_entry : \@bed_entry;
 503
 504     return wantarray ? @bed_entry : \@bed_entry if $cols == 4;
 505
 506     if ( exists $bp_record->{ "SCORE" } ) {
 507         $bed_entry[ score ] = $bp_record->{ "SCORE" };
 508     } else {
 509         return wantarray ? @bed_entry : \@bed_entry;
 510     }
 511
 512     return wantarray ? @bed_entry : \@bed_entry if $cols == 5;
 513
 514     if ( exists $bp_record->{ "STRAND" } ) {
 515         $bed_entry[ strand ] = $bp_record->{ "STRAND" };
 516     } else {
 517         return wantarray ? @bed_entry : \@bed_entry;
 518     }
 519
 520     return wantarray ? @bed_entry : \@bed_entry if $cols == 6;
 521
 522     if ( defined $bp_record->{ "THICK_BEG" }   and
 523          defined $bp_record->{ "THICK_END" } )
 524     {
 525         $bed_entry[ thickStart ] = $bp_record->{ "THICK_BEG" };
 526         $bed_entry[ thickEnd ]   = $bp_record->{ "THICK_END" };
 527     }
 528
 529     return wantarray ? @bed_entry : \@bed_entry if $cols == 8;
 530
 531     if ( exists $bp_record->{ "COLOR" } ) {
 532         $bed_entry[ itemRgb ] = $bp_record->{ "COLOR" };
 533     } else {
 534         return wantarray ? @bed_entry : \@bed_entry;
 535     }
 536
 537     return wantarray ? @bed_entry : \@bed_entry if $cols == 9;
 538
 539     if ( defined $bp_record->{ "BLOCK_COUNT" } and
 540          defined $bp_record->{ "BLOCK_LENS" }  and
 541          defined $bp_record->{ "Q_BEGS" } )
 542     {
 543         $bed_entry[ blockCount ]  = $bp_record->{ "BLOCK_COUNT" };
 544         $bed_entry[ blockSizes ]  = $bp_record->{ "BLOCK_LENS" };
 545         $bed_entry[ blockStarts ] = $bp_record->{ "Q_BEGS" };
 546         $bed_entry[ thickEnd ]++;
 547     }
 548
 549     return wantarray ? @bed_entry : \@bed_entry;
 550 }
 551
 552
 553 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> TAG CONTIGS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 554
 555
 556 sub tag_contigs_assemble
 557 {
 558     # Martin A. Hansen, November 2008.
 559
 560     # Returns a path to a BED file with tab contigs.
 561
 562     my ( $bed_file,   # path to BED file
 563          $chr,        # chromosome
 564          $strand,     # strand
 565          $dir,        # working directory
 566        ) = @_;
 567
 568     # Returns a string
 569
 570     my ( $fh_in, $fh_out, $array, $new_file, $pos, $id, $beg, $end, $score, $entry );
 571
 572     $fh_in = Maasha::Filesys::file_read_open( $bed_file );
 573
 574     $array = tag_contigs_assemble_array( $fh_in );
 575
 576     close $fh_in;
 577
 578     $new_file = "$bed_file.tag_contigs";
 579
 580     $fh_out = Maasha::Filesys::file_write_open( $new_file );
 581
 582     $pos = 0;
 583     $id  = 0;
 584
 585     while ( ( $beg, $end, $score ) = tag_contigs_scan( $array, $pos ) and $beg )
 586     {
 587         $entry->[ chrom ]      = $chr;
 588         $entry->[ chromStart ] = $beg;
 589         $entry->[ chromEnd ]   = $end;
 590         $entry->[ name ]       = sprintf( "TC%06d", $id );
 591         $entry->[ score ]      = $score;
 592         $entry->[ strand ]     = $strand;
 593
 594         bed_entry_put( $entry, $fh_out );
 595
 596         $pos = $end;
 597         $id++;
 598     }
 599
 600     close $fh_out;
 601
 602     return $new_file;
 603 }
 604
 605
 606 sub tag_contigs_assemble_array
 607 {
 608     # Martin A. Hansen, November 2008.
 609
 610     # Given a BED file with entries from only one
 611     # chromosome assembles tag contigs from these
 612     # ignoring strand information. Only tags with
 613     # a score higher than the clone count over
 614     # genomic loci (the score field) is included
 615     # in the tag contigs.
 616
 617     #       -----------              tags
 618     #          -------------
 619     #               ----------
 620     #                     ----------
 621     #       ========================  tag contig
 622
 623
 624     my ( $fh,   # file handle to BED file
 625        ) = @_;
 626
 627     # Returns an arrayref.
 628
 629     my ( $entry, $clones, $score, @array );
 630
 631     while ( $entry = bed_entry_get( $fh ) )
 632     {
 633         if ( $entry->[ name ] =~ /_(\d+)$/ )
 634         {
 635             $clones = $1;
 636
 637             if ( $entry->[ score ] )
 638             {
 639                 $score = int( $clones / $entry->[ score ] );
 640
 641                 map { $array[ $_ ] += $score } $entry->[ chromStart ] .. $entry->[ chromEnd ] - 1 if $score >= 1;
 642             }
 643         }
 644     }
 645
 646     return wantarray ? @array : \@array;
 647 }
 648
 649
 650 sub tag_contigs_scan
 651 {
 652     # Martin A. Hansen, November 2008.
 653
 654     # Scans an array with tag contigs and locates
 655     # the next contig from a given position. The
 656     # score of the tag contig is determined as the
 657     # maximum value of the tag contig. If a tag contig
 658     # is found a triple is returned with beg, end and score
 659     # otherwise an empty list is returned.
 660
 661     my ( $array,   # array to scan
 662          $beg,     # position to start scanning from
 663        ) = @_;
 664
 665     # Returns an arrayref.
 666
 667     my ( $end, $score );
 668
 669     $score = 0;
 670
 671     while ( $beg < scalar @{ $array } and not $array->[ $beg ] ) { $beg++ }
 672
 673     $end = $beg;
 674
 675     while ( $array->[ $end ] )
 676     {
 677         $score = Maasha::Calc::max( $score, $array->[ $end ] );
 678
 679         $end++
 680     }
 681
 682     if ( $score > 0 ) {
 683         return wantarray ? ( $beg, $end, $score ) : [ $beg, $end, $score ];
 684     } else {
 685         return wantarray ? () : [];
 686     }
 687 }
 688
 689
 690 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 691
 692
 693 1;
 694
 695
 696 __END__