forester/archive/perl/extractSWISS-PROT.pl

   1 #!/usr/bin/perl -W
   2
   3 # extractSWISS-PROT.pl
   4 # --------------------
   5 #
   6 # Copyright (C) 2001 Washington University School of Medicine
   7 # and Howard Hughes Medical Institute
   8 # All rights reserved
   9 #
  10 # Created: 09/25/01
  11 # Author: Christian M. Zmasek
  12 # zmasek@genetics.wustl.edu
  13 # http://www.genetics.wustl.edu/eddy/people/zmasek/
  14 #
  15 # Last modified 05/23/02
  16
  17 # Purpose. Extracts ID, DE, and species from a "sprot.dat" file.
  18 #          The output is used by "rio.pl".
  19 #          If a species list (format: SWISS-PROT-code=full name) is supplied:
  20 #          only sequences from species found in this list are written to
  21 #          outfile (recommended).
  22 #
  23 # Usage.   extractSWISS-PROT.pl <infile> <outfile> [species list]
  24
  25 # Remark.  Need to re-run this if species in species tree or species list
  26 #          are added/changed or if a new version of Pfam is used!!
  27
  28
  29 use strict;
  30
  31
  32 my $VERSION       = "1.001";
  33 my $infile        = "";
  34 my $outfile       = "";
  35 my $speciesfile   = "";
  36 my $return_line   = "";
  37 my $read          = 0;
  38 my $ac            = "";
  39 my $de            = "";
  40 my $os            = "";
  41 my %Species_names = (); # SWISS-PROT name -> "".
  42 my $i             = 0;
  43
  44 if ( @ARGV != 2 && @ARGV != 3 ) {
  45     &errorInCommandLine();
  46 }
  47
  48 $infile  = $ARGV[ 0 ];
  49 $outfile = $ARGV[ 1 ];
  50
  51 if ( @ARGV == 3 ) {
  52     $speciesfile = $ARGV[ 2 ];
  53     unless ( ( -s $speciesfile ) && ( -f $speciesfile ) && ( -T $speciesfile ) ) {
  54         die "\n$0: <<$speciesfile>> does not exist, is empty, or is not a plain textfile.\n\n";
  55     }
  56     &readSpeciesNamesFile( $speciesfile );
  57 }
  58
  59 if ( -e $outfile ) {
  60     die "\n$0: <<$outfile>> already exists.\n\n";
  61 }
  62 unless ( ( -s $infile ) && ( -f $infile ) && ( -T $infile ) ) {
  63     die "\n$0: <<$infile>> does not exist, is empty, or is not a plain textfile.\n\n";
  64 }
  65
  66 open( IN, "$infile" ) || die "\n$0: Cannot open file <<$infile>>: $!\n";
  67 open( OUT, ">$outfile" ) || die "\n$0: Cannot create file <<$outfile>>: $!\n";
  68
  69 print OUT "# extractTrembl.pl version: $VERSION\n";
  70 print OUT "# trembl.dat file: $infile\n";
  71 print OUT "# output file    : $outfile\n";
  72 print OUT "# species file   : $speciesfile\n";
  73 print OUT "# date           : ".`date`."\n\n";
  74
  75 $read = 0;
  76
  77 while ( $return_line = <IN> ) {
  78     if ( $return_line =~ /^ID\s+(\S+)/ ) {
  79         $ac = $1;
  80         $read = 1;
  81         if ( $ac =~ /[A-Z0-9]+_([A-Z0-9]+)/ ) {
  82             $os = $1;
  83         }
  84         else {
  85             die "\n$0: Unexpected format: $ac.\n\n";
  86         }
  87         if ( $speciesfile ne "" ) {
  88             unless ( exists( $Species_names{ $os } ) ) {
  89                 $read = 0;
  90                 $ac = "";
  91                 $de = "";
  92                 $os = "";
  93                 next;
  94             }
  95         }
  96     }
  97     elsif ( $return_line =~ /^DE\s+(.+)/ && $read == 1 ) {
  98         if ( $de ne "" ) {
  99             $de .= " ".$1;
 100         }
 101         else {
 102             $de = $1;
 103         }
 104     }
 105     elsif ( $return_line =~ /^\/\// && $read == 1 ) {
 106         $read = 0;
 107         print OUT "$ac;$de;$os\n";
 108         $ac = "";
 109         $de = "";
 110         $os = "";
 111         $i++;
 112     }
 113 }
 114
 115 close( IN );
 116
 117 print OUT "\n # $i entries.\n";
 118
 119 close( OUT );
 120
 121 exit( 0 );
 122
 123
 124
 125 # Reads in species file.
 126 # Format: SWISS-PROT=full name (e.g. "BACSU=Bacillus subtilis")
 127 # Lines beginning with "#" are ignored.
 128 # One argument: species file-name
 129 # Last modified: 04/24/01
 130 sub readSpeciesNamesFile {
 131     my $infile = $_[ 0 ];
 132     my $return_line = "";
 133     my $sp          = "";
 134     my $full        = "";
 135
 136     unless ( ( -s $infile ) && ( -f $infile ) && ( -T $infile ) ) {
 137         die "\n$0: readSpeciesNamesFile: <<$infile>> does not exist, is empty, or is not a plain textfile.\n";
 138     }
 139
 140     open( IN_RSNF, "$infile" ) || die "\n$0: Cannot open file <<$infile>>: $!\n";
 141     while ( $return_line = <IN_RSNF> ) {
 142         if ( $return_line !~ /^\s*#/ && $return_line =~ /(\S+)=(.+)/ ) {
 143             $sp   = $1;
 144             $full = $2;
 145             $full =~ s/^\s+//;
 146             $full =~ s/\s+$//;
 147             $Species_names{ $sp } = "";
 148         }
 149     }
 150     close( IN_RSNF );
 151
 152     return;
 153 }
 154
 155
 156
 157 sub errorInCommandLine {
 158     print "\n";
 159     print " extractSWISS-PROT.pl  $VERSION\n";
 160     print " --------------------\n";
 161     print "\n";
 162     print " Christian Zmasek (zmasek\@genetics.wustl.edu)\n";
 163     print "\n";
 164     print " Purpose. Extracts ID, DE, and species from a \"sprot.dat\" file.\n";
 165     print "          The resulting output is used by \"rio.pl\".\n";
 166     print "          If a species list (format: SWISS-PROT-code=full name) is supplied:\n";
 167     print "          only sequences from species found in this list are written to\n";
 168     print "          outfile (recommended).\n";
 169     print "\n";
 170     print " Usage.   extractSWISS-PROT.pl <infile> <outfile> [species list]\n";
 171     print "\n";
 172     print " Remark.  Need to re-run this if species in species tree or species list\n";
 173     print "          are added/changed or if a new version of Pfam is used!!\n";
 174     print "\n\n";
 175     exit( -1 );
 176 }