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