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