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_index_file
 113 {
 114     # Martin A. Hansen, June 2008.
 115
 116     # Create a index with line numbers of the different tables
 117     # in a soft file. The index is returned.
 118
 119     my ( $file,   # file to index
 120        ) = @_;
 121
 122     # Returns a list
 123
 124     my ( $fh, $line, $i, $c, @index, $first );
 125
 126     $fh = Maasha::Filesys::file_read_open( $file );
 127
 128     $first = 1;
 129
 130     $i     = 0;
 131     $c     = 0;
 132
 133     while ( $line = <$fh> )
 134     {
 135         chomp $line;
 136
 137         if ( $line =~ /^\^/ )
 138         {
 139             push @index, {
 140                 SECTION  => $line,
 141                 LINE_BEG => $i,
 142             };
 143
 144             if ( not $first )
 145             {
 146                 $index[ $c - 1 ]{ "LINE_END" } = $i - 1;
 147             }
 148             else
 149             {
 150                 $first = 0;
 151             }
 152
 153             $c++;
 154         }
 155
 156         $i++;
 157     }
 158
 159     $index[ $c - 1 ]{ "LINE_END" } = $i - 1;
 160
 161     close $fh;
 162
 163     return wantarray ? @index : \@index;
 164 }
 165
 166
 167 sub soft_get_platform
 168 {
 169     # Martin A. Hansen, June 2008.
 170
 171     # Given a filehandle to a SOFT file parses all the platform tables
 172     # which is returned.
 173
 174     my ( $fh,    # filehandle
 175          $beg,   # line number where platform tables begin
 176          $end,   # line number where platform tables end
 177        ) = @_;
 178
 179     # Returns hashref
 180
 181     my ( $line, @lines, $i, $c, @fields, %key_hash, %id_hash );
 182
 183     $i = 0;
 184
 185     while ( $line = <$fh> )
 186     {
 187         chomp $line;
 188
 189         push @lines, $line if $i >= $beg;
 190
 191         last if $i == $end;
 192
 193         $i++;
 194     }
 195
 196     $i = 0;
 197
 198     while ( $i < @lines )
 199     {
 200         if ( $lines[ $i ] =~ /^!platform_table_begin$/ )
 201         {
 202             @fields = split "\t", $lines[ $i + 1 ];
 203
 204             for ( $c = 0; $c < @fields; $c++ ) {
 205                 $key_hash{ $fields[ $c ] } = $c;
 206             }
 207
 208             $c = $i + 2;
 209
 210             while ( $lines[ $c ] !~ /^!platform_table_end$/ )
 211             {
 212                 @fields = split "\t", $lines[ $c ];
 213
 214                 if ( defined $key_hash{ "SEQUENCE" } ) {
 215                     $id_hash{ $fields[ $key_hash{ "ID" } ] } = $fields[ $key_hash{ "SEQUENCE" } ];
 216                 } else {
 217                     $id_hash{ $fields[ $key_hash{ "ID" } ] } = $fields[ $key_hash{ "ID" } ];
 218                 }
 219
 220                 $c++;
 221             }
 222
 223             $i = $c;
 224         }
 225
 226         $i++;
 227     }
 228
 229     return wantarray ? %id_hash : \%id_hash;
 230 }
 231
 232
 233 sub soft_get_sample
 234 {
 235     # Martin A. Hansen, June 2008.
 236
 237     # Given a filehandle to a SOFT file and a platform table
 238     # parses sample records which are returned.
 239
 240     my ( $fh,           # filehandle
 241          $plat_table,   # hashref with platform table
 242          $beg,          # line number where sample table begin
 243          $end,          # line number where sample table end
 244          $skip,         # flag indicating that this sample table should not be parsed.
 245        ) = @_;
 246
 247     # Returns list of hashref
 248
 249     my ( $line, @lines, $i, $c, $platform_id, @fields, %key_hash, $num, $sample_id, $sample_title, $id, $seq, $count, @records, $record );
 250
 251     $i = 0;
 252
 253     while ( $line = <$fh> )
 254     {
 255         if ( not $skip )
 256         {
 257             chomp $line;
 258
 259             push @lines, $line if $i >= $beg;
 260         }
 261
 262         last if $i == $end;
 263
 264         $i++;
 265     }
 266
 267     if ( $skip ) {
 268         return wantarray ? () : [];
 269     }
 270
 271     $i = 0;
 272
 273     $num = 1;
 274
 275     while ( $i < @lines )
 276     {
 277         if ( $lines[ $i ] =~ /^\^SAMPLE = (.+)/ )
 278         {
 279             $sample_id = $1;
 280         }
 281         elsif ( $lines[ $i ] =~ /!Sample_platform_id = (.+)/ )
 282         {
 283             $platform_id = $1;
 284         }
 285         elsif ( $lines[ $i ] =~ /^!Sample_title = (.+)/ )
 286         {
 287             $sample_title = $1;
 288         }
 289         elsif ( $lines[ $i ] =~ /^!sample_table_begin/ )
 290         {
 291             undef %key_hash;
 292
 293             @fields = split "\t", $lines[ $i + 1 ];
 294
 295             for ( $c = 0; $c < @fields; $c++ )
 296             {
 297                 $fields[ $c ] = "ID_REF" if $fields[ $c ] eq "SEQUENCE";
 298                 $fields[ $c ] = "VALUE"  if $fields[ $c ] eq "COUNT";
 299
 300                 $key_hash{ $fields[ $c ] } = $c;
 301             }
 302
 303             $c = $i + 2;
 304
 305             while ( $lines[ $c ] !~ /^!sample_table_end$/ )
 306             {
 307                 undef $record;
 308
 309                 @fields = split "\t", $lines[ $c ];
 310                 $id    = $fields[ $key_hash{ "ID_REF" } ];
 311                 $seq   = $plat_table->{ $id } || $id;
 312                 $count = $fields[ $key_hash{ "VALUE" } ];
 313
 314                 $seq =~ tr/./N/;
 315
 316                 $record->{ "SAMPLE_TITLE" } = $sample_title;
 317                 $record->{ "SEQ" }          = $seq;
 318                 $record->{ "SEQ_NAME" }     = join( "_", $platform_id, $sample_id, $num, $count );
 319
 320                 push @records, $record;
 321
 322                 $c++;
 323                 $num++;
 324             }
 325
 326             $i = $c;
 327
 328             $num = 1;
 329         }
 330
 331         $i++;
 332     }
 333
 334     return wantarray ? @records : \@records;
 335 }
 336
 337
 338 sub blast_index
 339 {
 340     # Martin A. Hansen, Juli 2008.
 341
 342     # Create a BLAST index of a given FASTA file
 343     # in a specified directory using a given name.
 344
 345     my ( $file,         # path to FASTA file
 346          $src_dir,      # source directory with FASTA file
 347          $dst_dir,      # destination directory for BLAST index
 348          $index_type,   # protein or nucleotide
 349          $index_name,   # name of index
 350        ) = @_;
 351
 352     # Returns nothing.
 353
 354     my ( $type );
 355
 356     if ( $index_type =~ /protein/i ) {
 357         $type = "T";
 358     } else {
 359         $type = "F";
 360     }
 361
 362     Maasha::Filesys::dir_create_if_not_exists( $dst_dir );
 363     Maasha::Common::run( "formatdb", "-p $type -i $src_dir/$file -t $index_name -l /dev/null" );
 364     Maasha::Common::run( "mv", "$src_dir/*hr $dst_dir" );
 365     Maasha::Common::run( "mv", "$src_dir/*in $dst_dir" );
 366     Maasha::Common::run( "mv", "$src_dir/*sq $dst_dir" );
 367 }
 368
 369
 370 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 371
 372
 373 __END__