code_perl/Maasha/UCSC/PSL.pm

   1 package Maasha::UCSC::PSL;
   2
   3 # Copyright (C) 2007-2008 Martin A. Hansen.
   4
   5 # This program is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU General Public License
   7 # as published by the Free Software Foundation; either version 2
   8 # of the License, or (at your option) any later version.
   9
  10 # This program is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 # GNU General Public License for more details.
  14
  15 # You should have received a copy of the GNU General Public License
  16 # along with this program; if not, write to the Free Software
  17 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18
  19 # http://www.gnu.org/copyleft/gpl.html
  20
  21
  22 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  23
  24
  25 # Routines for interaction with PSL entries and files.
  26
  27 # Read more about the PSL format here: http://genome.ucsc.edu/goldenPath/help/hgTracksHelp.html#PSL
  28
  29
  30 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  31
  32
  33 use warnings;
  34 use strict;
  35
  36 use Data::Dumper;
  37 use Maasha::Common;
  38 use Maasha::Filesys;
  39
  40 use vars qw( @ISA @EXPORT_OK );
  41
  42 require Exporter;
  43
  44 @ISA = qw( Exporter );
  45
  46 @EXPORT_OK = qw(
  47 );
  48
  49 use constant {
  50     matches     => 0,   # PSL field names
  51     misMatches  => 1,
  52     repMatches  => 2,
  53     nCount      => 3,
  54     qNumInsert  => 4,
  55     qBaseInsert => 5,
  56     tNumInsert  => 6,
  57     tBaseInsert => 7,
  58     strand      => 8,
  59     qName       => 9,
  60     qSize       => 10,
  61     qStart      => 11,
  62     qEnd        => 12,
  63     tName       => 13,
  64     tSize       => 14,
  65     tStart      => 15,
  66     tEnd        => 16,
  67     blockCount  => 17,
  68     blockSizes  => 18,
  69     qStarts     => 19,
  70     tStarts     => 20,
  71     MATCHES     => 0,    # Biopieces field names
  72     MISMATCHES  => 1,
  73     REPMATCHES  => 2,
  74     NCOUNT      => 3,
  75     QNUMINSERT  => 4,
  76     QBASEINSERT => 5,
  77     SNUMINSERT  => 6,
  78     SBASEINSERT => 7,
  79     STRAND      => 8,
  80     Q_ID        => 9,
  81     Q_LEN       => 10,
  82     Q_BEG       => 11,
  83     Q_END       => 12,
  84     S_ID        => 13,
  85     S_LEN       => 14,
  86     S_BEG       => 15,
  87     S_END       => 16,
  88     BLOCK_COUNT => 17,
  89     BLOCK_LENS  => 18,
  90     Q_BEGS      => 19,
  91     S_BEGS      => 20,
  92 };
  93
  94
  95 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> PSL format <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  96
  97
  98 sub psl_entry_get
  99 {
 100     # Martin A. Hansen, August 2008.
 101
 102     # Reads PSL next entry from a PSL file and returns an array ref.
 103
 104     my ( $fh,   # file handle of PSL filefull path to PSL file
 105        ) = @_;
 106
 107     # Returns hashref.
 108
 109     my ( $line, @fields );
 110
 111     while ( $line = <$fh> )
 112     {
 113         next if $line !~ /^\d/;
 114
 115         chomp $line;
 116
 117         @fields = split "\t", $line;
 118
 119         return wantarray ? @fields : \@fields;   # 21 fields
 120     }
 121
 122     return undef;
 123 }
 124
 125
 126 sub psl2biopiece
 127 {
 128     # Martin A. Hansen, November 2008
 129
 130     # Converts PSL record keys to Biopiece record keys.
 131
 132     my ( $psl,   # hashref
 133        ) = @_;
 134
 135     # returns hashref
 136
 137     my %record;
 138
 139     %record = (
 140         REC_TYPE     => "PSL",
 141         BLOCK_LENS   => $psl->[ blockSizes ],
 142         SNUMINSERT   => $psl->[ tNumInsert ],
 143         Q_END        => $psl->[ qEnd ] - 1,
 144         SBASEINSERT  => $psl->[ tBaseInsert ],
 145         S_END        => $psl->[ tEnd ] - 1,
 146         QBASEINSERT  => $psl->[ qBaseInsert ],
 147         REPMATCHES   => $psl->[ repMatches ],
 148         QNUMINSERT   => $psl->[ qNumInsert ],
 149         MISMATCHES   => $psl->[ misMatches ],
 150         BLOCK_COUNT  => $psl->[ blockCount ],
 151         Q_LEN        => $psl->[ qSize ],
 152         S_ID         => $psl->[ tName ],
 153         STRAND       => $psl->[ strand ],
 154         Q_ID         => $psl->[ qName ],
 155         MATCHES      => $psl->[ matches ],
 156         S_LEN        => $psl->[ tSize ],
 157         NCOUNT       => $psl->[ nCount ],
 158         Q_BEGS       => $psl->[ qStarts ],
 159         S_BEGS       => $psl->[ tStarts ],
 160         S_BEG        => $psl->[ tStart ],
 161         Q_BEG        => $psl->[ qStart ],
 162     );
 163
 164     $record{ "SCORE" } = $record{ "MATCHES" } + int( $record{ "REPMATCHES" } / 2 ) - $record{ "MISMATCHES" } - $record{ "QNUMINSERT" } - $record{ "SNUMINSERT" };
 165     $record{ "SPAN" }  = $record{ "Q_END" } - $record{ "Q_BEG" };
 166
 167     return wantarray ? %record : \%record;
 168 }
 169
 170
 171 sub psl_calc_score
 172 {
 173     # Martin A. Hansen, November 2008.
 174
 175     # Calculates the score for a PSL entry according to:
 176     # http://genome.ucsc.edu/FAQ/FAQblat#blat4
 177
 178     my ( $psl_entry,   # array ref
 179        ) = @_;
 180
 181     # Returns a float
 182
 183     my ( $score );
 184
 185     $score = $psl_entry->[ MATCHES ] +
 186              int( $psl_entry->[ REPMATCHES ] / 2 ) -
 187              $psl_entry->[ MISMATCHES ] -
 188              $psl_entry->[ QNUMINSERT ] -
 189              $psl_entry->[ SNUMINSERT ];
 190
 191     return $score;
 192 }
 193
 194
 195 sub psl_put_header
 196 {
 197     # Martin A. Hansen, September 2007.
 198
 199     # Write a PSL header to file.
 200
 201     my ( $fh,  # file handle  - OPTIONAL
 202        ) = @_;
 203
 204     # Returns nothing.
 205
 206     $fh = \*STDOUT if not $fh;
 207
 208     print $fh qq(psLayout version 3
 209 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
 210 ---------------------------------------------------------------------------------------------------------------------------------------------------------------
 211 );
 212 }
 213
 214
 215 sub psl_put_entry
 216 {
 217     # Martin A. Hansen, September 2007.
 218
 219     # Write a PSL entry to file.
 220
 221     my ( $record,       # hashref
 222          $fh,           # file handle  -  OPTIONAL
 223        ) = @_;
 224
 225     # Returns nothing.
 226
 227     $fh = \*STDOUT if not $fh;
 228
 229     my @output;
 230
 231     push @output, $record->{ "MATCHES" };
 232     push @output, $record->{ "MISMATCHES" };
 233     push @output, $record->{ "REPMATCHES" };
 234     push @output, $record->{ "NCOUNT" };
 235     push @output, $record->{ "QNUMINSERT" };
 236     push @output, $record->{ "QBASEINSERT" };
 237     push @output, $record->{ "SNUMINSERT" };
 238     push @output, $record->{ "SBASEINSERT" };
 239     push @output, $record->{ "STRAND" };
 240     push @output, $record->{ "Q_ID" };
 241     push @output, $record->{ "Q_LEN" };
 242     push @output, $record->{ "Q_BEG" };
 243     push @output, $record->{ "Q_END" } + 1;
 244     push @output, $record->{ "S_ID" };
 245     push @output, $record->{ "S_LEN" };
 246     push @output, $record->{ "S_BEG" };
 247     push @output, $record->{ "S_END" } + 1;
 248     push @output, $record->{ "BLOCK_COUNT" };
 249     push @output, $record->{ "BLOCK_LENS" };
 250     push @output, $record->{ "Q_BEGS" };
 251     push @output, $record->{ "S_BEGS" };
 252
 253     print $fh join( "\t", @output ), "\n";
 254 }
 255
 256
 257 sub psl_upload_to_ucsc
 258 {
 259     # Martin A. Hansen, September 2007.
 260
 261     # Upload a PSL file to the UCSC database.
 262
 263     my ( $file,      # file to upload,
 264          $options,   # argument hashref
 265          $append,    # flag indicating table should be appended
 266        ) = @_;
 267
 268     # Returns nothing.
 269
 270     my ( $args );
 271
 272     if ( $append ) {
 273         $args = join " ", $options->{ "database" }, "-table=$options->{ 'table' }", "-clientLoad", "-append", $file;
 274     } else {
 275         $args = join " ", $options->{ "database" }, "-table=$options->{ 'table' }", "-clientLoad", $file;
 276     }
 277
 278     Maasha::Common::run( "hgLoadPsl", "$args > /dev/null 2>&1" );
 279 }
 280
 281
 282 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 283
 284
 285 1;
 286
 287
 288 __END__