code_perl/Maasha/UCSC.pm

   1 package Maasha::UCSC;
   2
   3 # Copyright (C) 2007 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 # Stuff for interacting with UCSC genome browser
  26
  27
  28 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  29
  30
  31 use strict;
  32 use vars qw ( @ISA @EXPORT );
  33
  34 use Data::Dumper;
  35 use Time::HiRes qw( gettimeofday );
  36
  37 use Maasha::Common;
  38 use Maasha::Calc;
  39 use Maasha::Matrix;
  40
  41 use constant {
  42     CHR_BEG      => 0,
  43     NEXT_CHR_BEG => 1,
  44     CHR_END      => 2,
  45     INDEX_BEG    => 3,
  46     INDEX_LEN    => 4,
  47 };
  48
  49 @ISA = qw( Exporter );
  50
  51 my $TIME = gettimeofday();
  52
  53
  54 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> BED format <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  55
  56
  57 # http://genome.ucsc.edu/goldenPath/help/hgTracksHelp.html#BED
  58
  59
  60 sub bed_get_entry
  61 {
  62     # Martin A. Hansen, December 2007.
  63
  64     # Reads a bed entry given a filehandle.
  65
  66     my ( $fh,        # file handle
  67          $columns,   # number of BED columns to read  -  OPTIONAL
  68        ) = @_;
  69
  70     # Returns hashref.
  71
  72     my ( $line, @fields, %entry );
  73
  74     $line = <$fh>;
  75
  76     $line =~ tr/\n\r//d;    # some people have carriage returns in their BED files -> Grrrr
  77
  78     return if not defined $line;
  79
  80     @fields = split "\t", $line;
  81
  82     $columns ||= scalar @fields;
  83
  84     if ( $columns == 3 )
  85     {
  86         %entry = (
  87             "CHR"      => $fields[ 0 ],
  88             "CHR_BEG"  => $fields[ 1 ],
  89             "CHR_END"  => $fields[ 2 ] - 1,
  90         );
  91     }
  92     elsif ( $columns == 4 )
  93     {
  94         %entry = (
  95             "CHR"      => $fields[ 0 ],
  96             "CHR_BEG"  => $fields[ 1 ],
  97             "CHR_END"  => $fields[ 2 ] - 1,
  98             "Q_ID"     => $fields[ 3 ],
  99         );
 100     }
 101     elsif ( $columns == 5 )
 102     {
 103         %entry = (
 104             "CHR"      => $fields[ 0 ],
 105             "CHR_BEG"  => $fields[ 1 ],
 106             "CHR_END"  => $fields[ 2 ] - 1,
 107             "Q_ID"     => $fields[ 3 ],
 108             "SCORE"    => $fields[ 4 ],
 109         );
 110     }
 111     elsif ( $columns == 6 )
 112     {
 113         %entry = (
 114             "CHR"      => $fields[ 0 ],
 115             "CHR_BEG"  => $fields[ 1 ],
 116             "CHR_END"  => $fields[ 2 ] - 1,
 117             "Q_ID"     => $fields[ 3 ],
 118             "SCORE"    => $fields[ 4 ],
 119             "STRAND"   => $fields[ 5 ],
 120         );
 121     }
 122     elsif ( $columns == 12 )
 123     {
 124         %entry = (
 125             "CHR"        => $fields[ 0 ],
 126             "CHR_BEG"    => $fields[ 1 ],
 127             "CHR_END"    => $fields[ 2 ] - 1,
 128             "Q_ID"       => $fields[ 3 ],
 129             "SCORE"      => $fields[ 4 ],
 130             "STRAND"     => $fields[ 5 ],
 131             "THICK_BEG"  => $fields[ 6 ],
 132             "THICK_END"  => $fields[ 7 ] - 1,
 133             "ITEMRGB"    => $fields[ 8 ],
 134             "BLOCKCOUNT" => $fields[ 9 ],
 135             "BLOCKSIZES" => $fields[ 10 ],
 136             "Q_BEGS"     => $fields[ 11 ],
 137         );
 138     }
 139     else
 140     {
 141         Maasha::Common::error( qq(Bad BED format in line->$line<-) );
 142     }
 143
 144     $entry{ "REC_TYPE" } = "BED";
 145     $entry{ "BED_LEN" }  = $entry{ "CHR_END" } - $entry{ "CHR_BEG" } + 1;
 146     $entry{ "BED_COLS" } = $columns;
 147
 148     return wantarray ? %entry : \%entry;
 149 }
 150
 151
 152 sub bed_get_entries
 153 {
 154     # Martin A. Hansen, January 2008.
 155
 156     # Given a path to a BED file, read in all entries
 157     # and return.
 158
 159     my ( $path,     # full path to BED file
 160          $columns,  # number of columns in BED file - OPTIONAL (but is faster)
 161        ) = @_;
 162
 163     # Returns a list.
 164
 165     my ( $fh, $entry, @list );
 166
 167     $fh = Maasha::Common::read_open( $path );
 168
 169     while ( $entry = bed_get_entry( $fh ) ) {
 170         push @list, $entry;
 171     }
 172
 173     close $fh;
 174
 175     return wantarray ? @list : \@list;
 176 }
 177
 178
 179 sub bed_put_entry
 180 {
 181     # Martin A. Hansen, Septermber 2007.
 182
 183     # Writes a BED entry to file.
 184
 185     # NB, this could really be more robust!?
 186
 187     my ( $record,       # hashref
 188          $fh,           # file handle                   - OPTIONAL
 189          $columns,      # number of columns in BED file - OPTIONAL (but is faster)
 190        ) = @_;
 191
 192     # Returns nothing.
 193
 194     my ( @fields );
 195
 196     $columns ||= 12;   # max number of columns possible
 197
 198     if ( $columns == 3 )
 199     {
 200         push @fields, $record->{ "CHR" };
 201         push @fields, $record->{ "CHR_BEG" };
 202         push @fields, $record->{ "CHR_END" } + 1;
 203     }
 204     elsif ( $columns == 4 )
 205     {
 206         $record->{ "Q_ID" }  =~ s/\s+/_/g;
 207
 208         push @fields, $record->{ "CHR" };
 209         push @fields, $record->{ "CHR_BEG" };
 210         push @fields, $record->{ "CHR_END" } + 1;
 211         push @fields, $record->{ "Q_ID" };
 212     }
 213     elsif ( $columns == 5 )
 214     {
 215         $record->{ "Q_ID" }  =~ s/\s+/_/g;
 216         $record->{ "SCORE" } =~ s/\.\d*//;
 217
 218         push @fields, $record->{ "CHR" };
 219         push @fields, $record->{ "CHR_BEG" };
 220         push @fields, $record->{ "CHR_END" } + 1;
 221         push @fields, $record->{ "Q_ID" };
 222         push @fields, $record->{ "SCORE" };
 223     }
 224     elsif ( $columns == 6 )
 225     {
 226         $record->{ "Q_ID" }  =~ s/\s+/_/g;
 227         $record->{ "SCORE" } =~ s/\.\d*//;
 228
 229         push @fields, $record->{ "CHR" };
 230         push @fields, $record->{ "CHR_BEG" };
 231         push @fields, $record->{ "CHR_END" } + 1;
 232         push @fields, $record->{ "Q_ID" };
 233         push @fields, $record->{ "SCORE" };
 234         push @fields, $record->{ "STRAND" };
 235     }
 236     else
 237     {
 238         $record->{ "Q_ID" }  =~ s/\s+/_/g;
 239         $record->{ "SCORE" } =~ s/\.\d*//;
 240
 241         push @fields, $record->{ "CHR" };
 242         push @fields, $record->{ "CHR_BEG" };
 243         push @fields, $record->{ "CHR_END" } + 1;
 244         push @fields, $record->{ "Q_ID" };
 245         push @fields, $record->{ "SCORE" };
 246         push @fields, $record->{ "STRAND" };
 247         push @fields, $record->{ "THICK_BEG" }     if defined $record->{ "THICK_BEG" };
 248         push @fields, $record->{ "THICK_END" } + 1 if defined $record->{ "THICK_END" };
 249         push @fields, $record->{ "ITEMRGB" }       if defined $record->{ "ITEMRGB" };
 250         push @fields, $record->{ "BLOCKCOUNT" }    if defined $record->{ "BLOCKCOUNT" };
 251         push @fields, $record->{ "BLOCKSIZES" }    if defined $record->{ "BLOCKSIZES" };
 252         push @fields, $record->{ "Q_BEGS" }        if defined $record->{ "Q_BEGS" };
 253     }
 254
 255     if ( $fh ) {
 256         print $fh join( "\t", @fields ), "\n";
 257     } else {
 258         print join( "\t", @fields ), "\n";
 259     }
 260 }
 261
 262
 263 sub bed_put_entries
 264 {
 265     # Martin A. Hansen, January 2008.
 266
 267     # Write a list of BED entries.
 268
 269     my ( $entries,   # list of entries,
 270          $fh,        # file handle - OPTIONAL
 271        ) = @_;
 272
 273     # Returns nothing.
 274
 275     map { bed_put_entry( $_, $fh ) } @{ $entries };
 276 }
 277
 278
 279 sub bed_analyze
 280 {
 281     # Martin A. Hansen, March 2008.
 282
 283     # Given a bed record, analysis this to give information
 284     # about intron/exon sizes.
 285
 286     my ( $entry,   # BED entry
 287        ) = @_;
 288
 289     # Returns hashref.
 290
 291     my ( $i, @begs, @lens, $exon_max, $exon_min, $exon_len, $exon_tot, $intron_max, $intron_min, $intron_len, $intron_tot );
 292
 293     $exon_max   = 0;
 294     $exon_min   = 9999999999;
 295     $intron_max = 0;
 296     $intron_min = 9999999999;
 297
 298     $entry->{ "EXONS" }   = $entry->{ "BLOCKCOUNT" };
 299
 300     @begs = split /,/, $entry->{ "Q_BEGS" };
 301     @lens = split /,/, $entry->{ "BLOCKSIZES" };
 302
 303     for ( $i = 0; $i < $entry->{ "BLOCKCOUNT" }; $i++ )
 304     {
 305         $exon_len = @lens[ $i ];
 306
 307         $entry->{ "EXON_LEN_$i" } = $exon_len;
 308
 309         $exon_max = $exon_len if $exon_len > $exon_max;
 310         $exon_min = $exon_len if $exon_len < $exon_min;
 311
 312         $exon_tot += $exon_len;
 313     }
 314
 315     $entry->{ "EXON_LEN_-1" }   = $exon_len;
 316     $entry->{ "EXON_MAX_LEN" }  = $exon_max;
 317     $entry->{ "EXON_MIN_LEN" }  = $exon_min;
 318     $entry->{ "EXON_MEAN_LEN" } = int( $exon_tot / $entry->{ "EXONS" } );
 319
 320     $entry->{ "INTRONS" } = $entry->{ "BLOCKCOUNT" } - 1;
 321     $entry->{ "INTRONS" } = 0 if $entry->{ "INTRONS" } < 0;
 322
 323     if ( $entry->{ "INTRONS" } )
 324     {
 325         for ( $i = 1; $i < $entry->{ "BLOCKCOUNT" }; $i++ )
 326         {
 327             $intron_len = @begs[ $i ] - ( @begs[ $i - 1 ] + @lens[ $i - 1 ] );
 328
 329             $entry->{ "INTRON_LEN_" . ( $i - 1 ) } = $intron_len;
 330
 331             $intron_max = $intron_len if $intron_len > $intron_max;
 332             $intron_min = $intron_len if $intron_len < $intron_min;
 333
 334             $intron_tot += $intron_len;
 335         }
 336
 337         $entry->{ "INTRON_LEN_-1" }   = $intron_len;
 338         $entry->{ "INTRON_MAX_LEN" }  = $intron_max;
 339         $entry->{ "INTRON_MIN_LEN" }  = $intron_min;
 340         $entry->{ "INTRON_MEAN_LEN" } = int( $intron_tot / $entry->{ "INTRONS" } );
 341     }
 342
 343     return wantarray ? %{ $entry } : $entry;
 344 }
 345
 346
 347 sub bed_sort
 348 {
 349     # Martin A. Hansen, March 2008.
 350
 351     # Sort a potential huge BED file according to
 352     # CHR, CHR_BEG and optionally STRAND.
 353
 354     my ( $tmp_dir,   # temporary directory used for sorting
 355          $file,      # BED file to sort
 356          $strand,    # flag to sort on strand - OPTIONAL
 357        ) = @_;
 358
 359     # Returns nothing.
 360
 361     my ( $fh_in, $key, $fh_out, %fh_hash, $part_file, $entry, $entries );
 362
 363     $fh_in = Maasha::Common::read_open( $file );
 364
 365     while ( $entry = bed_get_entry( $fh_in ) )
 366     {
 367         if ( $strand ) {
 368             $key = join "_", $entry->{ "CHR" }, $entry->{ "STRAND" };
 369         } else {
 370             $key = $entry->{ "CHR" };
 371         }
 372
 373         $fh_hash{ $key } = Maasha::Common::write_open( "$tmp_dir/$key.sort" ) if not exists $fh_hash{ $key };
 374
 375         bed_put_entry( $entry, $fh_hash{ $key } );
 376     }
 377
 378     close $fh_in;
 379
 380     map { close $_ } keys %fh_hash;
 381
 382     $fh_out = Maasha::Common::write_open( "$tmp_dir/temp.sort" );
 383
 384     foreach $part_file ( sort keys %fh_hash )
 385     {
 386         $entries = bed_get_entries( "$tmp_dir/$part_file.sort" );
 387
 388         @{ $entries } = sort { $a->{ "CHR_BEG" } <=> $b->{ "CHR_BEG" } } @{ $entries };
 389
 390         map { bed_put_entry( $_, $fh_out ) } @{ $entries };
 391
 392         unlink "$tmp_dir/$part_file.sort";
 393     }
 394
 395     close $fh_out;
 396
 397     rename "$tmp_dir/temp.sort", $file;
 398 }
 399
 400
 401 sub bed_merge_entries
 402 {
 403     # Martin A. Hansen, February 2008.
 404
 405     # Merge a list of given BED entries in one big entry.
 406
 407     my ( $entries,     # list of BED entries to be merged
 408        ) = @_;
 409
 410     # Returns hash.
 411
 412     my ( $i, @q_ids, @q_begs, @blocksizes, @new_q_begs, @new_blocksizes, %new_entry );
 413
 414     @{ $entries } = sort { $a->{ "CHR_BEG" } <=> $b->{ "CHR_BEG" } } @{ $entries };
 415
 416     for ( $i = 0; $i < @{ $entries }; $i++ )
 417     {
 418         Maasha::Common::error( qq(Attempted merge of BED entries from different chromosomes) ) if $entries->[ 0 ]->{ "CHR" }    ne $entries->[ $i ]->{ "CHR" };
 419         Maasha::Common::error( qq(Attempted merge of BED entries from different strands) )     if $entries->[ 0 ]->{ "STRAND" } ne $entries->[ $i ]->{ "STRAND" };
 420
 421         push @q_ids, $entries->[ $i ]->{ "Q_ID" } || sprintf( "ID%06d", $i );
 422
 423         if ( exists $entries->[ $i ]->{ "Q_BEGS" } )
 424         {
 425             @q_begs     = split ",", $entries->[ $i ]->{ "Q_BEGS" };
 426             @blocksizes = split ",", $entries->[ $i ]->{ "BLOCKSIZES" };
 427         }
 428         else
 429         {
 430             @q_begs     = 0;
 431             @blocksizes = $entries->[ $i ]->{ "CHR_END" } - $entries->[ $i ]->{ "CHR_BEG" } + 1;
 432         }
 433
 434         map { $_ += $entries->[ $i ]->{ "CHR_BEG" } } @q_begs;
 435
 436         push @new_q_begs, @q_begs;
 437         push @new_blocksizes, @blocksizes;
 438     }
 439
 440     map { $_ -= $entries->[ 0 ]->{ "CHR_BEG" } } @new_q_begs;
 441
 442     %new_entry = (
 443         CHR         => $entries->[ 0 ]->{ "CHR" },
 444         CHR_BEG     => $entries->[ 0 ]->{ "CHR_BEG" },
 445         CHR_END     => $entries->[ -1 ]->{ "CHR_END" },
 446         REC_TYPE    => "BED",
 447         BED_LEN     => $entries->[ -1 ]->{ "CHR_END" } - $entries->[ 0 ]->{ "CHR_BEG" } + 1,
 448         BED_COLS    => 12,
 449         Q_ID        => join( ":", @q_ids ),
 450         SCORE       => 999,
 451         STRAND      => $entries->[ 0 ]->{ "STRAND" }     || "+",
 452         THICK_BEG   => $entries->[ 0 ]->{ "THICK_BEG" }  || $entries->[ 0 ]->{ "CHR_BEG" },
 453         THICK_END   => $entries->[ -1 ]->{ "THICK_END" } || $entries->[ -1 ]->{ "CHR_END" },
 454         ITEMRGB     => "0,0,0",
 455         BLOCKCOUNT  => scalar @new_q_begs,
 456         BLOCKSIZES  => join( ",", @new_blocksizes ),
 457         Q_BEGS      => join( ",", @new_q_begs ),
 458     );
 459
 460     return wantarray ? %new_entry : \%new_entry;
 461 }
 462
 463
 464 sub bed_split_entry
 465 {
 466     # Martin A. Hansen, February 2008.
 467
 468     # Splits a given BED entry into a list of blocks,
 469     # which are returned. A list of 6 column BED entry is returned.
 470
 471     my ( $entry,    # BED entry hashref
 472        ) = @_;
 473
 474     # Returns a list.
 475
 476     my ( @q_begs, @blocksizes, $block, @blocks, $i );
 477
 478     if ( exists $entry->{ "BLOCKCOUNT" } )
 479     {
 480         @q_begs     = split ",", $entry->{ "Q_BEGS" };
 481         @blocksizes = split ",", $entry->{ "BLOCKSIZES" };
 482
 483         for ( $i = 0; $i < @q_begs; $i++ )
 484         {
 485             undef $block;
 486
 487             $block->{ "CHR" }      = $entry->{ "CHR" };
 488             $block->{ "CHR_BEG" }  = $entry->{ "CHR_BEG" } + $q_begs[ $i ];
 489             $block->{ "CHR_END" }  = $entry->{ "CHR_BEG" } + $q_begs[ $i ] + $blocksizes[ $i ] - 1;
 490             $block->{ "Q_ID" }     = $entry->{ "Q_ID" } . sprintf( "_%03d", $i );
 491             $block->{ "SCORE" }    = $entry->{ "SCORE" };
 492             $block->{ "STRAND" }   = $entry->{ "STRAND" };
 493             $block->{ "BED_LEN" }  = $block->{ "CHR_END" } - $block->{ "CHR_BEG" } + 1,
 494             $block->{ "BED_COLS" } = 6;
 495             $block->{ "REC_TYPE" } = "BED";
 496
 497             push @blocks, $block;
 498         }
 499     }
 500     else
 501     {
 502         @blocks = @{ $entry };
 503     }
 504
 505     return wantarray ? @blocks : \@blocks;
 506 }
 507
 508
 509
 510 sub bed_overlap
 511 {
 512     # Martin A. Hansen, February 2008.
 513
 514     # Checks if two BED entries overlap and
 515     # return 1 if so - else 0;
 516
 517     my ( $entry1,      # hashref
 518          $entry2,      # hashref
 519          $no_strand,   # don't check strand flag - OPTIONAL
 520        ) = @_;
 521
 522     # Return bolean.
 523
 524     return 0 if $entry1->{ "CHR" }    ne $entry2->{ "CHR" };
 525     return 0 if $entry1->{ "STRAND" } ne $entry2->{ "STRAND" };
 526
 527     if ( $entry1->{ "CHR_END" } < $entry2->{ "CHR_BEG" } or $entry1->{ "CHR_BEG" } > $entry2->{ "CHR_END" } ) {
 528         return 0;
 529     } else {
 530         return 1;
 531     }
 532 }
 533
 534
 535 sub bed_upload_to_ucsc
 536 {
 537     # Martin A. Hansen, September 2007.
 538
 539     # Upload a BED file to the UCSC database.
 540
 541     my ( $tmp_dir,   # temporary directory
 542          $file,      # file to upload,
 543          $options,   # argument hashref
 544          $append,    # flag indicating table should be appended
 545        ) = @_;
 546
 547     # Returns nothing.
 548
 549     my ( $args, $table, $sql_file, $fh_out, $fh_in );
 550
 551     if ( $append ) {
 552         $args = join " ", $options->{ "database" }, $options->{ "table" }, "-tmpDir=$tmp_dir", "-oldTable", $file;
 553     } else {
 554         $args = join " ", $options->{ "database" }, $options->{ "table" }, "-tmpDir=$tmp_dir", $file;
 555     }
 556
 557     if ( $options->{ "sec_struct" } )
 558     {
 559         $table = $options->{ "table" };
 560
 561         Maasha::Common::error( "Attempt to load secondary structure track without 'rnaSecStr' in table name" ) if not $table =~ /rnaSecStr/;
 562
 563         $sql_file = "$tmp_dir/upload_RNA_SS.sql";
 564
 565         $fh_out   = Maasha::Common::write_open( $sql_file );
 566
 567         print $fh_out qq(
 568 CREATE TABLE $table (
 569     bin smallint not null,              # Bin number for browser speedup
 570     chrom varchar(255) not null,        # Chromosome or FPC contig
 571     chromStart int unsigned not null,   # Start position in chromosome
 572     chromEnd int unsigned not null,     # End position in chromosome
 573     name varchar(255) not null,         # Name of item
 574     score int unsigned not null,        # Score from 0-1000
 575     strand char(1) not null,            # + or -
 576     size int unsigned not null,         # Size of element.
 577     secStr longblob not null,           # Parentheses and '.'s which define the secondary structure
 578     conf longblob not null,             # Confidence of secondary-structure annotation per position (0.0-1.0).
 579     #Indices
 580     INDEX(name(16)),
 581     INDEX(chrom(8), bin),
 582     INDEX(chrom(8), chromStart)
 583 );
 584         );
 585
 586         close $fh_out;
 587
 588         Maasha::Common::run( "hgLoadBed", "-notItemRgb -sqlTable=$sql_file $options->{ 'database' } $options->{ 'table' } -tmpDir=$tmp_dir $file > /dev/null 2>&1" );
 589
 590         unlink $sql_file;
 591     }
 592     else
 593     {
 594         Maasha::Common::run( "hgLoadBed", "$args > /dev/null 2>&1" );
 595     }
 596 }
 597
 598
 599 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> PSL format <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 600
 601
 602 sub psl_get_entries
 603 {
 604     # Martin A. Hansen, February 2008.
 605
 606     # Reads PSL entries and returns a record.
 607
 608     my ( $path,   # full path to PSL file
 609        ) = @_;
 610
 611     # Returns hashref.
 612
 613     my ( $fh, @lines, @fields, $i, %record, @records );
 614
 615     $fh = Maasha::Common::read_open( $path );
 616
 617     @lines = <$fh>;
 618
 619     close $fh;
 620
 621     chomp @lines;
 622
 623     for ( $i = 5; $i < @lines; $i++ )
 624     {
 625         @fields = split "\t", $lines[ $i ];
 626
 627         Maasha::Common::error( qq(Bad PSL format in file "$path") ) if not @fields == 21;
 628
 629         undef %record;
 630
 631         %record = (
 632             REC_TYPE    => "PSL",
 633             MATCHES     => $fields[ 0 ],
 634             MISMATCHES  => $fields[ 1 ],
 635             REPMATCHES  => $fields[ 2 ],
 636             NCOUNT      => $fields[ 3 ],
 637             QNUMINSERT  => $fields[ 4 ],
 638             QBASEINSERT => $fields[ 5 ],
 639             SNUMINSERT  => $fields[ 6 ],
 640             SBASEINSERT => $fields[ 7 ],
 641             STRAND      => $fields[ 8 ],
 642             Q_ID        => $fields[ 9 ],
 643             Q_LEN       => $fields[ 10 ],
 644             Q_BEG       => $fields[ 11 ],
 645             Q_END       => $fields[ 12 ] - 1,
 646             S_ID        => $fields[ 13 ],
 647             S_LEN       => $fields[ 14 ],
 648             S_BEG       => $fields[ 15 ],
 649             S_END       => $fields[ 16 ] - 1,
 650             BLOCKCOUNT  => $fields[ 17 ],
 651             BLOCKSIZES  => $fields[ 18 ],
 652             Q_BEGS      => $fields[ 19 ],
 653             S_BEGS      => $fields[ 20 ],
 654         );
 655
 656         $record{ "SCORE" } = $record{ "MATCHES" } + int( $record{ "REPMATCHES" } / 2 ) - $record{ "MISMATCHES" } - $record{ "QNUMINSERT" } - $record{ "SNUMINSERT" };
 657
 658         push @records, { %record };
 659     }
 660
 661     return wantarray ? @records : \@records;
 662 }
 663
 664
 665 sub psl_put_header
 666 {
 667     # Martin A. Hansen, September 2007.
 668
 669     # Write a PSL header to file.
 670
 671     my ( $fh,  # file handle  - OPTIONAL
 672        ) = @_;
 673
 674     # Returns nothing.
 675
 676     $fh = \*STDOUT if not $fh;
 677
 678     print $fh qq(psLayout version 3
 679 match   mis-    rep.    N's     Q gap   Q gap   T gap   T gap   strand  Q               Q       Q       Q       T               T       T       T       block   blockSizes      qStart        match   match           count   bases   count   bases           name            size    start   end     name            size    start   end     count
 680 ---------------------------------------------------------------------------------------------------------------------------------------------------------------
 681 );
 682 }
 683
 684
 685 sub psl_put_entry
 686 {
 687     # Martin A. Hansen, September 2007.
 688
 689     # Write a PSL entry to file.
 690
 691     my ( $record,       # hashref
 692          $fh,           # file handle  -  OPTIONAL
 693        ) = @_;
 694
 695     # Returns nothing.
 696
 697     $fh = \*STDOUT if not $fh;
 698
 699     my @output;
 700
 701     push @output, $record->{ "MATCHES" };
 702     push @output, $record->{ "MISMATCHES" };
 703     push @output, $record->{ "REPMATCHES" };
 704     push @output, $record->{ "NCOUNT" };
 705     push @output, $record->{ "QNUMINSERT" };
 706     push @output, $record->{ "QBASEINSERT" };
 707     push @output, $record->{ "SNUMINSERT" };
 708     push @output, $record->{ "SBASEINSERT" };
 709     push @output, $record->{ "STRAND" };
 710     push @output, $record->{ "Q_ID" };
 711     push @output, $record->{ "Q_LEN" };
 712     push @output, $record->{ "Q_BEG" };
 713     push @output, $record->{ "Q_END" } + 1;
 714     push @output, $record->{ "S_ID" };
 715     push @output, $record->{ "S_LEN" };
 716     push @output, $record->{ "S_BEG" };
 717     push @output, $record->{ "S_END" } + 1;
 718     push @output, $record->{ "BLOCKCOUNT" };
 719     push @output, $record->{ "BLOCKSIZES" };
 720     push @output, $record->{ "Q_BEGS" };
 721     push @output, $record->{ "S_BEGS" };
 722
 723     print $fh join( "\t", @output ), "\n";
 724 }
 725
 726
 727 sub psl_upload_to_ucsc
 728 {
 729     # Martin A. Hansen, September 2007.
 730
 731     # Upload a PSL file to the UCSC database.
 732
 733     my ( $file,      # file to upload,
 734          $options,   # argument hashref
 735          $append,    # flag indicating table should be appended
 736        ) = @_;
 737
 738     # Returns nothing.
 739
 740     my ( $args );
 741
 742     if ( $append ) {
 743         $args = join " ", $options->{ "database" }, "-table=$options->{ 'table' }", "-clientLoad", "-append", $file;
 744     } else {
 745         $args = join " ", $options->{ "database" }, "-table=$options->{ 'table' }", "-clientLoad", $file;
 746     }
 747
 748     Maasha::Common::run( "hgLoadPsl", "$args > /dev/null 2>&1" );
 749 }
 750
 751
 752 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> TRACK FILE <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 753
 754
 755 sub update_my_tracks
 756 {
 757     # Martin A. Hansen, September 2007.
 758
 759     # Update the /home/user/ucsc/my_tracks.ra file and executes makeCustomTracks.pl
 760
 761     my ( $options,   # hashref
 762          $type,      # track type
 763        ) = @_;
 764
 765     # Returns nothing.
 766
 767     my ( $file, $fh_in, $fh_out, $line, $time );
 768
 769     $file = $ENV{ "HOME" } . "/ucsc/my_tracks.ra";
 770
 771     # ---- create a backup ----
 772
 773     $fh_in  = Maasha::Common::read_open( $file );
 774     $fh_out = Maasha::Common::write_open( "$file~" );
 775
 776     while ( $line = <$fh_in> ) {
 777         print $fh_out $line;
 778     }
 779
 780     close $fh_in;
 781     close $fh_out;
 782
 783     # ---- append track ----
 784
 785     $time = Maasha::Common::time_stamp();
 786
 787     $fh_out = Maasha::Common::append_open( $file );
 788
 789     if ( $type eq "sec_struct" )
 790     {
 791         print $fh_out "\n\n# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n";
 792
 793         print $fh_out "\n# Track added by 'upload_to_ucsc' $time\n\n";
 794
 795         print $fh_out "# Database $options->{ 'database' }\n\n";
 796
 797         print $fh_out "track $options->{ 'table' }\n";
 798         print $fh_out "shortLabel $options->{ 'short_label' }\n";
 799         print $fh_out "longLabel $options->{ 'long_label' }\n";
 800         print $fh_out "group $options->{ 'group' }\n";
 801         print $fh_out "priority $options->{ 'priority' }\n";
 802         print $fh_out "visibility $options->{ 'visibility' }\n";
 803         print $fh_out "color $options->{ 'color' }\n";
 804         print $fh_out "type bed 6 +\n";
 805         print $fh_out "mafTrack multiz17way\n";
 806
 807         print $fh_out "\n# //\n";
 808     }
 809     else
 810     {
 811         print $fh_out "\n\n# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n";
 812
 813         print $fh_out "\n# Track added by 'upload_to_ucsc' $time\n\n";
 814
 815         print $fh_out "# Database $options->{ 'database' }\n\n";
 816
 817         print $fh_out "track $options->{ 'table' }\n";
 818         print $fh_out "shortLabel $options->{ 'short_label' }\n";
 819         print $fh_out "longLabel $options->{ 'long_label' }\n";
 820         print $fh_out "group $options->{ 'group' }\n";
 821         print $fh_out "priority $options->{ 'priority' }\n";
 822         print $fh_out "useScore 1\n" if $options->{ 'use_score' };
 823         print $fh_out "visibility $options->{ 'visibility' }\n";
 824         print $fh_out "maxHeightPixels 50:50:11\n" if $type eq "wig 0";
 825         print $fh_out "color $options->{ 'color' }\n";
 826         print $fh_out "type $type\n";
 827
 828         print $fh_out "\n# //\n";
 829     }
 830
 831     close $fh_out;
 832
 833     Maasha::Common::run( "ucscMakeTracks.pl", "-b > /dev/null 2>&1" );
 834 }
 835
 836
 837 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> fixedStep format <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 838
 839
 840 sub fixedstep_get_entry
 841 {
 842     # Martin A. Hansen, December 2007.
 843
 844     # Given a file handle to a PhastCons file get the
 845     # next entry which is all the lines after a "fixedStep"
 846     # line and until the next "fixedStep" line or EOF.
 847
 848     my ( $fh,   # filehandle
 849        ) = @_;
 850
 851     # Returns a list of lines
 852
 853     my ( $entry, @lines );
 854
 855     local $/ = "\nfixedStep ";
 856
 857     $entry = <$fh>;
 858
 859     chomp $entry;
 860
 861     @lines = split "\n", $entry;
 862
 863     return if @lines == 0;
 864
 865     $lines[ 0 ] =~ s/fixedStep?\s*//;
 866
 867     return wantarray ? @lines : \@lines;
 868 }
 869
 870
 871 sub fixedstep_index_create
 872 {
 873     # Martin A. Hansen, January 2008.
 874
 875     # Indexes a concatenated fixedStep file.
 876     # The index consists of a hash with chromosomes as keys,
 877     # and a list of [ chr_beg, next_chr_beg, chr_end, index_beg, index_len ] as values.
 878
 879     my ( $path,   # path to fixedStep file
 880        ) = @_;
 881
 882     # Returns a hashref
 883
 884     my ( $fh, $pos, $index_beg, $index_len, $entry, $locator, $chr, $step, $beg, $end, $len, %index, $i );
 885
 886     $fh = Maasha::Common::read_open( $path );
 887
 888     $pos = 0;
 889
 890     while ( $entry = Maasha::UCSC::fixedstep_get_entry( $fh ) )
 891     {
 892         $locator = shift @{ $entry };
 893
 894         if ( $locator =~ /chrom=([^ ]+) start=(\d+) step=(\d+)/ )
 895         {
 896             $chr  = $1;
 897             $beg  = $2 - 1;  #  fixedStep files are 1-based
 898             $step = $3;
 899         }
 900         else
 901         {
 902             Maasha::Common::error( qq(Could not parse locator: $locator) );
 903         }
 904
 905         $pos += length( $locator ) + 11;
 906
 907         $index_beg = $pos;
 908
 909 #        map { $pos += length( $_ ) + 1 } @{ $entry };
 910
 911         $pos += 6 * scalar @{ $entry };
 912
 913         $index_len = $pos - $index_beg;
 914
 915         push @{ $index{ $chr } }, [ $beg, undef, $beg + scalar @{ $entry } - 1, $index_beg, $index_len ];
 916     }
 917
 918     close $fh;
 919
 920     foreach $chr ( keys %index )
 921     {
 922         for ( $i = 0; $i < @{ $index{ $chr } } - 1; $i++ ) {
 923             $index{ $chr }->[ $i ]->[ NEXT_CHR_BEG ] = $index{ $chr }->[ $i + 1 ]->[ 0 ];
 924         }
 925
 926         $index{ $chr }->[ -1 ]->[ NEXT_CHR_BEG ] = $index{ $chr }->[ -1 ]->[ CHR_END ] + 1;
 927     }
 928
 929     return wantarray ? %index : \%index;
 930 }
 931
 932
 933 sub fixedstep_index_store
 934 {
 935     # Martin A. Hansen, January 2008.
 936
 937     # Writes a fixedStep index to binary file.
 938
 939     my ( $path,   # full path to file
 940          $index,  # list with index
 941        ) = @_;
 942
 943     # returns nothing
 944
 945     Maasha::Common::file_store( $path, $index );
 946 }
 947
 948
 949 sub fixedstep_index_retrieve
 950 {
 951     # Martin A. Hansen, January 2008.
 952
 953     # Retrieves a fixedStep index from binary file.
 954
 955     my ( $path,   # full path to file
 956        ) = @_;
 957
 958     # returns list
 959
 960     my $index;
 961
 962     $index = Maasha::Common::file_retrieve( $path );
 963
 964     return wantarray ? %{ $index } : $index;
 965 }
 966
 967
 968 sub fixedStep_index_lookup
 969 {
 970     # Martin A. Hansen, January 2008.
 971
 972     # Retrieve fixedStep scores from a indexed
 973     # fixedStep file given a chromosome and
 974     # begin and end positions.
 975
 976     my ( $index,     # data structure
 977          $fh,        # filehandle to datafile
 978          $chr,       # chromosome
 979          $chr_beg,   # chromosome beg
 980          $chr_end,   # chromosome end
 981          $flank,     # include flanking region - OPTIONAL
 982        ) = @_;
 983
 984     # Returns a list
 985
 986     my ( $index_beg, $index_end, $i, $c, $beg, $end, @vals, $scores );
 987
 988     $flank ||= 0;
 989
 990     $chr_beg -= $flank;
 991     $chr_end += $flank;
 992
 993 #    print "chr_beg->$chr_beg   chr_end->$chr_end   flank->$flank\n";
 994
 995     if ( exists $index->{ $chr } )
 996     {
 997         $index_beg = Maasha::Matrix::interval_search( $index->{ $chr }, 0, 1, $chr_beg );
 998
 999         if ( $index_beg < 0 ) {
1000             Maasha::Common::error( qq(Index search failed - begin index position doesn't exists: $chr_beg) );
1001         }
1002
1003         if ( $chr_end < $index->{ $chr }->[ $index_beg ]->[ 1 ] )
1004         {
1005             $index_end = $index_beg;
1006         }
1007         else
1008         {
1009             $index_end = Maasha::Matrix::interval_search( $index->{ $chr }, 0, 1, $chr_end );
1010
1011             if ( $index_end < 0 ) {
1012                 Maasha::Common::error( qq(Index search failed - end index position doesn't exists: $chr_end) );
1013             }
1014         }
1015
1016         map { $scores->[ $_ ] = 0 } 0 .. $chr_end - $chr_beg;
1017
1018         if ( $index_beg == $index_end )
1019         {
1020             $beg = Maasha::Calc::max( $chr_beg, $index->{ $chr }->[ $index_beg ]->[ CHR_BEG ] );
1021             $end = Maasha::Calc::min( $chr_end, $index->{ $chr }->[ $index_end ]->[ CHR_END ] );
1022
1023             if ( $beg <= $index->{ $chr }->[ $index_beg ]->[ CHR_END ] and $end >= $index->{ $chr }->[ $index_beg ]->[ CHR_BEG ] )
1024             {
1025                 @vals = split "\n", Maasha::Common::file_read(
1026                     $fh,
1027                     $index->{ $chr }->[ $index_beg ]->[ INDEX_BEG ] + 6 * ( $beg - $index->{ $chr }->[ $index_beg ]->[ CHR_BEG ] ),
1028                     6 * ( $end - $beg + 1 ),
1029                 );
1030             }
1031
1032             for ( $c = 0; $c < @vals; $c++ ) {
1033                 $scores->[ $c + $beg - $chr_beg ] = $vals[ $c ];
1034             }
1035         }
1036         else
1037         {
1038             $beg = Maasha::Calc::max( $chr_beg, $index->{ $chr }->[ $index_beg ]->[ CHR_BEG ] );
1039
1040 #            print Dumper( $beg, $index->{ $chr }->[ $index_beg ] );
1041 #            print Dumper( "next", $index->{ $chr }->[ $index_beg ]->[ NEXT_CHR_BEG ] );
1042
1043             #      beg         next
1044             #      v           v
1045             #  |||||||||.......
1046
1047             if ( $beg <= $index->{ $chr }->[ $index_beg ]->[ CHR_END ] )
1048             {
1049                 @vals = split "\n", Maasha::Common::file_read(
1050                     $fh,
1051                     $index->{ $chr }->[ $index_beg ]->[ INDEX_BEG ] + 6 * ( $beg - $index->{ $chr }->[ $index_beg ]->[ CHR_BEG ] ),
1052                     6 * ( $index->{ $chr }->[ $index_beg ]->[ CHR_END ] - $beg + 1 ),
1053                 );
1054
1055                 for ( $c = 0; $c < @vals; $c++ ) {
1056                     $scores->[ $c + $beg - $chr_beg ] = $vals[ $c ];
1057                 }
1058             }
1059
1060             $end = Maasha::Calc::min( $chr_end, $index->{ $chr }->[ $index_end ]->[ CHR_END ] );
1061
1062             if ( $end <= $index->{ $chr }->[ $index_end ]->[ CHR_END ] )
1063             {
1064                 @vals = split "\n", Maasha::Common::file_read(
1065                     $fh,
1066                     $index->{ $chr }->[ $index_end ]->[ INDEX_BEG ],
1067                     6 * ( $end - $index->{ $chr }->[ $index_end ]->[ CHR_BEG ] + 1 ),
1068                 );
1069
1070                 for ( $c = 0; $c < @vals; $c++ ) {
1071                     $scores->[ $c + $index->{ $chr }->[ $index_end ]->[ CHR_BEG ] - $chr_beg ] = $vals[ $c ];
1072                 }
1073             }
1074
1075             for ( $i = $index_beg + 1; $i <= $index_end - 1; $i++ )
1076             {
1077                 @vals = split "\n", Maasha::Common::file_read(
1078                     $fh,
1079                     $index->{ $chr }->[ $i ]->[ INDEX_BEG ],
1080                     6 * ( $index->{ $chr }->[ $i ]->[ CHR_END ] - $index->{ $chr }->[ $i ]->[ CHR_BEG ] + 1 ),
1081                 );
1082
1083                 for ( $c = 0; $c < @vals; $c++ ) {
1084                     $scores->[ $c + $index->{ $chr }->[ $i ]->[ CHR_BEG ] - $chr_beg ] = $vals[ $c ];
1085                 }
1086             }
1087         }
1088     }
1089     else
1090     {
1091         Maasha::Common::error( qq(Chromosome "$chr" was not found in index) );
1092     }
1093
1094     return wantarray ? @{ $scores } : $scores;
1095 }
1096
1097
1098 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> PhastCons format <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1099
1100
1101 sub phastcons_index
1102 {
1103     # Martin A. Hansen, July 2008
1104
1105     # Create a fixedStep index for PhastCons data.
1106
1107     my ( $file,   # file to index
1108          $dir,    # dir with file
1109        ) = @_;
1110
1111     # Returns nothing.
1112
1113     my ( $index );
1114
1115     $index = fixedstep_index_create( "$dir/$file" );
1116
1117     fixedstep_index_store( "$dir/$file.index", $index );
1118 }
1119
1120
1121 sub phastcons_parse_entry
1122 {
1123     # Martin A. Hansen, December 2007.
1124
1125     # Given a PhastCons entry converts this to a
1126     # list of super blocks.
1127
1128     my ( $lines,   # list of lines
1129          $args,    # argument hash
1130        ) = @_;
1131
1132     # Returns
1133
1134     my ( $info, $chr, $beg, $step, $i, $c, $j, @blocks, @super_blocks, @entries, $super_block, $block, @lens, @begs );
1135
1136     $info = shift @{ $lines };
1137
1138     if ( $info =~ /^chrom=([^ ]+) start=(\d+) step=(\d+)$/ )
1139     {
1140         $chr  = $1;
1141         $beg  = $2;
1142         $step = $3;
1143
1144         die qq(ERROR: step size $step != 1 -> problem!\n) if $step != 1; # in an ideal world should would be fixed ...
1145     }
1146
1147     $i = 0;
1148
1149     while ( $i < @{ $lines } )
1150     {
1151         if ( $lines->[ $i ] >= $args->{ "threshold" } )
1152         {
1153             $c = $i + 1;
1154
1155             while ( $c < @{ $lines } )
1156             {
1157                 if ( $lines->[ $c ] < $args->{ "threshold" } )
1158                 {
1159                     $j = $c + 1;
1160
1161                     while ( $j < @{ $lines } and $lines->[ $j ] < $args->{ "threshold" } ) {
1162                         $j++;
1163                     }
1164
1165                     if ( $j - $c > $args->{ "gap" } )
1166                     {
1167                         if ( $c - $i >= $args->{ "min" } )
1168                         {
1169                             push @blocks, {
1170                                 CHR     => $chr,
1171                                 CHR_BEG => $beg + $i - 1,
1172                                 CHR_END => $beg + $c - 2,
1173                                 CHR_LEN => $c - $i,
1174                             };
1175                         }
1176
1177                         $i = $j;
1178
1179                         last;
1180                     }
1181
1182                     $c = $j
1183                 }
1184                 else
1185                 {
1186                     $c++;
1187                 }
1188             }
1189
1190             if ( $c - $i >= $args->{ "min" } )
1191             {
1192                 push @blocks, {
1193                     CHR     => $chr,
1194                     CHR_BEG => $beg + $i - 1,
1195                     CHR_END => $beg + $c - 2,
1196                     CHR_LEN => $c - $i,
1197                 };
1198             }
1199
1200             $i = $c;
1201         }
1202         else
1203         {
1204             $i++;
1205         }
1206     }
1207
1208     $i = 0;
1209
1210     while ( $i < @blocks )
1211     {
1212         $c = $i + 1;
1213
1214         while ( $c < @blocks and $blocks[ $c ]->{ "CHR_BEG" } - $blocks[ $c - 1 ]->{ "CHR_END" } <= $args->{ "dist" } )
1215         {
1216             $c++;
1217         }
1218
1219         push @super_blocks, [ @blocks[ $i .. $c - 1 ] ];
1220
1221         $i = $c;
1222     }
1223
1224     foreach $super_block ( @super_blocks )
1225     {
1226         foreach $block ( @{ $super_block } )
1227         {
1228             push @begs, $block->{ "CHR_BEG" } - $super_block->[ 0 ]->{ "CHR_BEG" };
1229             push @lens, $block->{ "CHR_LEN" } - 1;
1230         }
1231
1232         $lens[ -1 ]++;
1233
1234         push @entries, {
1235             CHR        => $super_block->[ 0 ]->{ "CHR" },
1236             CHR_BEG    => $super_block->[ 0 ]->{ "CHR_BEG" },
1237             CHR_END    => $super_block->[ -1 ]->{ "CHR_END" },
1238             Q_ID       => "Q_ID",
1239             SCORE      => 100,
1240             STRAND     => "+",
1241             THICK_BEG  => $super_block->[ 0 ]->{ "CHR_BEG" },
1242             THICK_END  => $super_block->[ -1 ]->{ "CHR_END" } + 1,
1243             ITEMRGB    => "0,200,100",
1244             BLOCKCOUNT => scalar @{ $super_block },
1245             BLOCKSIZES => join( ",", @lens ),
1246             Q_BEGS     => join( ",", @begs ),
1247         };
1248
1249         undef @begs;
1250         undef @lens;
1251     }
1252
1253     return wantarray ? @entries : \@entries;
1254 }
1255
1256
1257 sub phastcons_normalize
1258 {
1259     # Martin A. Hansen, January 2008.
1260
1261     # Normalizes a list of lists with PhastCons scores,
1262     # in such a way that each list contains the same number
1263     # or PhastCons scores.
1264
1265     my ( $AoA,    # AoA with PhastCons scores
1266        ) = @_;
1267
1268     # Returns AoA.
1269
1270     my ( $list, $max, $min, $mean, $diff );
1271
1272     $min = 99999999;
1273     $max = 0;
1274
1275     foreach $list ( @{ $AoA } )
1276     {
1277         $min = scalar @{ $list } if scalar @{ $list } < $min;
1278         $max = scalar @{ $list } if scalar @{ $list } > $max;
1279     }
1280
1281     $mean = int( ( $min + $max ) / 2 );
1282
1283 #    print STDERR "min->$min   max->$max   mean->$mean\n";
1284
1285     foreach $list ( @{ $AoA } )
1286     {
1287         $diff = scalar @{ $list } - $mean;
1288
1289         phastcons_list_inflate( $list, abs( $diff ) ) if $diff < 0;
1290         phastcons_list_deflate( $list, $diff )        if $diff > 0;
1291     }
1292
1293     return wantarray ? @{ $AoA } : $AoA;
1294 }
1295
1296
1297 sub phastcons_list_inflate
1298 {
1299     # Martin A. Hansen, January 2008.
1300
1301     # Inflates a list with a given number of elements
1302     # in such a way that the extra elements are introduced
1303     # evenly over the entire length of the list. The value
1304     # of the extra elements is based on a mean of the
1305     # adjacent elements.
1306
1307     my ( $list,   # list of elements
1308          $diff,   # number of elements to introduce
1309        ) = @_;
1310
1311     # Returns nothing
1312
1313     my ( $len, $space, $i, $pos );
1314
1315     $len = scalar @{ $list };
1316
1317     $space = $len / $diff;
1318
1319     for ( $i = 0; $i < $diff; $i++ )
1320     {
1321         $pos = int( ( $space / 2 ) + $i * $space );
1322
1323         splice @{ $list }, $pos, 0, ( $list->[ $pos - 1 ] + $list->[ $pos + 1 ] ) / 2;
1324         # splice @{ $list }, $pos, 0, "X";
1325     }
1326
1327     die qq(ERROR: bad inflate\n) if scalar @{ $list } != $len + $diff;
1328 }
1329
1330
1331 sub phastcons_list_deflate
1332 {
1333     # Martin A. Hansen, January 2008.
1334
1335     # Deflates a list by removing a given number of elements
1336     # evenly distributed over the entire list.
1337
1338     my ( $list,   # list of elements
1339          $diff,   # number of elements to remove
1340        ) = @_;
1341
1342     # Returns nothing
1343
1344     my ( $len, $space, $i, $pos );
1345
1346     $len = scalar @{ $list };
1347
1348     $space = ( $len - $diff ) / $diff;
1349
1350     for ( $i = 0; $i < $diff; $i++ )
1351     {
1352         $pos = int( ( $space / 2 ) + $i * $space );
1353
1354         splice @{ $list }, $pos, 1;
1355     }
1356
1357     die qq(ERROR: bad deflate\n) if scalar @{ $list } != $len - $diff;
1358 }
1359
1360
1361 sub phastcons_mean
1362 {
1363     # Martin A. Hansen, January 2008.
1364
1365     # Given a normalized PhastCons matrix in an AoA,
1366     # calculate the mean for each column and return as a list.
1367
1368     my ( $AoA,    # AoA with normalized PhastCons scores
1369        ) = @_;
1370
1371     # Returns a list
1372
1373     my ( @list );
1374
1375     $AoA = Maasha::Matrix::matrix_flip( $AoA );
1376
1377     map { push @list, Maasha::Calc::mean( $_ ) } @{ $AoA };
1378
1379     return wantarray ? @list : \@list;
1380 }
1381
1382
1383 sub phastcons_median
1384 {
1385     # Martin A. Hansen, January 2008.
1386
1387     # Given a normalized PhastCons matrix in an AoA,
1388     # calculate the median for each column and return as a list.
1389
1390     my ( $AoA,    # AoA with normalized PhastCons scores
1391        ) = @_;
1392
1393     # Returns a list
1394
1395     my ( @list );
1396
1397     $AoA = Maasha::Matrix::matrix_flip( $AoA );
1398
1399     map { push @list, Maasha::Calc::median( $_ ) } @{ $AoA };
1400
1401     return wantarray ? @list : \@list;
1402 }
1403
1404
1405 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> MULTIPLE ALIGNMENT FILES <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1406
1407
1408 sub maf_extract
1409 {
1410     # Martin A. Hansen, April 2008.
1411
1412     # Executes mafFrag to extract a subalignment from a multiz track
1413     # in the UCSC genome browser database.
1414
1415     my ( $tmp_dir,    # temporary directory
1416          $database,   # genome database
1417          $table,      # table with the multiz track
1418          $chr,        # chromosome
1419          $beg,        # begin position
1420          $end,        # end position
1421          $strand,     # strand
1422        ) = @_;
1423
1424     # Returns a list of record
1425
1426     my ( $tmp_file, $align );
1427
1428     $tmp_file = "$tmp_dir/maf_extract.maf";
1429
1430     Maasha::Common::run( "mafFrag", "$database $table $chr $beg $end $strand $tmp_file" );
1431
1432     $align = maf_parse( $tmp_file );
1433
1434     unlink $tmp_file;
1435
1436     return wantarray ? @{ $align } : $align;
1437 }
1438
1439
1440 sub maf_parse
1441 {
1442     # Martin A. Hansen, April 2008.
1443
1444
1445     my ( $path,   # full path to MAF file
1446        ) = @_;
1447
1448     # Returns a list of record.
1449
1450     my ( $fh, $line, @fields, @align );
1451
1452     $fh = Maasha::Common::read_open( $path );
1453
1454     while ( $line = <$fh> )
1455     {
1456         chomp $line;
1457
1458         if ( $line =~ /^s/ )
1459         {
1460             @fields = split / /, $line;
1461
1462             push @align, {
1463                 SEQ_NAME  => $fields[ 1 ],
1464                 SEQ       => $fields[ -1 ],
1465                 ALIGN     => 1,
1466                 ALIGN_LEN => length $fields[ -1 ],
1467             }
1468         }
1469     }
1470
1471     close $fh;
1472
1473     return wantarray ? @align : \@align;
1474 }
1475
1476
1477 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> WIGGLE FORMAT <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1478
1479
1480 sub fixedstep_put_entry
1481 {
1482     # Martin A. Hansen, April 2008.
1483
1484     # Outputs a block of fixedStep values.
1485     # Used for outputting wiggle data.
1486
1487     my ( $chr,      # chromosome
1488          $beg,      # start position
1489          $block,    # list of scores
1490          $fh,       # filehandle - OPTIONAL
1491          $log10,    # flag indicating that log10 scores should be used
1492        ) = @_;
1493
1494     # Returns nothing.
1495
1496     $beg += 1;   # fixedStep format is 1 based.
1497
1498     $fh ||= \*STDOUT;
1499
1500     print $fh "fixedStep chrom=$chr start=$beg step=1\n";
1501
1502     if ( $log10 ) {
1503         map { printf( $fh "%f\n", Maasha::Calc::log10( $_ + 1 ) ) } @{ $block };
1504     } else {
1505         map { printf( $fh "%d\n", ( $_ + 1 ) ) } @{ $block };
1506     }
1507 }
1508
1509
1510 sub wiggle_upload_to_ucsc
1511 {
1512     # Martin A. Hansen, May 2008.
1513
1514     # Upload a wiggle file to the UCSC database.
1515
1516     my ( $tmp_dir,    # temporary directory
1517          $wib_dir,    # wib directory
1518          $wig_file,   # file to upload,
1519          $options,    # argument hashref
1520        ) = @_;
1521
1522     # Returns nothing.
1523
1524     my ( $args );
1525
1526 #    $args = join " ", "-tmpDir=$tmp_dir", "-pathPrefix=$wib_dir", $options->{ "database" }, $options->{ 'table' }, $wig_file;
1527
1528 #    Maasha::Common::run( "hgLoadWiggle", "$args > /dev/null 2>&1" );
1529
1530     `cd $tmp_dir && hgLoadWiggle -tmpDir=$tmp_dir -pathPrefix=$wib_dir $options->{ 'database' } $options->{ 'table' } $wig_file > /dev/null 2>&1`;
1531 }
1532
1533
1534 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> MySQL CONF <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1535
1536
1537 sub ucsc_get_user
1538 {
1539     # Martin A. Hansen, May 2008
1540
1541     # Fetches the MySQL database user name from the
1542     # .hg.conf file in the users home directory.
1543
1544     # Returns a string.
1545
1546     my ( $fh, $line, $user );
1547
1548     $fh = Maasha::Common::read_open( "$ENV{ 'HOME' }/.hg.conf" );
1549
1550     while ( $line = <$fh> )
1551     {
1552         chomp $line;
1553
1554         if ( $line =~ /^db\.user=(.+)/ )
1555         {
1556             $user = $1;
1557
1558             last;
1559         }
1560     }
1561
1562     close $fh;
1563
1564     return $user;
1565 }
1566
1567
1568 sub ucsc_get_password
1569 {
1570     # Martin A. Hansen, May 2008
1571
1572     # Fetches the MySQL database password from the
1573     # .hg.conf file in the users home directory.
1574
1575     # Returns a string.
1576
1577     my ( $fh, $line, $password );
1578
1579     $fh = Maasha::Common::read_open( "$ENV{ 'HOME' }/.hg.conf" );
1580
1581     while ( $line = <$fh> )
1582     {
1583         chomp $line;
1584
1585         if ( $line =~ /^db\.password=(.+)/ )
1586         {
1587             $password = $1;
1588
1589             last;
1590         }
1591     }
1592
1593     close $fh;
1594
1595     return $password;
1596 }
1597
1598
1599 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1600
1601
1602 __END__
1603