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 =~ s/(\n|\r)$//g;    # 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 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> PhastCons format <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 838
 839
 840 sub phastcons_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 phastcons_parse_entry
 872 {
 873     # Martin A. Hansen, December 2007.
 874
 875     # Given a PhastCons entry converts this to a
 876     # list of super blocks.
 877
 878     my ( $lines,   # list of lines
 879          $args,    # argument hash
 880        ) = @_;
 881
 882     # Returns
 883
 884     my ( $info, $chr, $beg, $step, $i, $c, $j, @blocks, @super_blocks, @entries, $super_block, $block, @lens, @begs );
 885
 886     $info = shift @{ $lines };
 887
 888     if ( $info =~ /^chrom=([^ ]+) start=(\d+) step=(\d+)$/ )
 889     {
 890         $chr  = $1;
 891         $beg  = $2;
 892         $step = $3;
 893
 894         die qq(ERROR: step size $step != 1 -> problem!\n) if $step != 1; # in an ideal world should would be fixed ...
 895     }
 896
 897     $i = 0;
 898
 899     while ( $i < @{ $lines } )
 900     {
 901         if ( $lines->[ $i ] >= $args->{ "threshold" } )
 902         {
 903             $c = $i + 1;
 904
 905             while ( $c < @{ $lines } )
 906             {
 907                 if ( $lines->[ $c ] < $args->{ "threshold" } )
 908                 {
 909                     $j = $c + 1;
 910
 911                     while ( $j < @{ $lines } and $lines->[ $j ] < $args->{ "threshold" } ) {
 912                         $j++;
 913                     }
 914
 915                     if ( $j - $c > $args->{ "gap" } )
 916                     {
 917                         if ( $c - $i >= $args->{ "min" } )
 918                         {
 919                             push @blocks, {
 920                                 CHR     => $chr,
 921                                 CHR_BEG => $beg + $i - 1,
 922                                 CHR_END => $beg + $c - 2,
 923                                 CHR_LEN => $c - $i,
 924                             };
 925                         }
 926
 927                         $i = $j;
 928
 929                         last;
 930                     }
 931
 932                     $c = $j
 933                 }
 934                 else
 935                 {
 936                     $c++;
 937                 }
 938             }
 939
 940             if ( $c - $i >= $args->{ "min" } )
 941             {
 942                 push @blocks, {
 943                     CHR     => $chr,
 944                     CHR_BEG => $beg + $i - 1,
 945                     CHR_END => $beg + $c - 2,
 946                     CHR_LEN => $c - $i,
 947                 };
 948             }
 949
 950             $i = $c;
 951         }
 952         else
 953         {
 954             $i++;
 955         }
 956     }
 957
 958     $i = 0;
 959
 960     while ( $i < @blocks )
 961     {
 962         $c = $i + 1;
 963
 964         while ( $c < @blocks and $blocks[ $c ]->{ "CHR_BEG" } - $blocks[ $c - 1 ]->{ "CHR_END" } <= $args->{ "dist" } )
 965         {
 966             $c++;
 967         }
 968
 969         push @super_blocks, [ @blocks[ $i .. $c - 1 ] ];
 970
 971         $i = $c;
 972     }
 973
 974     foreach $super_block ( @super_blocks )
 975     {
 976         foreach $block ( @{ $super_block } )
 977         {
 978             push @begs, $block->{ "CHR_BEG" } - $super_block->[ 0 ]->{ "CHR_BEG" };
 979             push @lens, $block->{ "CHR_LEN" } - 1;
 980         }
 981
 982         $lens[ -1 ]++;
 983
 984         push @entries, {
 985             CHR        => $super_block->[ 0 ]->{ "CHR" },
 986             CHR_BEG    => $super_block->[ 0 ]->{ "CHR_BEG" },
 987             CHR_END    => $super_block->[ -1 ]->{ "CHR_END" },
 988             Q_ID       => "Q_ID",
 989             SCORE      => 100,
 990             STRAND     => "+",
 991             THICK_BEG  => $super_block->[ 0 ]->{ "CHR_BEG" },
 992             THICK_END  => $super_block->[ -1 ]->{ "CHR_END" } + 1,
 993             ITEMRGB    => "0,200,100",
 994             BLOCKCOUNT => scalar @{ $super_block },
 995             BLOCKSIZES => join( ",", @lens ),
 996             Q_BEGS     => join( ",", @begs ),
 997         };
 998
 999         undef @begs;
1000         undef @lens;
1001     }
1002
1003     return wantarray ? @entries : \@entries;
1004 }
1005
1006
1007 sub phastcons_index_create
1008 {
1009     # Martin A. Hansen, January 2008.
1010
1011     # Indexes a concatenated PhastCons file.
1012     # The index consists of a hash with chromosomes as keys,
1013     # and a list of [ chr_beg, next_chr_beg, chr_end, index_beg, index_len ] as values.
1014
1015     my ( $path,   # path to PhastCons file
1016        ) = @_;
1017
1018     # Returns a hashref
1019
1020     my ( $fh, $pos, $index_beg, $index_len, $entry, $locator, $chr, $step, $beg, $end, $len, %index, $i );
1021
1022     $fh = &Maasha::Common::read_open( $path );
1023
1024     $pos = 0;
1025
1026     while ( $entry = &Maasha::UCSC::phastcons_get_entry( $fh ) )
1027     {
1028         $locator = shift @{ $entry };
1029
1030         if ( $locator =~ /chrom=([^ ]+) start=(\d+) step=(\d+)/ )
1031         {
1032             $chr  = $1;
1033             $beg  = $2 - 1;  #  phastcons files are 1-based
1034             $step = $3;
1035         }
1036         else
1037         {
1038             &Maasha::Common::error( qq(Could not parse PhastCons locator: $locator) );
1039         }
1040
1041         $pos += length( $locator ) + 11;
1042
1043         $index_beg = $pos;
1044
1045 #        map { $pos += length( $_ ) + 1 } @{ $entry };
1046
1047         $pos += 6 * scalar @{ $entry };
1048
1049         $index_len = $pos - $index_beg;
1050
1051         push @{ $index{ $chr } }, [ $beg, undef, $beg + scalar @{ $entry } - 1, $index_beg, $index_len ];
1052     }
1053
1054     close $fh;
1055
1056     foreach $chr ( keys %index )
1057     {
1058         for ( $i = 0; $i < @{ $index{ $chr } } - 1; $i++ ) {
1059             $index{ $chr }->[ $i ]->[ NEXT_CHR_BEG ] = $index{ $chr }->[ $i + 1 ]->[ 0 ];
1060         }
1061
1062         $index{ $chr }->[ -1 ]->[ NEXT_CHR_BEG ] = $index{ $chr }->[ -1 ]->[ CHR_END ] + 1;
1063     }
1064
1065     return wantarray ? %index : \%index;
1066 }
1067
1068
1069 sub phastcons_index_store
1070 {
1071     # Martin A. Hansen, January 2008.
1072
1073     # Writes a PhastCons index to binary file.
1074
1075     my ( $path,   # full path to file
1076          $index,  # list with index
1077        ) = @_;
1078
1079     # returns nothing
1080
1081     &Maasha::Common::file_store( $path, $index );
1082 }
1083
1084
1085 sub phastcons_index_retrieve
1086 {
1087     # Martin A. Hansen, January 2008.
1088
1089     # Retrieves a PhastCons index from binary file.
1090
1091     my ( $path,   # full path to file
1092        ) = @_;
1093
1094     # returns list
1095
1096     my $index;
1097
1098     $index = &Maasha::Common::file_retrieve( $path );
1099
1100     return wantarray ? %{ $index } : $index;
1101 }
1102
1103
1104 sub phastcons_index_lookup
1105 {
1106     # Martin A. Hansen, January 2008.
1107
1108     # Retrieve PhastCons scores from a indexed
1109     # Phastcons file given a chromosome and
1110     # begin and end positions.
1111
1112     my ( $index,     # data structure
1113          $fh,        # filehandle to datafile
1114          $chr,       # chromosome
1115          $chr_beg,   # chromosome beg
1116          $chr_end,   # chromosome end
1117          $flank,     # include flanking region - OPTIONAL
1118        ) = @_;
1119
1120     # Returns a list
1121
1122     my ( $index_beg, $index_end, $i, $c, $beg, $end, @vals, $scores );
1123
1124     $flank ||= 0;
1125
1126     $chr_beg -= $flank;
1127     $chr_end += $flank;
1128
1129 #    print "chr_beg->$chr_beg   chr_end->$chr_end   flank->$flank\n";
1130
1131     if ( exists $index->{ $chr } )
1132     {
1133         $index_beg = &Maasha::Matrix::interval_search( $index->{ $chr }, 0, 1, $chr_beg );
1134
1135         if ( $index_beg < 0 ) {
1136             &Maasha::Common::error( qq(Index search failed - begin index position doesn't exists: $chr_beg) );
1137         }
1138
1139         if ( $chr_end < $index->{ $chr }->[ $index_beg ]->[ 1 ] )
1140         {
1141             $index_end = $index_beg;
1142         }
1143         else
1144         {
1145             $index_end = &Maasha::Matrix::interval_search( $index->{ $chr }, 0, 1, $chr_end );
1146
1147             if ( $index_end < 0 ) {
1148                 &Maasha::Common::error( qq(Index search failed - end index position doesn't exists: $chr_end) );
1149             }
1150         }
1151
1152         map { $scores->[ $_ ] = 0 } 0 .. $chr_end - $chr_beg;
1153
1154         if ( $index_beg == $index_end )
1155         {
1156             $beg = &Maasha::Calc::max( $chr_beg, $index->{ $chr }->[ $index_beg ]->[ CHR_BEG ] );
1157             $end = &Maasha::Calc::min( $chr_end, $index->{ $chr }->[ $index_end ]->[ CHR_END ] );
1158
1159             if ( $beg <= $index->{ $chr }->[ $index_beg ]->[ CHR_END ] and $end >= $index->{ $chr }->[ $index_beg ]->[ CHR_BEG ] )
1160             {
1161                 @vals = split "\n", &Maasha::Common::file_read(
1162                     $fh,
1163                     $index->{ $chr }->[ $index_beg ]->[ INDEX_BEG ] + 6 * ( $beg - $index->{ $chr }->[ $index_beg ]->[ CHR_BEG ] ),
1164                     6 * ( $end - $beg + 1 ),
1165                 );
1166             }
1167
1168             for ( $c = 0; $c < @vals; $c++ ) {
1169                 $scores->[ $c + $beg - $chr_beg ] = $vals[ $c ];
1170             }
1171         }
1172         else
1173         {
1174             $beg = &Maasha::Calc::max( $chr_beg, $index->{ $chr }->[ $index_beg ]->[ CHR_BEG ] );
1175
1176 #            print Dumper( $beg, $index->{ $chr }->[ $index_beg ] );
1177 #            print Dumper( "next", $index->{ $chr }->[ $index_beg ]->[ NEXT_CHR_BEG ] );
1178
1179             #      beg         next
1180             #      v           v
1181             #  |||||||||.......
1182
1183             if ( $beg <= $index->{ $chr }->[ $index_beg ]->[ CHR_END ] )
1184             {
1185                 @vals = split "\n", &Maasha::Common::file_read(
1186                     $fh,
1187                     $index->{ $chr }->[ $index_beg ]->[ INDEX_BEG ] + 6 * ( $beg - $index->{ $chr }->[ $index_beg ]->[ CHR_BEG ] ),
1188                     6 * ( $index->{ $chr }->[ $index_beg ]->[ CHR_END ] - $beg + 1 ),
1189                 );
1190
1191                 for ( $c = 0; $c < @vals; $c++ ) {
1192                     $scores->[ $c + $beg - $chr_beg ] = $vals[ $c ];
1193                 }
1194             }
1195
1196             $end = &Maasha::Calc::min( $chr_end, $index->{ $chr }->[ $index_end ]->[ CHR_END ] );
1197
1198             if ( $end <= $index->{ $chr }->[ $index_end ]->[ CHR_END ] )
1199             {
1200                 @vals = split "\n", &Maasha::Common::file_read(
1201                     $fh,
1202                     $index->{ $chr }->[ $index_end ]->[ INDEX_BEG ],
1203                     6 * ( $end - $index->{ $chr }->[ $index_end ]->[ CHR_BEG ] + 1 ),
1204                 );
1205
1206                 for ( $c = 0; $c < @vals; $c++ ) {
1207                     $scores->[ $c + $index->{ $chr }->[ $index_end ]->[ CHR_BEG ] - $chr_beg ] = $vals[ $c ];
1208                 }
1209             }
1210
1211             for ( $i = $index_beg + 1; $i <= $index_end - 1; $i++ )
1212             {
1213                 @vals = split "\n", &Maasha::Common::file_read(
1214                     $fh,
1215                     $index->{ $chr }->[ $i ]->[ INDEX_BEG ],
1216                     6 * ( $index->{ $chr }->[ $i ]->[ CHR_END ] - $index->{ $chr }->[ $i ]->[ CHR_BEG ] + 1 ),
1217                 );
1218
1219                 for ( $c = 0; $c < @vals; $c++ ) {
1220                     $scores->[ $c + $index->{ $chr }->[ $i ]->[ CHR_BEG ] - $chr_beg ] = $vals[ $c ];
1221                 }
1222             }
1223         }
1224     }
1225     else
1226     {
1227         &Maasha::Common::error( qq(Chromosome "$chr" was not found in index) );
1228     }
1229
1230     return wantarray ? @{ $scores } : $scores;
1231 }
1232
1233
1234 sub phastcons_normalize
1235 {
1236     # Martin A. Hansen, January 2008.
1237
1238     # Normalizes a list of lists with PhastCons scores,
1239     # in such a way that each list contains the same number
1240     # or PhastCons scores.
1241
1242     my ( $AoA,    # AoA with PhastCons scores
1243        ) = @_;
1244
1245     # Returns AoA.
1246
1247     my ( $list, $max, $min, $mean, $diff );
1248
1249     $min = 99999999;
1250     $max = 0;
1251
1252     foreach $list ( @{ $AoA } )
1253     {
1254         $min = scalar @{ $list } if scalar @{ $list } < $min;
1255         $max = scalar @{ $list } if scalar @{ $list } > $max;
1256     }
1257
1258     $mean = int( ( $min + $max ) / 2 );
1259
1260 #    print STDERR "min->$min   max->$max   mean->$mean\n";
1261
1262     foreach $list ( @{ $AoA } )
1263     {
1264         $diff = scalar @{ $list } - $mean;
1265
1266         &phastcons_list_inflate( $list, abs( $diff ) ) if $diff < 0;
1267         &phastcons_list_deflate( $list, $diff )        if $diff > 0;
1268     }
1269
1270     return wantarray ? @{ $AoA } : $AoA;
1271 }
1272
1273
1274 sub phastcons_list_inflate
1275 {
1276     # Martin A. Hansen, January 2008.
1277
1278     # Inflates a list with a given number of elements
1279     # in such a way that the extra elements are introduced
1280     # evenly over the entire length of the list. The value
1281     # of the extra elements is based on a mean of the
1282     # adjacent elements.
1283
1284     my ( $list,   # list of elements
1285          $diff,   # number of elements to introduce
1286        ) = @_;
1287
1288     # Returns nothing
1289
1290     my ( $len, $space, $i, $pos );
1291
1292     $len = scalar @{ $list };
1293
1294     $space = $len / $diff;
1295
1296     for ( $i = 0; $i < $diff; $i++ )
1297     {
1298         $pos = int( ( $space / 2 ) + $i * $space );
1299
1300         splice @{ $list }, $pos, 0, ( $list->[ $pos - 1 ] + $list->[ $pos + 1 ] ) / 2;
1301         # splice @{ $list }, $pos, 0, "X";
1302     }
1303
1304     die qq(ERROR: bad inflate\n) if scalar @{ $list } != $len + $diff;
1305 }
1306
1307
1308 sub phastcons_list_deflate
1309 {
1310     # Martin A. Hansen, January 2008.
1311
1312     # Deflates a list by removing a given number of elements
1313     # evenly distributed over the entire list.
1314
1315     my ( $list,   # list of elements
1316          $diff,   # number of elements to remove
1317        ) = @_;
1318
1319     # Returns nothing
1320
1321     my ( $len, $space, $i, $pos );
1322
1323     $len = scalar @{ $list };
1324
1325     $space = ( $len - $diff ) / $diff;
1326
1327     for ( $i = 0; $i < $diff; $i++ )
1328     {
1329         $pos = int( ( $space / 2 ) + $i * $space );
1330
1331         splice @{ $list }, $pos, 1;
1332     }
1333
1334     die qq(ERROR: bad deflate\n) if scalar @{ $list } != $len - $diff;
1335 }
1336
1337
1338 sub phastcons_mean
1339 {
1340     # Martin A. Hansen, January 2008.
1341
1342     # Given a normalized PhastCons matrix in an AoA,
1343     # calculate the mean for each column and return as a list.
1344
1345     my ( $AoA,    # AoA with normalized PhastCons scores
1346        ) = @_;
1347
1348     # Returns a list
1349
1350     my ( @list );
1351
1352     $AoA = &Maasha::Matrix::matrix_flip( $AoA );
1353
1354     map { push @list, &Maasha::Calc::mean( $_ ) } @{ $AoA };
1355
1356     return wantarray ? @list : \@list;
1357 }
1358
1359
1360 sub phastcons_median
1361 {
1362     # Martin A. Hansen, January 2008.
1363
1364     # Given a normalized PhastCons matrix in an AoA,
1365     # calculate the median for each column and return as a list.
1366
1367     my ( $AoA,    # AoA with normalized PhastCons scores
1368        ) = @_;
1369
1370     # Returns a list
1371
1372     my ( @list );
1373
1374     $AoA = &Maasha::Matrix::matrix_flip( $AoA );
1375
1376     map { push @list, &Maasha::Calc::median( $_ ) } @{ $AoA };
1377
1378     return wantarray ? @list : \@list;
1379 }
1380
1381
1382 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> MULTIPLE ALIGNMENT FILES <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1383
1384
1385 sub maf_extract
1386 {
1387     # Martin A. Hansen, April 2008.
1388
1389     # Executes mafFrag to extract a subalignment from a multiz track
1390     # in the UCSC genome browser database.
1391
1392     my ( $tmp_dir,    # temporary directory
1393          $database,   # genome database
1394          $table,      # table with the multiz track
1395          $chr,        # chromosome
1396          $beg,        # begin position
1397          $end,        # end position
1398          $strand,     # strand
1399        ) = @_;
1400
1401     # Returns a list of record
1402
1403     my ( $tmp_file, $align );
1404
1405     $tmp_file = "$tmp_dir/maf_extract.maf";
1406
1407     &Maasha::Common::run( "mafFrag", "$database $table $chr $beg $end $strand $tmp_file" );
1408
1409     $align = &maf_parse( $tmp_file );
1410
1411     unlink $tmp_file;
1412
1413     return wantarray ? @{ $align } : $align;
1414 }
1415
1416
1417 sub maf_parse
1418 {
1419     # Martin A. Hansen, April 2008.
1420
1421
1422     my ( $path,   # full path to MAF file
1423        ) = @_;
1424
1425     # Returns a list of record.
1426
1427     my ( $fh, $line, @fields, @align );
1428
1429     $fh = &Maasha::Common::read_open( $path );
1430
1431     while ( $line = <$fh> )
1432     {
1433         chomp $line;
1434
1435         if ( $line =~ /^s/ )
1436         {
1437             @fields = split / /, $line;
1438
1439             push @align, {
1440                 SEQ_NAME  => $fields[ 1 ],
1441                 SEQ       => $fields[ -1 ],
1442                 ALIGN     => 1,
1443                 ALIGN_LEN => length $fields[ -1 ],
1444             }
1445         }
1446     }
1447
1448     close $fh;
1449
1450     return wantarray ? @align : \@align;
1451 }
1452
1453
1454 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> WIGGLE FORMAT <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1455
1456
1457 sub fixedstep_put_entry
1458 {
1459     # Martin A. Hansen, April 2008.
1460
1461     # Outputs a block of fixedStep values.
1462     # Used for outputting wiggle data.
1463
1464     my ( $chr,      # chromosome
1465          $beg,      # start position
1466          $block,    # list of scores
1467          $fh,       # filehandle - OPTIONAL
1468        ) = @_;
1469
1470     # Returns nothing.
1471
1472     $beg += 1;   # fixedStep format is 1 based.
1473
1474     if ( $fh )
1475     {
1476         print $fh "fixedStep chrom=$chr start=$beg step=1\n";
1477
1478         map { printf( $fh "%d\n", ( $_ + 1 ) ) } @{ $block };
1479     }
1480     else
1481     {
1482         print "fixedStep chrom=$chr start=$beg step=1\n";
1483
1484         map { printf( "%d\n", ( $_ + 1 ) ) } @{ $block };
1485     }
1486 }
1487
1488
1489 sub wiggle_upload_to_ucsc
1490 {
1491     # Martin A. Hansen, May 2008.
1492
1493     # Upload a wiggle file to the UCSC database.
1494
1495     my ( $tmp_dir,    # temporary directory
1496          $wib_dir,    # wib directory
1497          $wig_file,   # file to upload,
1498          $options,    # argument hashref
1499        ) = @_;
1500
1501     # Returns nothing.
1502
1503     my ( $args );
1504
1505 #    $args = join " ", "-tmpDir=$tmp_dir", "-pathPrefix=$wib_dir", $options->{ "database" }, $options->{ 'table' }, $wig_file;
1506
1507 #    &Maasha::Common::run( "hgLoadWiggle", "$args > /dev/null 2>&1" );
1508
1509     `cd $tmp_dir && hgLoadWiggle -tmpDir=$tmp_dir -pathPrefix=$wib_dir $options->{ 'database' } $options->{ 'table' } $wig_file > /dev/null 2>&1`;
1510 }
1511
1512
1513 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> MySQL CONF <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1514
1515
1516 sub ucsc_get_user
1517 {
1518     # Martin A. Hansen, May 2008
1519
1520     # Fetches the MySQL database user name from the
1521     # .hg.conf file in the users home directory.
1522
1523     # Returns a string.
1524
1525     my ( $fh, $line, $user );
1526
1527     $fh = &Maasha::Common::read_open( "$ENV{ 'HOME' }/.hg.conf" );
1528
1529     while ( $line = <$fh> )
1530     {
1531         chomp $line;
1532
1533         if ( $line =~ /^db\.user=(.+)/ )
1534         {
1535             $user = $1;
1536
1537             last;
1538         }
1539     }
1540
1541     close $fh;
1542
1543     return $user;
1544 }
1545
1546
1547 sub ucsc_get_password
1548 {
1549     # Martin A. Hansen, May 2008
1550
1551     # Fetches the MySQL database password from the
1552     # .hg.conf file in the users home directory.
1553
1554     # Returns a string.
1555
1556     my ( $fh, $line, $password );
1557
1558     $fh = &Maasha::Common::read_open( "$ENV{ 'HOME' }/.hg.conf" );
1559
1560     while ( $line = <$fh> )
1561     {
1562         chomp $line;
1563
1564         if ( $line =~ /^db\.password=(.+)/ )
1565         {
1566             $password = $1;
1567
1568             last;
1569         }
1570     }
1571
1572     close $fh;
1573
1574     return $password;
1575 }
1576
1577
1578 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1579
1580
1581 __END__
1582