code_perl/Maasha/NCBI.pm

   1 package Maasha::NCBI;
   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 NCBI Entrez
  26
  27
  28 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  29
  30
  31 use strict;
  32 use Data::Dumper;
  33 use LWP::Simple;
  34 use Maasha::Common;
  35 use Maasha::Filesys;
  36 use Maasha::Seq;
  37
  38 use vars qw( @ISA @EXPORT );
  39
  40 @ISA = qw( Exporter );
  41
  42
  43 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  44
  45
  46 sub get_entry
  47 {
  48     # Martin A. Hansen, March 2007.
  49
  50     # connects to the ncbi website and retrieves a genbank record,
  51     # which is returned.
  52
  53     my ( $db,    # database <nucleotide|protein>
  54          $id,    # genbank id
  55          $type,  # retrieval type <gb|gp>
  56        ) = @_;
  57
  58     # returns string
  59
  60     my ( $content, @lines, $i, $seq );
  61
  62     $content = get "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=$db&id=$id&rettype=$type";
  63
  64     return $content;
  65 }
  66
  67
  68 sub get_seq
  69 {
  70     # Martin A. Hansen, March 2007.
  71
  72     # connects to the ncbi website and retrieves a genbank record,
  73     # from which the sequence is parsed and returned.
  74
  75     my ( $db,    # database <nucleotide|protein>
  76          $id,    # genbank id
  77          $type,  # retrieval type <gb|gp>
  78        ) = @_;
  79
  80     # returns string
  81
  82     my ( $content, @lines, $i, $seq );
  83
  84     $content = get "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=$db&id=$id&rettype=$type";
  85
  86     @lines = split "\n", $content;
  87
  88     $i = 0;
  89
  90     while ( $lines[ $i ] !~ /^ORIGIN/ ) {
  91         $i++
  92     }
  93
  94     $i++;
  95
  96     while ( $lines[ $i ] !~ /^\/\// )
  97     {
  98         $lines[ $i ] =~ s/^\s*\d+//;
  99
 100         $seq .= $lines[ $i ];
 101
 102         $i++;
 103     }
 104
 105     $seq =~ tr/ //d;
 106
 107     return $seq;
 108 }
 109
 110
 111 sub soft_parse
 112 {
 113     # Martin A. Hansen, February 2008.
 114
 115     # !!! NOT USED ANYMORE !!! #
 116
 117     # Reads in and parses a file in SOFT format.
 118
 119     my ( $path,   # full path to SOFT file
 120        ) = @_;
 121
 122     # Returns a list.
 123
 124     my ( $fh, @lines, $i, $c, $num, %key_hash, @fields, %id_hash, $id, $seq, $count, $record, @records, $platform_id, $sample_id, $sample_title );
 125
 126     $fh = Maasha::Common::read_open( $path );
 127
 128     @lines = <$fh>;
 129
 130     close $fh;
 131
 132     chomp @lines;
 133
 134     $i = 0;
 135
 136     $num = 1;
 137
 138     while ( $i < @lines )
 139     {
 140         if ( $lines[ $i ] =~ /^\^PLATFORM = (.+)/ )
 141         {
 142             $platform_id = $1;
 143         }
 144         elsif ( $lines[ $i ] =~ /^!platform_table_begin$/ )
 145         {
 146             @fields = split "\t", $lines[ $i + 1 ];
 147
 148             for ( $c = 0; $c < @fields; $c++ ) {
 149                 $key_hash{ $fields[ $c ] } = $c;
 150             }
 151
 152             $c = $i + 2;
 153
 154             while ( $lines[ $c ] !~ /^!platform_table_end$/ )
 155             {
 156                 @fields = split "\t", $lines[ $c ];
 157
 158                 $id_hash{ $fields[ $key_hash{ "ID" } ] } = $fields[ $key_hash{ "SEQUENCE" } ];
 159
 160                 $c++;
 161             }
 162
 163             $i = $c;
 164         }
 165         elsif ( $lines[ $i ] =~ /^\^SAMPLE = (.+)/ )
 166         {
 167             $sample_id = $1;
 168         }
 169         elsif ( $lines[ $i ] =~ /^!Sample_title = (.+)/ )
 170         {
 171             $sample_title = $1;
 172         }
 173         elsif ( $lines[ $i ] =~ /^!sample_table_begin/ )
 174         {
 175             undef %key_hash;
 176
 177             @fields = split "\t", $lines[ $i + 1 ];
 178
 179             for ( $c = 0; $c < @fields; $c++ ) {
 180                 $key_hash{ $fields[ $c ] } = $c;
 181             }
 182
 183             $c = $i + 2;
 184
 185             while ( $lines[ $c ] !~ /^!sample_table_end$/ )
 186             {
 187                 undef $record;
 188
 189                 @fields = split "\t", $lines[ $c ];
 190
 191                 $id    = $fields[ $key_hash{ "ID_REF" } ];
 192                 $seq   = $id_hash{ $id };
 193                 $count = $fields[ $key_hash{ "VALUE" } ];
 194
 195                 $seq =~ tr/./N/;
 196
 197                 $record->{ "SAMPLE_TITLE" } = $sample_title;
 198                 $record->{ "SEQ" }          = $seq;
 199                 $record->{ "SEQ_NAME" }     = join( "_", $platform_id, $sample_id, $num, $count );
 200
 201                 push @records, $record;
 202
 203                 $c++;
 204                 $num++;
 205             }
 206
 207             $i = $c;
 208
 209             $num = 1;
 210         }
 211
 212         $i++;
 213     }
 214
 215     return wantarray ? @records : \@records;
 216 }
 217
 218
 219 sub soft_index_file
 220 {
 221     # Martin A. Hansen, June 2008.
 222
 223     # Create a index with linenumbers of the different tables
 224     # in a soft file. The index is returned.
 225
 226     my ( $file,   # file to index
 227        ) = @_;
 228
 229     # Returns
 230
 231     my ( $fh, $line, $i, $c, @index, $first );
 232
 233     $fh = Maasha::Common::read_open( $file );
 234
 235     $first = 1;
 236
 237     $i      = 0;
 238     $c      = 0;
 239
 240     while ( $line = <$fh> )
 241     {
 242         chomp $line;
 243
 244         if ( $line =~ /^\^/ )
 245         {
 246             push @index, {
 247                 SECTION  => $line,
 248                 LINE_BEG => $i,
 249             };
 250
 251             if ( not $first )
 252             {
 253                 $index[ $c - 1 ]{ "LINE_END" } = $i - 1;
 254             }
 255             else
 256             {
 257                 $first = 0;
 258             }
 259
 260             $c++;
 261         }
 262
 263         $i++;
 264     }
 265
 266     $index[ $c - 1 ]{ "LINE_END" } = $i - 1;
 267
 268     close $fh;
 269
 270     return wantarray ? @index : \@index;
 271 }
 272
 273
 274 sub soft_get_platform
 275 {
 276     # Martin A. Hansen, June 2008.
 277
 278     # Given a filehandle to a SOFT file parses all the platform tables
 279     # which is returned.
 280
 281     my ( $fh,    # filehandle
 282          $beg,   # line number where platform tables begin
 283          $end,   # line number where platform tables end
 284        ) = @_;
 285
 286     # Returns hashref
 287
 288     my ( $line, @lines, $i, $c, @fields, %key_hash, %id_hash );
 289
 290     $i = 0;
 291
 292     while ( $line = <$fh> )
 293     {
 294         chomp $line;
 295
 296         push @lines, $line if $i >= $beg;
 297
 298         last if $i == $end;
 299
 300         $i++;
 301     }
 302
 303     $i = 0;
 304
 305     while ( $i < @lines )
 306     {
 307         if ( $lines[ $i ] =~ /^!platform_table_begin$/ )
 308         {
 309             @fields = split "\t", $lines[ $i + 1 ];
 310
 311             for ( $c = 0; $c < @fields; $c++ ) {
 312                 $key_hash{ $fields[ $c ] } = $c;
 313             }
 314
 315             $c = $i + 2;
 316
 317             while ( $lines[ $c ] !~ /^!platform_table_end$/ )
 318             {
 319                 @fields = split "\t", $lines[ $c ];
 320
 321                 $id_hash{ $fields[ $key_hash{ "ID" } ] } = $fields[ $key_hash{ "SEQUENCE" } ] || $fields[ $key_hash{ "ID" } ];
 322
 323                 $c++;
 324             }
 325
 326             $i = $c;
 327         }
 328
 329         $i++;
 330     }
 331
 332     return wantarray ? %id_hash : \%id_hash;
 333 }
 334
 335
 336 sub soft_get_sample
 337 {
 338     # Martin A. Hansen, June 2008.
 339
 340     # Given a filehandle to a SOFT file and a platform table
 341     # parses sample records which are returned.
 342
 343     my ( $fh,           # filehandle
 344          $plat_table,   # hashref with platform table
 345          $beg,          # line number where sample table begin
 346          $end,          # line number where sample table end
 347          $skip,         # flag indicating that this sample table should not be parsed.
 348        ) = @_;
 349
 350     # Returns list of hashref
 351
 352     my ( $line, @lines, $i, $c, $platform_id, @fields, %key_hash, $num, $sample_id, $sample_title, $id, $seq, $count, @records, $record );
 353
 354     $i = 0;
 355
 356     while ( $line = <$fh> )
 357     {
 358         if ( not $skip )
 359         {
 360             chomp $line;
 361
 362             push @lines, $line if $i >= $beg;
 363         }
 364
 365         last if $i == $end;
 366
 367         $i++;
 368     }
 369
 370     if ( $skip ) {
 371         return wantarray ? () : [];
 372     }
 373
 374     $i = 0;
 375
 376     $num = 1;
 377
 378     while ( $i < @lines )
 379     {
 380         if ( $lines[ $i ] =~ /^\^SAMPLE = (.+)/ )
 381         {
 382             $sample_id = $1;
 383         }
 384         elsif ( $lines[ $i ] =~ /!Sample_platform_id = (.+)/ )
 385         {
 386             $platform_id = $1;
 387         }
 388         elsif ( $lines[ $i ] =~ /^!Sample_title = (.+)/ )
 389         {
 390             $sample_title = $1;
 391         }
 392         elsif ( $lines[ $i ] =~ /^!sample_table_begin/ )
 393         {
 394             undef %key_hash;
 395
 396             @fields = split "\t", $lines[ $i + 1 ];
 397
 398             for ( $c = 0; $c < @fields; $c++ ) {
 399                 $key_hash{ $fields[ $c ] } = $c;
 400             }
 401
 402             $c = $i + 2;
 403
 404             while ( $lines[ $c ] !~ /^!sample_table_end$/ )
 405             {
 406                 undef $record;
 407
 408                 @fields = split "\t", $lines[ $c ];
 409
 410                 $id    = $fields[ $key_hash{ "ID_REF" } ];
 411                 $seq   = $plat_table->{ $id };
 412                 $count = $fields[ $key_hash{ "VALUE" } ];
 413
 414                 $seq =~ tr/./N/;
 415
 416                 $record->{ "SAMPLE_TITLE" } = $sample_title;
 417                 $record->{ "SEQ" }          = $seq;
 418                 $record->{ "SEQ_NAME" }     = join( "_", $platform_id, $sample_id, $num, $count );
 419
 420                 push @records, $record;
 421
 422                 $c++;
 423                 $num++;
 424             }
 425
 426             $i = $c;
 427
 428             $num = 1;
 429         }
 430
 431         $i++;
 432     }
 433
 434     return wantarray ? @records : \@records;
 435 }
 436
 437
 438 sub blast_index
 439 {
 440     # Martin A. Hansen, Juli 2008.
 441
 442     # Create a BLAST index of a given FASTA file
 443     # in a specified directory using a given name.
 444
 445     my ( $file,         # path to FASTA file
 446          $src_dir,      # source directory with FASTA file
 447          $dst_dir,      # destination directory for BLAST index
 448          $index_type,   # protein or nucleotide
 449          $index_name,   # name of index
 450        ) = @_;
 451
 452     # Returns nothing.
 453
 454     my ( $type );
 455
 456     if ( $index_type =~ /protein/i ) {
 457         $type = "T";
 458     } else {
 459         $type = "F";
 460     }
 461
 462     Maasha::Filesys::dir_create_if_not_exists( $dst_dir );
 463     Maasha::Common::run( "formatdb", "-p $type -i $src_dir/$file -t $index_name -l /dev/null" );
 464     Maasha::Common::run( "mv", "$src_dir/*hr $dst_dir" );
 465     Maasha::Common::run( "mv", "$src_dir/*in $dst_dir" );
 466     Maasha::Common::run( "mv", "$src_dir/*sq $dst_dir" );
 467 }
 468
 469
 470 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 471
 472
 473 __END__