forester/archive/perl/extractTrembl.pl

   1 #!/usr/bin/perl -W
   2
   3 # extractTrembl.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: 04/24/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 AC, DE, and OS from a "trembl.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 and their full species names replaced with their SWISS-PROT
  22 #          code (recommended).
  23 #
  24 # Usage.   extractTrembl.pl <infile> <outfile> [species list]
  25
  26 # Remark.  Need to re-run this if species in species tree or species list
  27 #          are added/changed or if a new version of Pfam is used!!
  28
  29 # Some "heuristic" is required for Synechococcus, Synechocystis, Anabaena:
  30 # see below.
  31
  32 use strict;
  33
  34
  35 my $VERSION       = "1.001";
  36 my $infile        = "";
  37 my $outfile       = "";
  38 my $speciesfile   = "";
  39 my $return_line   = "";
  40 my $read          = 0;
  41 my $ac            = "";
  42 my $de            = "";
  43 my $os            = "";
  44 my %Species_names = (); # full name -> SWISS-PROT name.
  45 my $i             = 0;
  46
  47 if ( @ARGV != 2 && @ARGV != 3 ) {
  48     &errorInCommandLine();
  49 }
  50
  51 $infile  = $ARGV[ 0 ];
  52 $outfile = $ARGV[ 1 ];
  53
  54 if ( @ARGV == 3 ) {
  55     $speciesfile = $ARGV[ 2 ];
  56     unless ( ( -s $speciesfile ) && ( -f $speciesfile ) && ( -T $speciesfile ) ) {
  57         die "\n$0: <<$speciesfile>> does not exist, is empty, or is not a plain textfile.\n\n";
  58     }
  59     &readSpeciesNamesFile( $speciesfile );
  60 }
  61
  62 if ( -e $outfile ) {
  63     die "\n$0: <<$outfile>> already exists.\n\n";
  64 }
  65 unless ( ( -s $infile ) && ( -f $infile ) && ( -T $infile ) ) {
  66     die "\n$0: <<$infile>> does not exist, is empty, or is not a plain textfile.\n\n";
  67 }
  68
  69 open( IN, "$infile" ) || die "\n$0: Cannot open file <<$infile>>: $!\n";
  70 open( OUT, ">$outfile" ) || die "\n$0: Cannot create file <<$outfile>>: $!\n";
  71
  72 print OUT "# extractTrembl.pl version: $VERSION\n";
  73 print OUT "# trembl.dat file: $infile\n";
  74 print OUT "# output file    : $outfile\n";
  75 print OUT "# species file   : $speciesfile\n";
  76 print OUT "# date           : ".`date`."\n\n";
  77
  78 $read = 0;
  79
  80 while ( $return_line = <IN> ) {
  81     if ( $return_line =~ /^AC\s+(\S+);/ ) {
  82         $ac = $1;
  83         $read = 1;
  84     }
  85     elsif ( $return_line =~ /^DE\s+(.+)/ && $read == 1 ) {
  86         if ( $de ne "" ) {
  87             $de .= " ".$1;
  88         }
  89         else {
  90             $de = $1;
  91         }
  92     }
  93
  94     elsif ( $return_line =~ /^OS\s+(.+)\.\s*$/ && $read == 1 ) {
  95         $os = $1;
  96         if ( $speciesfile ne ""
  97         && ( $os =~ /Synechococcus/
  98         ||   $os =~ /Synechocystis/
  99         ||   $os =~ /Anabaena/ ) ) {
 100             if ( $os =~ /PCC\s*(\d+)/ ) {
 101                 $os = "PCC ".$1;
 102             }
 103             else {
 104                 $read = 0;
 105                 $ac = "";
 106                 $de = "";
 107                 $os = "";
 108                 next;
 109             }
 110         }
 111         else {
 112             $os =~ s/\(.+\)//g;
 113         }
 114         $os =~ s/^\s+//;
 115         $os =~ s/\s+$//;
 116         if ( $speciesfile ne "" ) {
 117             unless ( exists( $Species_names{ $os } ) ) {
 118                 $read = 0;
 119                 $ac = "";
 120                 $de = "";
 121                 $os = "";
 122                 next;
 123             }
 124             $os = $Species_names{ $os };
 125         }
 126     }
 127     elsif ( $return_line =~ /^\/\// && $read == 1 ) {
 128         $read = 0;
 129         print OUT "$ac;$de;$os\n";
 130         $ac = "";
 131         $de = "";
 132         $os = "";
 133         $i++;
 134     }
 135 }
 136
 137 close( IN );
 138
 139 print OUT "\n # $i entries.\n";
 140
 141 close( OUT );
 142
 143 exit( 0 );
 144
 145
 146
 147 # Reads in species file.
 148 # Format: SWISS-PROT=full name (e.g. "BACSU=Bacillus subtilis")
 149 # Lines beginning with "#" are ignored.
 150 # One argument: species file-name
 151 # Last modified: 04/24/01
 152 sub readSpeciesNamesFile {
 153     my $infile = $_[ 0 ];
 154     my $return_line = "";
 155     my $sp          = "";
 156     my $full        = "";
 157
 158     unless ( ( -s $infile ) && ( -f $infile ) && ( -T $infile ) ) {
 159         die "\n$0: readSpeciesNamesFile: <<$infile>> does not exist, is empty, or is not a plain textfile.\n";
 160     }
 161
 162     open( IN_RSNF, "$infile" ) || die "\n$0: Cannot open file <<$infile>>: $!\n";
 163     while ( $return_line = <IN_RSNF> ) {
 164         if ( $return_line !~ /^\s*#/ && $return_line =~ /(\S+)=(.+)/ ) {
 165             $sp   = $1;
 166             $full = $2;
 167             $full =~ s/^\s+//;
 168             $full =~ s/\s+$//;
 169             $Species_names{ $full } = $sp;
 170         }
 171     }
 172     close( IN_RSNF );
 173
 174     return;
 175 }
 176
 177
 178
 179 sub errorInCommandLine {
 180     print "\n";
 181     print " extractTrembl.pl  $VERSION\n";
 182     print " ----------------\n";
 183     print "\n";
 184     print " Christian Zmasek (zmasek\@genetics.wustl.edu)\n";
 185     print "\n";
 186     print " Purpose. Extracts AC, DE, and OS from a \"trembl.dat\" file.\n";
 187     print "          The resulting output is used by \"rio.pl\".\n";
 188     print "          If a species list (format: SWISS-PROT-code=full name) is supplied:\n";
 189     print "          only sequences from species found in this list are written to\n";
 190     print "          outfile and their full species names replaced with their SWISS-PROT\n";
 191     print "          code (recommended).\n";
 192     print "\n";
 193     print " Usage.   extractTrembl.pl <infile> <outfile> [species list]\n";
 194     print "\n";
 195     print " Remark.  Need to re-run this if species in species tree or species list\n";
 196     print "          are added/changed or if a new version of Pfam is used!!\n";
 197     print "\n\n";
 198     exit( -1 );
 199 }