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 ( $entry, $cols, @block_sizes, @block_starts );
 192
 193     $entry = join "\t", @{ $bed };
 194     $cols  = scalar @{ $bed };
 195
 196     if ( $cols < 3 ) {
 197         Maasha::Common::error( qq(Bad BED entry "$entry" - must contain at least 3 columns - not $cols) );
 198     }
 199
 200     if ( $cols > 12 ) {
 201         Maasha::Common::error( qq(Bad BED entry "$entry" - must contains no more than 12 columns - not $cols) );
 202     }
 203
 204     if ( $bed->[ chrom ] =~ /\s/ ) {
 205         Maasha::Common::error( qq(Bad BED entry "$entry" - no white space allowed in chrom field: "$bed->[ chrom ]") );
 206     }
 207
 208     if ( $bed->[ chromStart ] =~ /\D/ ) {
 209         Maasha::Common::error( qq(Bad BED entry "$entry" - chromStart must be a whole number - not "$bed->[ chromStart ]") );
 210     }
 211
 212     if ( $bed->[ chromEnd ] =~ /\D/ ) {
 213         Maasha::Common::error( qq(Bad BED entry "$entry" - chromEnd must be a whole number - not "$bed->[ chromEnd ]") );
 214     }
 215
 216     if ( $bed->[ chromEnd ] < $bed->[ chromStart ] ) {
 217         Maasha::Common::error( qq(Bad BED entry "$entry" - chromEnd must be greater than chromStart - not "$bed->[ chromStart ] > $bed->[ chromEnd ]") );
 218     }
 219
 220     return if $cols == 3;
 221
 222     if ( $bed->[ name ] =~ /\s/ ) {
 223         Maasha::Common::error( qq(Bad BED entry "$entry" - no white space allowed in name field: "$bed->[ name ]") );
 224     }
 225
 226     return if $cols == 4;
 227
 228     if ( $bed->[ score ] =~ /\D/ ) {
 229         Maasha::Common::error( qq(Bad BED entry "$entry" - score must be a whole number - not "$bed->[ score ]") );
 230     }
 231
 232     # if ( $bed->[ score ] < 0 or $bed->[ score ] > 1000 ) { # disabled - too restrictive !
 233     if ( $bed->[ score ] < 0 ) {
 234         Maasha::Common::error( qq(Bad BED entry "$entry" - score must be between 0 and 1000 - not "$bed->[ score ]") );
 235     }
 236
 237     return if $cols == 5;
 238
 239     if ( $bed->[ strand ] ne '+' and $bed->[ strand ] ne '-' ) {
 240         Maasha::Common::error( qq(Bad BED entry "$entry" - strand must be + or - not "$bed->[ strand ]") );
 241     }
 242
 243     return if $cols == 6;
 244
 245     if ( $bed->[ thickStart ] =~ /\D/ ) {
 246         Maasha::Common::error( qq(Bad BED entry "$entry" - thickStart must be a whole number - not "$bed->[ thickStart ]") );
 247     }
 248
 249     if ( $bed->[ thickEnd ] =~ /\D/ ) {
 250         Maasha::Common::error( qq(Bad BED entry "$entry" - thickEnd must be a whole number - not "$bed->[ thickEnd ]") );
 251     }
 252
 253     if ( $bed->[ thickEnd ] < $bed->[ thickStart ] ) {
 254         Maasha::Common::error( qq(Bad BED entry "$entry" - thickEnd must be greater than thickStart - not "$bed->[ thickStart ] > $bed->[ thickEnd ]") );
 255     }
 256
 257     if ( $bed->[ thickStart ] < $bed->[ chromStart ] ) {
 258         Maasha::Common::error( qq(Bad BED entry "$entry" - thickStart must be greater than chromStart - not "$bed->[ thickStart ] < $bed->[ chromStart ]") );
 259     }
 260
 261     if ( $bed->[ thickStart ] > $bed->[ chromEnd ] ) {
 262         Maasha::Common::error( qq(Bad BED entry "$entry" - thickStart must be less than chromEnd - not "$bed->[ thickStart ] > $bed->[ chromEnd ]") );
 263     }
 264
 265     if ( $bed->[ thickEnd ] < $bed->[ chromStart ] ) {
 266         Maasha::Common::error( qq(Bad BED entry "$entry" - thickEnd must be greater than chromStart - not "$bed->[ thickEnd ] < $bed->[ chromStart ]") );
 267     }
 268
 269     if ( $bed->[ thickEnd ] > $bed->[ chromEnd ] ) {
 270         Maasha::Common::error( qq(Bad BED entry "$entry" - thickEnd must be less than chromEnd - not "$bed->[ thickEnd ] > $bed->[ chromEnd ]") );
 271     }
 272
 273     return if $cols == 8;
 274
 275     if ( $bed->[ itemRgb ] !~ /^(0|\d{1,3},\d{1,3},\d{1,3},?)$/ ) {
 276         Maasha::Common::error( qq(Bad BED entry "$entry" - itemRgb must be 0 or in the form of 255,0,0 - not "$bed->[ itemRgb ]") );
 277     }
 278
 279     return if $cols == 9;
 280
 281     if ( $bed->[ blockCount ] =~ /\D/ ) {
 282         Maasha::Common::error( qq(Bad BED entry "$entry" - blockCount must be a whole number - not "$bed->[ blockCount ]") );
 283     }
 284
 285     @block_sizes = split ",", $bed->[ blockSizes ];
 286
 287     if ( grep /\D/, @block_sizes ) {
 288         Maasha::Common::error( qq(Bad BED entry "$entry" - blockSizes must be whole numbers - not "$bed->[ blockSizes ]") );
 289     }
 290
 291     if ( $bed->[ blockCount ] != scalar @block_sizes ) {
 292         Maasha::Common::error( qq(Bad BED entry "$entry" - blockSizes "$bed->[ blockSizes ]" must match blockCount "$bed->[ blockCount ]") );
 293     }
 294
 295     @block_starts = split ",", $bed->[ blockStarts ];
 296
 297     if ( grep /\D/, @block_starts ) {
 298         Maasha::Common::error( qq(Bad BED entry "$entry" - blockStarts must be whole numbers - not "$bed->[ blockStarts ]") );
 299     }
 300
 301     if ( $bed->[ blockCount ] != scalar @block_starts ) {
 302         Maasha::Common::error( qq(Bad BED entry "$entry" - blockStarts "$bed->[ blockStarts ]" must match blockCount "$bed->[ blockCount ]") );
 303     }
 304
 305     if ( $bed->[ chromStart ] + $block_starts[ -1 ] + $block_sizes[ -1 ] != $bed->[ chromEnd ] ) {
 306         Maasha::Common::error( qq(Bad BED entry "$entry" - chromStart + blockStarts[last] + blockSizes[last] must equal chromEnd: ) .
 307                                qq($bed->[ chromStart ] + $block_starts[ -1 ] + $block_sizes[ -1 ] != $bed->[ chromEnd ]) );
 308     }
 309 }
 310
 311
 312 sub bed_sort
 313 {
 314     # Martin A. hansen, November 2008.
 315
 316     # Sorts a given BED file according to a given
 317     # sorting scheme:
 318     # 1: chr AND chr_beg.
 319     # 2: chr AND strand AND chr_beg.
 320     # 3: chr_beg.
 321     # 4: strand AND chr_beg.
 322
 323     # Currently 'bed_sort' is used for sorting = a standalone c program
 324     # that uses qsort (unstable sort).
 325
 326     my ( $file_in,  # path to BED file.
 327          $scheme,   # sort scheme
 328          $cols,     # number of columns in BED file
 329        ) = @_;
 330
 331     # Returns nothing.
 332
 333     my ( $file_out );
 334
 335     $file_out = "$file_in.sort";
 336
 337     Maasha::Common::run( "bed_sort", "--sort $scheme --cols $cols $file_in > $file_out" );
 338
 339     if ( Maasha::Filesys::file_size( $file_in ) !=  Maasha::Filesys::file_size( $file_out ) ) {
 340         Maasha::Common::error( qq(bed_sort of file "$file_in" failed) );
 341     }
 342
 343     rename $file_out, $file_in;
 344 }
 345
 346
 347 sub bed_file_split_on_chr
 348 {
 349     # Martin A. Hansen, November 2008.
 350
 351     # Given a path to a BED file splits
 352     # this file into one file per chromosome
 353     # in a temporary directory. Returns a hash
 354     # with chromosome name as key and the corresponding
 355     # file as value.
 356
 357     my ( $file,   # path to BED file
 358          $dir,    # working directory
 359          $cols,   # number of BED columns to read - OPTIONAL
 360        ) = @_;
 361
 362     # Returns a hashref
 363
 364     my ( $fh_in, $fh_out, $entry, %fh_hash, %file_hash );
 365
 366     $fh_in = Maasha::Filesys::file_read_open( $file );
 367
 368     while ( $entry = bed_entry_get( $fh_in, $cols ) )
 369     {
 370         if ( not exists $file_hash{ $entry->[ chrom ] } )
 371         {
 372             $fh_hash{ $entry->[ chrom ] }   = Maasha::Filesys::file_write_open( "$dir/$entry->[ chrom ]" );
 373             $file_hash{ $entry->[ chrom ] } = "$dir/$entry->[ chrom ]";
 374         }
 375
 376         $fh_out = $fh_hash{ $entry->[ chrom ] };
 377
 378         Maasha::UCSC::BED::bed_entry_put( $entry, $fh_out, $cols );
 379     }
 380
 381     map { close $_ } keys %fh_hash;
 382
 383     return wantarray ? %file_hash : \%file_hash;
 384 }
 385
 386
 387 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> BIOPIECES <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 388
 389
 390 sub bed2biopiece
 391 {
 392     # Martin A. Hansen, November 2008.
 393
 394     # Converts a BED entry given as an arrayref
 395     # to a Biopiece record which is returned as
 396     # a hashref.
 397
 398     # IMPORTANT! The BED_END and THICK_END positions
 399     # will be the exact position in contrary to the
 400     # UCSC scheme.
 401
 402     my ( $bed_entry,   # BED entry as arrayref
 403        ) = @_;
 404
 405     # Returns a hashref
 406
 407     my ( $cols, %bp_record );
 408
 409     $cols = scalar @{ $bed_entry };
 410
 411     if ( not defined $bed_entry->[ chrom ] and
 412          not defined $bed_entry->[ chromStart ] and
 413          not defined $bed_entry->[ chromEnd ] )
 414     {
 415         return 0;
 416     }
 417
 418     $bp_record{ "REC_TYPE" } = "BED";
 419     $bp_record{ "BED_COLS" } = $cols;
 420     $bp_record{ "CHR" }      = $bed_entry->[ chrom ];
 421     $bp_record{ "CHR_BEG" }  = $bed_entry->[ chromStart ];
 422     $bp_record{ "CHR_END" }  = $bed_entry->[ chromEnd ] - 1;
 423     $bp_record{ "BED_LEN" }  = $bed_entry->[ chromEnd ] - $bed_entry->[ chromStart ];
 424
 425     return wantarray ? %bp_record : \%bp_record if $cols == 3;
 426
 427     $bp_record{ "Q_ID" } = $bed_entry->[ name ];
 428
 429     return wantarray ? %bp_record : \%bp_record if $cols == 4;
 430
 431     $bp_record{ "SCORE" } = $bed_entry->[ score ];
 432
 433     return wantarray ? %bp_record : \%bp_record if $cols == 5;
 434
 435     $bp_record{ "STRAND" } = $bed_entry->[ strand ];
 436
 437     return wantarray ? %bp_record : \%bp_record if $cols == 6;
 438
 439     $bp_record{ "THICK_BEG" }   = $bed_entry->[ thickStart ];
 440     $bp_record{ "THICK_END" }   = $bed_entry->[ thickEnd ] - 1;
 441
 442     return wantarray ? %bp_record : \%bp_record if $cols == 8;
 443
 444     $bp_record{ "COLOR" }       = $bed_entry->[ itemRgb ];
 445
 446     return wantarray ? %bp_record : \%bp_record if $cols == 9;
 447
 448     $bp_record{ "BLOCK_COUNT" } = $bed_entry->[ blockCount ];
 449     $bp_record{ "BLOCK_LENS" }  = $bed_entry->[ blockSizes ];
 450     $bp_record{ "Q_BEGS" }      = $bed_entry->[ blockStarts ];
 451
 452     return wantarray ? %bp_record : \%bp_record;
 453 }
 454
 455
 456 sub biopiece2bed
 457 {
 458     # Martin A. Hansen, November 2008.
 459
 460     # Converts a Biopieces record given as a hashref
 461     # to a BED record which is returned as
 462     # an arrayref. As much as possible of the Biopiece
 463     # record is converted and undef is returned if
 464     # convertion failed.
 465
 466     # IMPORTANT! The chromEnd and thickEnd positions
 467     # will be the inexact position used in the
 468     # UCSC scheme (+1 to all ends).
 469
 470     my ( $bp_record,   # hashref
 471          $cols,        # number of columns in BED entry - OPTIONAL (but faster)
 472        ) = @_;
 473
 474     # Returns arrayref.
 475
 476     my ( @bed_entry, @begs );
 477
 478     $cols ||= 12;   # max number of columns possible
 479
 480     $bed_entry[ chrom ] = $bp_record->{ "CHR" }  ||
 481                           $bp_record->{ "S_ID" } ||
 482                           return undef;
 483
 484     if ( defined $bp_record->{ "CHR_BEG" } ) {
 485         $bed_entry[ chromStart ] = $bp_record->{ "CHR_BEG" };
 486     } elsif ( defined $bp_record->{ "S_BEG" } ) {
 487         $bed_entry[ chromStart ] = $bp_record->{ "S_BEG" };
 488     } else {
 489         return undef;
 490     }
 491
 492     if ( defined $bp_record->{ "CHR_END" } ) {
 493         $bed_entry[ chromEnd ] = $bp_record->{ "CHR_END" } + 1;
 494     } elsif ( defined $bp_record->{ "S_END" }) {
 495         $bed_entry[ chromEnd ] = $bp_record->{ "S_END" } + 1;
 496     } else {
 497         return undef;
 498     }
 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" } + 1;
 527     }
 528     elsif ( defined $bp_record->{ "BLOCK_COUNT" } )
 529     {
 530         $bed_entry[ thickStart ] = $bed_entry[ chromStart ];
 531         $bed_entry[ thickEnd ]   = $bed_entry[ chromEnd ];
 532     }
 533
 534     return wantarray ? @bed_entry : \@bed_entry if $cols == 8;
 535
 536     if ( defined $bp_record->{ "COLOR" } )
 537     {
 538         $bed_entry[ itemRgb ] = $bp_record->{ "COLOR" };
 539     }
 540     elsif ( defined $bp_record->{ "BLOCK_COUNT" } )
 541     {
 542         $bed_entry[ itemRgb ] = 0;
 543     }
 544
 545     return wantarray ? @bed_entry : \@bed_entry if $cols == 9;
 546
 547     if ( defined $bp_record->{ "BLOCK_COUNT" } and
 548          defined $bp_record->{ "BLOCK_LENS" }  and
 549          defined $bp_record->{ "S_BEGS" } )
 550     {
 551         @begs = split ",", $bp_record->{ "S_BEGS" };
 552         map { $_ -= $bed_entry[ chromStart ] } @begs;
 553
 554         $bed_entry[ blockCount ]  = $bp_record->{ "BLOCK_COUNT" };
 555         $bed_entry[ blockSizes ]  = $bp_record->{ "BLOCK_LENS" };
 556         $bed_entry[ blockStarts ] = join ",", @begs;
 557         # $bed_entry[ thickEnd ]    = $bp_record->{ "THICK_END" } + 1;
 558     }
 559     elsif ( defined $bp_record->{ "BLOCK_COUNT" } and
 560             defined $bp_record->{ "BLOCK_LENS" }  and
 561             defined $bp_record->{ "Q_BEGS" } )
 562     {
 563         $bed_entry[ blockCount ]  = $bp_record->{ "BLOCK_COUNT" };
 564         $bed_entry[ blockSizes ]  = $bp_record->{ "BLOCK_LENS" };
 565         $bed_entry[ blockStarts ] = $bp_record->{ "Q_BEGS" };
 566         # $bed_entry[ thickEnd  ]   = $bp_record->{ "THICK_END" } + 1;
 567     }
 568
 569     return wantarray ? @bed_entry : \@bed_entry;
 570 }
 571
 572
 573 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> TAG CONTIGS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 574
 575
 576 sub tag_contigs_assemble
 577 {
 578     # Martin A. Hansen, November 2008.
 579
 580     # Returns a path to a BED file with tab contigs.
 581
 582     my ( $bed_file,   # path to BED file
 583          $chr,        # chromosome
 584          $strand,     # strand
 585          $dir,        # working directory
 586        ) = @_;
 587
 588     # Returns a string
 589
 590     my ( $fh_in, $fh_out, $array, $new_file, $pos, $id, $beg, $end, $score, $entry );
 591
 592     $fh_in = Maasha::Filesys::file_read_open( $bed_file );
 593
 594     $array = tag_contigs_assemble_array( $fh_in );
 595
 596     close $fh_in;
 597
 598     $new_file = "$bed_file.tag_contigs";
 599
 600     $fh_out = Maasha::Filesys::file_write_open( $new_file );
 601
 602     $pos = 0;
 603     $id  = 0;
 604
 605     while ( ( $beg, $end, $score ) = tag_contigs_scan( $array, $pos ) and $beg )
 606     {
 607         $entry->[ chrom ]      = $chr;
 608         $entry->[ chromStart ] = $beg;
 609         $entry->[ chromEnd ]   = $end;
 610         $entry->[ name ]       = sprintf( "TC%06d", $id );
 611         $entry->[ score ]      = $score;
 612         $entry->[ strand ]     = $strand;
 613
 614         bed_entry_put( $entry, $fh_out );
 615
 616         $pos = $end;
 617         $id++;
 618     }
 619
 620     close $fh_out;
 621
 622     return $new_file;
 623 }
 624
 625
 626 sub tag_contigs_assemble_array
 627 {
 628     # Martin A. Hansen, November 2008.
 629
 630     # Given a BED file with entries from only one
 631     # chromosome assembles tag contigs from these
 632     # ignoring strand information. Only tags with
 633     # a score higher than the clone count over
 634     # genomic loci (the score field) is included
 635     # in the tag contigs.
 636
 637     #       -----------              tags
 638     #          -------------
 639     #               ----------
 640     #                     ----------
 641     #       ========================  tag contig
 642
 643
 644     my ( $fh,   # file handle to BED file
 645        ) = @_;
 646
 647     # Returns an arrayref.
 648
 649     my ( $entry, $clones, $score, @array );
 650
 651     while ( $entry = bed_entry_get( $fh ) )
 652     {
 653         if ( $entry->[ name ] =~ /_(\d+)$/ )
 654         {
 655             $clones = $1;
 656
 657             if ( $entry->[ score ] )
 658             {
 659                 $score = int( $clones / $entry->[ score ] );
 660
 661                 map { $array[ $_ ] += $score } $entry->[ chromStart ] .. $entry->[ chromEnd ] - 1 if $score >= 1;
 662             }
 663         }
 664     }
 665
 666     return wantarray ? @array : \@array;
 667 }
 668
 669
 670 sub tag_contigs_scan
 671 {
 672     # Martin A. Hansen, November 2008.
 673
 674     # Scans an array with tag contigs and locates
 675     # the next contig from a given position. The
 676     # score of the tag contig is determined as the
 677     # maximum value of the tag contig. If a tag contig
 678     # is found a triple is returned with beg, end and score
 679     # otherwise an empty list is returned.
 680
 681     my ( $array,   # array to scan
 682          $beg,     # position to start scanning from
 683        ) = @_;
 684
 685     # Returns an arrayref.
 686
 687     my ( $end, $score );
 688
 689     $score = 0;
 690
 691     while ( $beg < scalar @{ $array } and not $array->[ $beg ] ) { $beg++ }
 692
 693     $end = $beg;
 694
 695     while ( $array->[ $end ] )
 696     {
 697         $score = Maasha::Calc::max( $score, $array->[ $end ] );
 698
 699         $end++
 700     }
 701
 702     if ( $score > 0 ) {
 703         return wantarray ? ( $beg, $end, $score ) : [ $beg, $end, $score ];
 704     } else {
 705         return wantarray ? () : [];
 706     }
 707 }
 708
 709
 710 sub bed_upload_to_ucsc
 711 {
 712     # Martin A. Hansen, September 2007.
 713
 714     # Upload a BED file to the UCSC database.
 715
 716     my ( $tmp_dir,   # temporary directory
 717          $file,      # file to upload,
 718          $options,   # argument hashref
 719          $append,    # flag indicating table should be appended
 720        ) = @_;
 721
 722     # Returns nothing.
 723
 724     my ( $args, $table, $sql_file, $fh_out, $fh_in );
 725
 726     if ( $append ) {
 727         $args = join " ", $options->{ "database" }, $options->{ "table" }, "-tmpDir=$tmp_dir", "-oldTable", $file;
 728     } else {
 729         $args = join " ", $options->{ "database" }, $options->{ "table" }, "-tmpDir=$tmp_dir", $file;
 730     }
 731
 732     if ( $options->{ "table" } =~ /rnaSecStr/ )
 733     {
 734         $table = $options->{ "table" };
 735
 736         print qq(uploading secondary structure table:"$table"\n) if $options->{ "verbose" };
 737
 738         $sql_file = "$tmp_dir/upload_RNA_SS.sql";
 739
 740         $fh_out   = Maasha::Filesys::file_write_open( $sql_file );
 741
 742         print $fh_out qq(
 743 CREATE TABLE $table (
 744     bin smallint not null,              # Bin number for browser speedup
 745     chrom varchar(255) not null,        # Chromosome or FPC contig
 746     chromStart int unsigned not null,   # Start position in chromosome
 747     chromEnd int unsigned not null,     # End position in chromosome
 748     name varchar(255) not null,         # Name of item
 749     score int unsigned not null,        # Score from 0-1000
 750     strand char(1) not null,            # + or -
 751     size int unsigned not null,         # Size of element.
 752     secStr longblob not null,           # Parentheses and '.'s which define the secondary structure
 753     conf longblob not null,             # Confidence of secondary-structure annotation per position (0.0-1.0).
 754     #Indices
 755     INDEX(name(16)),
 756     INDEX(chrom(8), bin),
 757     INDEX(chrom(8), chromStart)
 758 );
 759         );
 760
 761         close $fh_out;
 762
 763         Maasha::Common::run( "hgLoadBed", "-notItemRgb -sqlTable=$sql_file $options->{ 'database' } $options->{ 'table' } -tmpDir=$tmp_dir $file > /dev/null 2>&1" );
 764
 765         unlink $sql_file;
 766     }
 767     else
 768     {
 769         Maasha::Common::run( "hgLoadBed", "$args > /dev/null 2>&1" );
 770     }
 771 }
 772
 773
 774 sub bed_analyze
 775 {
 776     # Martin A. Hansen, March 2008.
 777
 778     # Given a bed record, analysis this to give information
 779     # about intron/exon sizes.
 780
 781     my ( $entry,   # BED entry
 782        ) = @_;
 783
 784     # Returns hashref.
 785
 786     my ( $i, @begs, @lens, $exon_max, $exon_min, $exon_len, $exon_tot, $intron_max, $intron_min, $intron_len, $intron_tot );
 787
 788     $exon_max   = 0;
 789     $exon_min   = 9999999999;
 790     $intron_max = 0;
 791     $intron_min = 9999999999;
 792
 793     $entry->{ "EXONS" }   = $entry->{ "BLOCK_COUNT" };
 794
 795     @begs = split /,/, $entry->{ "Q_BEGS" };
 796     @lens = split /,/, $entry->{ "BLOCK_LENS" };
 797
 798     for ( $i = 0; $i < $entry->{ "BLOCK_COUNT" }; $i++ )
 799     {
 800         $exon_len = $lens[ $i ];
 801
 802         $entry->{ "EXON_LEN_$i" } = $exon_len;
 803
 804         $exon_max = $exon_len if $exon_len > $exon_max;
 805         $exon_min = $exon_len if $exon_len < $exon_min;
 806
 807         $exon_tot += $exon_len;
 808     }
 809
 810     $entry->{ "EXON_LEN_-1" }   = $exon_len;
 811     $entry->{ "EXON_MAX_LEN" }  = $exon_max;
 812     $entry->{ "EXON_MIN_LEN" }  = $exon_min;
 813     $entry->{ "EXON_MEAN_LEN" } = int( $exon_tot / $entry->{ "EXONS" } );
 814
 815     $entry->{ "INTRONS" } = $entry->{ "BLOCK_COUNT" } - 1;
 816     $entry->{ "INTRONS" } = 0 if $entry->{ "INTRONS" } < 0;
 817
 818     if ( $entry->{ "INTRONS" } )
 819     {
 820         for ( $i = 1; $i < $entry->{ "BLOCK_COUNT" }; $i++ )
 821         {
 822             $intron_len = $begs[ $i ] - ( $begs[ $i - 1 ] + $lens[ $i - 1 ] );
 823
 824             $entry->{ "INTRON_LEN_" . ( $i - 1 ) } = $intron_len;
 825
 826             $intron_max = $intron_len if $intron_len > $intron_max;
 827             $intron_min = $intron_len if $intron_len < $intron_min;
 828
 829             $intron_tot += $intron_len;
 830         }
 831
 832         $entry->{ "INTRON_LEN_-1" }   = $intron_len;
 833         $entry->{ "INTRON_MAX_LEN" }  = $intron_max;
 834         $entry->{ "INTRON_MIN_LEN" }  = $intron_min;
 835         $entry->{ "INTRON_MEAN_LEN" } = int( $intron_tot / $entry->{ "INTRONS" } );
 836     }
 837
 838     return wantarray ? %{ $entry } : $entry;
 839 }
 840
 841
 842 sub bed_merge_entries
 843 {
 844     # Martin A. Hansen, February 2008.
 845
 846     # Merge a list of given BED entries in one big entry.
 847
 848     my ( $entries,     # list of BED entries to be merged
 849        ) = @_;
 850
 851     # Returns hash.
 852
 853     my ( $i, @q_ids, @q_begs, @blocksizes, @new_q_begs, @new_blocksizes, %new_entry );
 854
 855     @{ $entries } = sort { $a->{ "CHR_BEG" } <=> $b->{ "CHR_BEG" } } @{ $entries };
 856
 857     for ( $i = 0; $i < @{ $entries }; $i++ )
 858     {
 859         Maasha::Common::error( qq(Attempted merge of BED entries from different chromosomes) ) if $entries->[ 0 ]->{ "CHR" }    ne $entries->[ $i ]->{ "CHR" };
 860         Maasha::Common::error( qq(Attempted merge of BED entries from different strands) )     if $entries->[ 0 ]->{ "STRAND" } ne $entries->[ $i ]->{ "STRAND" };
 861
 862         push @q_ids, $entries->[ $i ]->{ "Q_ID" } || sprintf( "ID%06d", $i );
 863
 864         if ( exists $entries->[ $i ]->{ "Q_BEGS" } )
 865         {
 866             @q_begs     = split ",", $entries->[ $i ]->{ "Q_BEGS" };
 867             @blocksizes = split ",", $entries->[ $i ]->{ "BLOCK_LENS" };
 868         }
 869         else
 870         {
 871             @q_begs     = 0;
 872             @blocksizes = $entries->[ $i ]->{ "CHR_END" } - $entries->[ $i ]->{ "CHR_BEG" } + 1;
 873         }
 874
 875         map { $_ += $entries->[ $i ]->{ "CHR_BEG" } } @q_begs;
 876
 877         push @new_q_begs, @q_begs;
 878         push @new_blocksizes, @blocksizes;
 879     }
 880
 881     map { $_ -= $entries->[ 0 ]->{ "CHR_BEG" } } @new_q_begs;
 882
 883     %new_entry = (
 884         CHR          => $entries->[ 0 ]->{ "CHR" },
 885         CHR_BEG      => $entries->[ 0 ]->{ "CHR_BEG" },
 886         CHR_END      => $entries->[ -1 ]->{ "CHR_END" },
 887         REC_TYPE     => "BED",
 888         BED_LEN      => $entries->[ -1 ]->{ "CHR_END" } - $entries->[ 0 ]->{ "CHR_BEG" } + 1,
 889         BED_COLS     => 12,
 890         Q_ID         => join( ":", @q_ids ),
 891         SCORE        => 999,
 892         STRAND       => $entries->[ 0 ]->{ "STRAND" }     || "+",
 893         THICK_BEG    => $entries->[ 0 ]->{ "THICK_BEG" }  || $entries->[ 0 ]->{ "CHR_BEG" },
 894         THICK_END    => $entries->[ -1 ]->{ "THICK_END" } || $entries->[ -1 ]->{ "CHR_END" },
 895         COLOR        => 0,
 896         BLOCK_COUNT  => scalar @new_q_begs,
 897         BLOCK_LENS   => join( ",", @new_blocksizes ),
 898         Q_BEGS       => join( ",", @new_q_begs ),
 899     );
 900
 901     return wantarray ? %new_entry : \%new_entry;
 902 }
 903
 904
 905 sub bed_split_entry
 906 {
 907     # Martin A. Hansen, February 2008.
 908
 909     # Splits a given BED entry into a list of blocks,
 910     # which are returned. A list of 6 column BED entry is returned.
 911
 912     my ( $entry,    # BED entry hashref
 913        ) = @_;
 914
 915     # Returns a list.
 916
 917     my ( @q_begs, @blocksizes, $block, @blocks, $i );
 918
 919     if ( exists $entry->{ "BLOCK_COUNT" } )
 920     {
 921         @q_begs     = split ",", $entry->{ "Q_BEGS" };
 922         @blocksizes = split ",", $entry->{ "BLOCK_LENS" };
 923
 924         for ( $i = 0; $i < @q_begs; $i++ )
 925         {
 926             undef $block;
 927
 928             $block->{ "CHR" }      = $entry->{ "CHR" };
 929             $block->{ "CHR_BEG" }  = $entry->{ "CHR_BEG" } + $q_begs[ $i ];
 930             $block->{ "CHR_END" }  = $entry->{ "CHR_BEG" } + $q_begs[ $i ] + $blocksizes[ $i ] - 1;
 931             $block->{ "Q_ID" }     = $entry->{ "Q_ID" } . sprintf( "_%03d", $i );
 932             $block->{ "SCORE" }    = $entry->{ "SCORE" };
 933             $block->{ "STRAND" }   = $entry->{ "STRAND" };
 934             $block->{ "BED_LEN" }  = $block->{ "CHR_END" } - $block->{ "CHR_BEG" } + 1,
 935             $block->{ "BED_COLS" } = 6;
 936             $block->{ "REC_TYPE" } = "BED";
 937
 938             push @blocks, $block;
 939         }
 940     }
 941     else
 942     {
 943         @blocks = @{ $entry };
 944     }
 945
 946     return wantarray ? @blocks : \@blocks;
 947 }
 948
 949
 950 sub bed_overlap
 951 {
 952     # Martin A. Hansen, February 2008.
 953
 954     # Checks if two BED entries overlap and
 955     # return 1 if so - else 0;
 956
 957     my ( $entry1,      # hashref
 958          $entry2,      # hashref
 959          $no_strand,   # don't check strand flag - OPTIONAL
 960        ) = @_;
 961
 962     # Return bolean.
 963
 964     return 0 if $entry1->{ "CHR" }    ne $entry2->{ "CHR" };
 965     return 0 if $entry1->{ "STRAND" } ne $entry2->{ "STRAND" };
 966
 967     if ( $entry1->{ "CHR_END" } < $entry2->{ "CHR_BEG" } or $entry1->{ "CHR_BEG" } > $entry2->{ "CHR_END" } ) {
 968         return 0;
 969     } else {
 970         return 1;
 971     }
 972 }
 973
 974
 975
 976
 977 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 978
 979
 980 1;
 981
 982
 983 __END__