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     } elsif ( exists $bp_record->{ "BIT_SCORE" } ) {
 509         $bed_entry[ score ] = $bp_record->{ "BIT_SCORE" };
 510     } else {
 511         return wantarray ? @bed_entry : \@bed_entry;
 512     }
 513
 514     return wantarray ? @bed_entry : \@bed_entry if $cols == 5;
 515
 516     if ( exists $bp_record->{ "STRAND" } ) {
 517         $bed_entry[ strand ] = $bp_record->{ "STRAND" };
 518     } else {
 519         return wantarray ? @bed_entry : \@bed_entry;
 520     }
 521
 522     return wantarray ? @bed_entry : \@bed_entry if $cols == 6;
 523
 524     if ( defined $bp_record->{ "THICK_BEG" }   and
 525          defined $bp_record->{ "THICK_END" } )
 526     {
 527         $bed_entry[ thickStart ] = $bp_record->{ "THICK_BEG" };
 528         $bed_entry[ thickEnd ]   = $bp_record->{ "THICK_END" } + 1;
 529     }
 530     elsif ( defined $bp_record->{ "BLOCK_COUNT" } )
 531     {
 532         $bed_entry[ thickStart ] = $bed_entry[ chromStart ];
 533         $bed_entry[ thickEnd ]   = $bed_entry[ chromEnd ];
 534     }
 535
 536     return wantarray ? @bed_entry : \@bed_entry if $cols == 8;
 537
 538     if ( defined $bp_record->{ "COLOR" } )
 539     {
 540         $bed_entry[ itemRgb ] = $bp_record->{ "COLOR" };
 541     }
 542     elsif ( defined $bp_record->{ "BLOCK_COUNT" } )
 543     {
 544         $bed_entry[ itemRgb ] = 0;
 545     }
 546
 547     return wantarray ? @bed_entry : \@bed_entry if $cols == 9;
 548
 549     if ( defined $bp_record->{ "BLOCK_COUNT" } and
 550          defined $bp_record->{ "BLOCK_LENS" }  and
 551          defined $bp_record->{ "S_BEGS" } )
 552     {
 553         @begs = split ",", $bp_record->{ "S_BEGS" };
 554         map { $_ -= $bed_entry[ chromStart ] } @begs;
 555
 556         $bed_entry[ blockCount ]  = $bp_record->{ "BLOCK_COUNT" };
 557         $bed_entry[ blockSizes ]  = $bp_record->{ "BLOCK_LENS" };
 558         $bed_entry[ blockStarts ] = join ",", @begs;
 559         # $bed_entry[ thickEnd ]    = $bp_record->{ "THICK_END" } + 1;
 560     }
 561     elsif ( defined $bp_record->{ "BLOCK_COUNT" } and
 562             defined $bp_record->{ "BLOCK_LENS" }  and
 563             defined $bp_record->{ "Q_BEGS" } )
 564     {
 565         $bed_entry[ blockCount ]  = $bp_record->{ "BLOCK_COUNT" };
 566         $bed_entry[ blockSizes ]  = $bp_record->{ "BLOCK_LENS" };
 567         $bed_entry[ blockStarts ] = $bp_record->{ "Q_BEGS" };
 568         # $bed_entry[ thickEnd  ]   = $bp_record->{ "THICK_END" } + 1;
 569     }
 570
 571     return wantarray ? @bed_entry : \@bed_entry;
 572 }
 573
 574
 575 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> TAG CONTIGS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 576
 577
 578 sub tag_contigs_assemble
 579 {
 580     # Martin A. Hansen, November 2008.
 581
 582     # Returns a path to a BED file with tab contigs.
 583
 584     my ( $bed_file,   # path to BED file
 585          $chr,        # chromosome
 586          $strand,     # strand
 587          $dir,        # working directory
 588        ) = @_;
 589
 590     # Returns a string
 591
 592     my ( $fh_in, $fh_out, $array, $new_file, $pos, $id, $beg, $end, $score, $entry );
 593
 594     $fh_in = Maasha::Filesys::file_read_open( $bed_file );
 595
 596     $array = tag_contigs_assemble_array( $fh_in );
 597
 598     close $fh_in;
 599
 600     $new_file = "$bed_file.tag_contigs";
 601
 602     $fh_out = Maasha::Filesys::file_write_open( $new_file );
 603
 604     $pos = 0;
 605     $id  = 0;
 606
 607     while ( ( $beg, $end, $score ) = tag_contigs_scan( $array, $pos ) and $beg )
 608     {
 609         $entry->[ chrom ]      = $chr;
 610         $entry->[ chromStart ] = $beg;
 611         $entry->[ chromEnd ]   = $end;
 612         $entry->[ name ]       = sprintf( "TC%06d", $id );
 613         $entry->[ score ]      = $score;
 614         $entry->[ strand ]     = $strand;
 615
 616         bed_entry_put( $entry, $fh_out );
 617
 618         $pos = $end;
 619         $id++;
 620     }
 621
 622     close $fh_out;
 623
 624     return $new_file;
 625 }
 626
 627
 628 sub tag_contigs_assemble_array
 629 {
 630     # Martin A. Hansen, November 2008.
 631
 632     # Given a BED file with entries from only one
 633     # chromosome assembles tag contigs from these
 634     # ignoring strand information. Only tags with
 635     # a score higher than the clone count over
 636     # genomic loci (the score field) is included
 637     # in the tag contigs.
 638
 639     #       -----------              tags
 640     #          -------------
 641     #               ----------
 642     #                     ----------
 643     #       ========================  tag contig
 644
 645
 646     my ( $fh,   # file handle to BED file
 647        ) = @_;
 648
 649     # Returns an arrayref.
 650
 651     my ( $entry, $clones, $score, @array );
 652
 653     while ( $entry = bed_entry_get( $fh ) )
 654     {
 655         if ( $entry->[ name ] =~ /_(\d+)$/ )
 656         {
 657             $clones = $1;
 658
 659             if ( $entry->[ score ] )
 660             {
 661                 $score = int( $clones / $entry->[ score ] );
 662
 663                 map { $array[ $_ ] += $score } $entry->[ chromStart ] .. $entry->[ chromEnd ] - 1 if $score >= 1;
 664             }
 665         }
 666     }
 667
 668     return wantarray ? @array : \@array;
 669 }
 670
 671
 672 sub tag_contigs_scan
 673 {
 674     # Martin A. Hansen, November 2008.
 675
 676     # Scans an array with tag contigs and locates
 677     # the next contig from a given position. The
 678     # score of the tag contig is determined as the
 679     # maximum value of the tag contig. If a tag contig
 680     # is found a triple is returned with beg, end and score
 681     # otherwise an empty list is returned.
 682
 683     my ( $array,   # array to scan
 684          $beg,     # position to start scanning from
 685        ) = @_;
 686
 687     # Returns an arrayref.
 688
 689     my ( $end, $score );
 690
 691     $score = 0;
 692
 693     while ( $beg < scalar @{ $array } and not $array->[ $beg ] ) { $beg++ }
 694
 695     $end = $beg;
 696
 697     while ( $array->[ $end ] )
 698     {
 699         $score = Maasha::Calc::max( $score, $array->[ $end ] );
 700
 701         $end++
 702     }
 703
 704     if ( $score > 0 ) {
 705         return wantarray ? ( $beg, $end, $score ) : [ $beg, $end, $score ];
 706     } else {
 707         return wantarray ? () : [];
 708     }
 709 }
 710
 711
 712 sub bed_upload_to_ucsc
 713 {
 714     # Martin A. Hansen, September 2007.
 715
 716     # Upload a BED file to the UCSC database.
 717
 718     my ( $tmp_dir,   # temporary directory
 719          $file,      # file to upload,
 720          $options,   # argument hashref
 721          $append,    # flag indicating table should be appended
 722        ) = @_;
 723
 724     # Returns nothing.
 725
 726     my ( $args, $table, $sql_file, $fh_out, $fh_in );
 727
 728     if ( $append ) {
 729         $args = join " ", $options->{ "database" }, $options->{ "table" }, "-tmpDir=$tmp_dir", "-oldTable", $file;
 730     } else {
 731         $args = join " ", $options->{ "database" }, $options->{ "table" }, "-tmpDir=$tmp_dir", $file;
 732     }
 733
 734     if ( $options->{ "table" } =~ /rnaSecStr/ )
 735     {
 736         $table = $options->{ "table" };
 737
 738         print qq(uploading secondary structure table:"$table"\n) if $options->{ "verbose" };
 739
 740         $sql_file = "$tmp_dir/upload_RNA_SS.sql";
 741
 742         $fh_out   = Maasha::Filesys::file_write_open( $sql_file );
 743
 744         print $fh_out qq(
 745 CREATE TABLE $table (
 746     bin smallint not null,              # Bin number for browser speedup
 747     chrom varchar(255) not null,        # Chromosome or FPC contig
 748     chromStart int unsigned not null,   # Start position in chromosome
 749     chromEnd int unsigned not null,     # End position in chromosome
 750     name varchar(255) not null,         # Name of item
 751     score int unsigned not null,        # Score from 0-1000
 752     strand char(1) not null,            # + or -
 753     size int unsigned not null,         # Size of element.
 754     secStr longblob not null,           # Parentheses and '.'s which define the secondary structure
 755     conf longblob not null,             # Confidence of secondary-structure annotation per position (0.0-1.0).
 756     #Indices
 757     INDEX(name(16)),
 758     INDEX(chrom(8), bin),
 759     INDEX(chrom(8), chromStart)
 760 );
 761         );
 762
 763         close $fh_out;
 764
 765         Maasha::Common::run( "hgLoadBed", "-notItemRgb -sqlTable=$sql_file $options->{ 'database' } $options->{ 'table' } -tmpDir=$tmp_dir $file > /dev/null 2>&1" );
 766
 767         unlink $sql_file;
 768     }
 769     else
 770     {
 771         Maasha::Common::run( "hgLoadBed", "$args > /dev/null 2>&1" );
 772     }
 773 }
 774
 775
 776 sub bed_analyze
 777 {
 778     # Martin A. Hansen, March 2008.
 779
 780     # Given a bed record, analysis this to give information
 781     # about intron/exon sizes.
 782
 783     my ( $entry,   # BED entry
 784        ) = @_;
 785
 786     # Returns hashref.
 787
 788     my ( $i, @begs, @lens, $exon_max, $exon_min, $exon_len, $exon_tot, $intron_max, $intron_min, $intron_len, $intron_tot );
 789
 790     $exon_max   = 0;
 791     $exon_min   = 9999999999;
 792     $intron_max = 0;
 793     $intron_min = 9999999999;
 794
 795     $entry->{ "EXONS" }   = $entry->{ "BLOCK_COUNT" };
 796
 797     @begs = split /,/, $entry->{ "Q_BEGS" };
 798     @lens = split /,/, $entry->{ "BLOCK_LENS" };
 799
 800     for ( $i = 0; $i < $entry->{ "BLOCK_COUNT" }; $i++ )
 801     {
 802         $exon_len = $lens[ $i ];
 803
 804         $entry->{ "EXON_LEN_$i" } = $exon_len;
 805
 806         $exon_max = $exon_len if $exon_len > $exon_max;
 807         $exon_min = $exon_len if $exon_len < $exon_min;
 808
 809         $exon_tot += $exon_len;
 810     }
 811
 812     $entry->{ "EXON_LEN_-1" }   = $exon_len;
 813     $entry->{ "EXON_MAX_LEN" }  = $exon_max;
 814     $entry->{ "EXON_MIN_LEN" }  = $exon_min;
 815     $entry->{ "EXON_MEAN_LEN" } = int( $exon_tot / $entry->{ "EXONS" } );
 816
 817     $entry->{ "INTRONS" } = $entry->{ "BLOCK_COUNT" } - 1;
 818     $entry->{ "INTRONS" } = 0 if $entry->{ "INTRONS" } < 0;
 819
 820     if ( $entry->{ "INTRONS" } )
 821     {
 822         for ( $i = 1; $i < $entry->{ "BLOCK_COUNT" }; $i++ )
 823         {
 824             $intron_len = $begs[ $i ] - ( $begs[ $i - 1 ] + $lens[ $i - 1 ] );
 825
 826             $entry->{ "INTRON_LEN_" . ( $i - 1 ) } = $intron_len;
 827
 828             $intron_max = $intron_len if $intron_len > $intron_max;
 829             $intron_min = $intron_len if $intron_len < $intron_min;
 830
 831             $intron_tot += $intron_len;
 832         }
 833
 834         $entry->{ "INTRON_LEN_-1" }   = $intron_len;
 835         $entry->{ "INTRON_MAX_LEN" }  = $intron_max;
 836         $entry->{ "INTRON_MIN_LEN" }  = $intron_min;
 837         $entry->{ "INTRON_MEAN_LEN" } = int( $intron_tot / $entry->{ "INTRONS" } );
 838     }
 839
 840     return wantarray ? %{ $entry } : $entry;
 841 }
 842
 843
 844 sub bed_merge_entries
 845 {
 846     # Martin A. Hansen, February 2008.
 847
 848     # Merge a list of given BED entries in one big entry.
 849
 850     my ( $entries,     # list of BED entries to be merged
 851        ) = @_;
 852
 853     # Returns hash.
 854
 855     my ( $i, @q_ids, @q_begs, @blocksizes, @new_q_begs, @new_blocksizes, %new_entry );
 856
 857     @{ $entries } = sort { $a->{ "CHR_BEG" } <=> $b->{ "CHR_BEG" } } @{ $entries };
 858
 859     for ( $i = 0; $i < @{ $entries }; $i++ )
 860     {
 861         Maasha::Common::error( qq(Attempted merge of BED entries from different chromosomes) ) if $entries->[ 0 ]->{ "CHR" }    ne $entries->[ $i ]->{ "CHR" };
 862         Maasha::Common::error( qq(Attempted merge of BED entries from different strands) )     if $entries->[ 0 ]->{ "STRAND" } ne $entries->[ $i ]->{ "STRAND" };
 863
 864         push @q_ids, $entries->[ $i ]->{ "Q_ID" } || sprintf( "ID%06d", $i );
 865
 866         if ( exists $entries->[ $i ]->{ "Q_BEGS" } )
 867         {
 868             @q_begs     = split ",", $entries->[ $i ]->{ "Q_BEGS" };
 869             @blocksizes = split ",", $entries->[ $i ]->{ "BLOCK_LENS" };
 870         }
 871         else
 872         {
 873             @q_begs     = 0;
 874             @blocksizes = $entries->[ $i ]->{ "CHR_END" } - $entries->[ $i ]->{ "CHR_BEG" } + 1;
 875         }
 876
 877         map { $_ += $entries->[ $i ]->{ "CHR_BEG" } } @q_begs;
 878
 879         push @new_q_begs, @q_begs;
 880         push @new_blocksizes, @blocksizes;
 881     }
 882
 883     map { $_ -= $entries->[ 0 ]->{ "CHR_BEG" } } @new_q_begs;
 884
 885     %new_entry = (
 886         CHR          => $entries->[ 0 ]->{ "CHR" },
 887         CHR_BEG      => $entries->[ 0 ]->{ "CHR_BEG" },
 888         CHR_END      => $entries->[ -1 ]->{ "CHR_END" },
 889         REC_TYPE     => "BED",
 890         BED_LEN      => $entries->[ -1 ]->{ "CHR_END" } - $entries->[ 0 ]->{ "CHR_BEG" } + 1,
 891         BED_COLS     => 12,
 892         Q_ID         => join( ":", @q_ids ),
 893         SCORE        => 999,
 894         STRAND       => $entries->[ 0 ]->{ "STRAND" }     || "+",
 895         THICK_BEG    => $entries->[ 0 ]->{ "THICK_BEG" }  || $entries->[ 0 ]->{ "CHR_BEG" },
 896         THICK_END    => $entries->[ -1 ]->{ "THICK_END" } || $entries->[ -1 ]->{ "CHR_END" },
 897         COLOR        => 0,
 898         BLOCK_COUNT  => scalar @new_q_begs,
 899         BLOCK_LENS   => join( ",", @new_blocksizes ),
 900         Q_BEGS       => join( ",", @new_q_begs ),
 901     );
 902
 903     return wantarray ? %new_entry : \%new_entry;
 904 }
 905
 906
 907 sub bed_split_entry
 908 {
 909     # Martin A. Hansen, February 2008.
 910
 911     # Splits a given BED entry into a list of blocks,
 912     # which are returned. A list of 6 column BED entry is returned.
 913
 914     my ( $entry,    # BED entry hashref
 915        ) = @_;
 916
 917     # Returns a list.
 918
 919     my ( @q_begs, @blocksizes, $block, @blocks, $i );
 920
 921     if ( exists $entry->{ "BLOCK_COUNT" } )
 922     {
 923         @q_begs     = split ",", $entry->{ "Q_BEGS" };
 924         @blocksizes = split ",", $entry->{ "BLOCK_LENS" };
 925
 926         for ( $i = 0; $i < @q_begs; $i++ )
 927         {
 928             undef $block;
 929
 930             $block->{ "CHR" }      = $entry->{ "CHR" };
 931             $block->{ "CHR_BEG" }  = $entry->{ "CHR_BEG" } + $q_begs[ $i ];
 932             $block->{ "CHR_END" }  = $entry->{ "CHR_BEG" } + $q_begs[ $i ] + $blocksizes[ $i ] - 1;
 933             $block->{ "Q_ID" }     = $entry->{ "Q_ID" } . sprintf( "_%03d", $i );
 934             $block->{ "SCORE" }    = $entry->{ "SCORE" };
 935             $block->{ "STRAND" }   = $entry->{ "STRAND" };
 936             $block->{ "BED_LEN" }  = $block->{ "CHR_END" } - $block->{ "CHR_BEG" } + 1,
 937             $block->{ "BED_COLS" } = 6;
 938             $block->{ "REC_TYPE" } = "BED";
 939
 940             push @blocks, $block;
 941         }
 942     }
 943     else
 944     {
 945         @blocks = @{ $entry };
 946     }
 947
 948     return wantarray ? @blocks : \@blocks;
 949 }
 950
 951
 952 sub bed_overlap
 953 {
 954     # Martin A. Hansen, February 2008.
 955
 956     # Checks if two BED entries overlap and
 957     # return 1 if so - else 0;
 958
 959     my ( $entry1,      # hashref
 960          $entry2,      # hashref
 961          $no_strand,   # don't check strand flag - OPTIONAL
 962        ) = @_;
 963
 964     # Return bolean.
 965
 966     return 0 if $entry1->{ "CHR" }    ne $entry2->{ "CHR" };
 967     return 0 if $entry1->{ "STRAND" } ne $entry2->{ "STRAND" };
 968
 969     if ( $entry1->{ "CHR_END" } < $entry2->{ "CHR_BEG" } or $entry1->{ "CHR_BEG" } > $entry2->{ "CHR_END" } ) {
 970         return 0;
 971     } else {
 972         return 1;
 973     }
 974 }
 975
 976
 977
 978
 979 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 980
 981
 982 1;
 983
 984
 985 __END__