[Bio] / FigKernelScripts / check_sims_basic.pl Repository:
ViewVC logotype

Diff of /FigKernelScripts/check_sims_basic.pl

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1, Sun May 30 06:22:40 2004 UTC revision 1.6, Fri Jul 8 18:02:10 2005 UTC
# Line 1  Line 1 
1  $usage = "usage: check_sims_basic NR < sims > checked.sims 2> errors";  # -*- perl -*-
2    
3  ($nr = shift @ARGV)  $SIG{HUP} = 'ignore';
4    
5    use FIG;
6    use File::Path;
7    use File::Basename;
8    
9    $usage = "usage: check_sims_basic [-delint_dir=Dir] [-logfile=log] [-synonyms=peg_synonyms_file] NR [ < sims | - | SimsDir | Sims1 Sims2 Sims3 ...] > checked.sims [2> errors (recommended if a logfile isn't specified)]";
10    
11    $outdir  = "";
12    $logfile = "";
13    $synfile = "";
14    $trouble = 0;
15    for ($i=0; $i < @ARGV; )
16    {
17        if ($ARGV[$i] =~ m/-delint_dir=(\S+)/)
18        {
19            $outdir = $1;
20            splice @ARGV, $i, 1;
21            if (-d $outdir)
22            {
23                $trouble = 1;
24                warn "$outdir exists";
25            } else {
26                mkpath($outdir, 0, 0777) || die "Could not create $outdir";
27            }
28        }
29        elsif ($ARGV[$i] =~ m/-logfile=(\S+)/)
30        {
31            $logfile = $1;
32            splice @ARGV, $i, 1;
33            open($logfh, ">$logfile") || die "Could not write-open $logfile";
34        }
35        elsif ($ARGV[$i] =~ m/-synonyms=(\S+)/)
36        {
37            $synfile = $1;
38            splice @ARGV, $i, 1;
39            if (-s $synfile)
40            {
41                open(TMP, "<$synfile") || die "Could not read-open $synfile";
42                while (defined($entry = <TMP>))
43                {
44                    chomp $entry;
45                    $entry =~ m/^([^,]+),\d+(\S+)$/o;
46                    ($major_syn, $syns) = ($1, $2);
47                    @syns = map { m/^([^,]+)/; $1 } split /;/, $syns;
48                    foreach $syn (@syns) { $major{$syn} = $major; }
49                }
50                close(TMP) || die "Could not close $synfile";
51            }
52        }
53        elsif (-s $ARGV[$i]) {
54            ++$i;
55        }
56        else {
57            $trouble = 1;
58            print STDERR "Invalid arg $ARGV[$i]\n";
59            ++$i;
60        }
61    }
62    die "aborting due to invalid args" if ($trouble);
63    
64    (($nr = shift @ARGV) && (-s $nr))
65      || die $usage;      || die $usage;
66    
67  open(NR,"<$nr") || die $usage;  if (@ARGV == 0)
68    {
69        if (-t STDIN)
70        {
71            push @ARGV, '-';
72        }
73        else
74        {
75            print STDERR "No arguments given --- checking $FIG_Config::data/Sims by default\n";
76            push @ARGV, "$FIG_Config::data/Sims";
77        }
78    }
79    
80  $/ = "\n>";  if ((@ARGV == 1) && (-d $ARGV[0]))
 while (defined($_ = <NR>))  
81  {  {
82      chomp;      $sims_dir = shift @ARGV;
83      if ($_ =~ /^>?(\S+)[^\n]*\n(.*)/s)      opendir(SIMS, $sims_dir) || die "Could not open $sims_dir";
84        @ARGV = grep !/^\./, readdir(SIMS);
85        @ARGV = map  { $_ = "$sims_dir/$_" } @ARGV;
86        closedir(SIMS) || die "Could not close $sims_dir";
87    }
88    
89    $trouble = 0;
90    foreach $file (@ARGV)
91      {      {
92          $id  =  $1;      next if ($file eq '-');
93          $seq =  $2;      if (!-e $file) { print STDERR "Simfile $file does not exist"; $trouble = 1; }
         $seq =~ s/\s//gs;  
         $ln{$id} = length($seq);  
94      }      }
95    die "There were nonexistent input files" if $trouble;
96    
97    unless ($logfile) { $logfh = \*STDERR; }
98    unless ($outdir)  { $outfh = \*STDOUT; }
99    
100    opendir(ORGS, "$FIG_Config::organisms") || die "Could not open dir $FIG_Config::organisms";
101    @env = grep s{^(9999999\.\d+)}{$FIG_Config::organisms/$1/Features/peg/fasta}, readdir(ORGS);
102    closedir(ORGS) || die "Could not close dir $FIG_Config::organisms";
103    
104    foreach $file ($nr, @env)
105    {
106        open(TMP, "<$file") || die "Could not read-open $file";
107        print STDERR "Loading lengths from $file ...\n" if $ENV{FIG_VERBOSE};
108        while (($id, $seqP) = &FIG::read_fasta_record(\*TMP))
109        {
110            $ln{$id} = length($$seqP);
111        }
112    }
113    
114    
115    foreach $simfile (@ARGV)
116    {
117        print STDERR "Processing $simfile\n" if $ENV{FIG_VERBOSE};
118    
119        open(SIMFILE, "<$simfile") || die "Could not open $simfile";
120        if ($outdir)
121        {
122            $outfile  = "$outdir/" . basename($simfile);
123            open($outfh, ">$outfile") || die "could not write-open $outfile";
124  }  }
 $/ = "\n";  
 close(NR);  
125    
126  while (defined($_ = <STDIN>))      while (defined($sim = <SIMFILE>))
127  {  {
128      chomp;    # $/)          chomp $sim;
129      if ($_ =~ m/^\S+\t\S+\t(\d+|\d+\.\d+)(\t\d+){7}\t(\d+(\.\d*)?e[-+]?\d+|\d+\.\d+)\t(\d\.\d*e[-+]?\d+|\d+\.\d+|\d+)/)          $sim =~ s/\t\t/\t/go;
130    
131            if ($sim =~ m/^(\S+)\t(\S+)\t(\d+|\d+\.\d+)\t\d+\t\d+\t\d+\t\d+\t\d+\t\d+\t\d+\t(\d+(\.\d*)?e[-+]?\d+|\d+\.\d+)\t(\d\.\d*e[-+]?\d+|\d+\.\d+|\d+)\t(\d+)\t(\d+)/o)
132      {      {
133          # print STDERR "$1\t$2\t$3\t$4\t$5\t$6\n";              # die "$1, $2, $3, $4, $5, $6, $7, $8, $9\n";
134          ($id1,$id2,$iden,$ali_ln,$mis,$gaps,$b1,$e1,$b2,$e2,$psc,$bsc,$ln1,$ln2) =              ($id1, $id2, $ln1, $ln2) = ($1, $2, $7, $8);
135              split(/\t/,$_);              # die "$id1, $id2, $ln1, $ln2";
136    
137          if ($ln{$id1} && $ln{$id2} && ($ln{$id1} == $ln1) && ($ln{$id2} == $ln2))          if ($ln{$id1} && $ln{$id2} && ($ln{$id1} == $ln1) && ($ln{$id2} == $ln2))
138          {          {
139              print "$_\n";                  print $outfh "$sim\n";   #...print valid sims to OUTPUT
140          }          }
141          else          else
142          {          {
143              if ($ln{$id1})              if ($ln{$id1})
144              {              {
145                  if ($ln{$id1} != $ln1) { print STDERR "badlen1\t$id1\t$ln{$id1}\t$ln1\t$_\n"; }                      if ($ln{$id1} != $ln1) { print $logfh "badlen1\t$simfile, $.:\t$id1\t$ln{$id1}\t$ln1\t$sim\n"; }
146                    }
147                    else
148                    {
149                        if ($synfile)
150                        {
151                            if ($major{$id1})
152                            {
153                                print $logfh "synref1\t$simfile, $.:\t$id1\t\t\t$sim\n";
154              }              }
155              else              else
156              {              {
157                  print STDERR "undef1\t$id1\t\t\t$_\n";                              print $logfh "undef1\t$simfile, $.:\t$id1\t\t\t$sim\n";
158                            }
159                        }
160                        else
161                        {
162                            print $logfh "undef1\t$simfile, $.:\t$id1\t\t\t$sim\n";
163                        }
164              }              }
165    
166              if ($ln{$id2})              if ($ln{$id2})
167              {              {
168                  if ($ln{$id2} != $ln2) { print STDERR "badlen2\t$id2\t$ln{$id2}\t$ln2\t$_\n"; }                      if ($ln{$id2} != $ln2) { print $logfh "badlen2\t$simfile, $.:\t$id2\t$ln{$id2}\t$ln2\t$sim\n"; }
169              }              }
170              else              else
171              {              {
172                  print STDERR "undef2\t$id2\t\t\t$_\n";                      if ($synfile)
173                        {
174                            if ($major{$id2})
175                            {
176                                print $logfh "synref2\t$simfile, $.:\t$id2\t\t\t$sim\n";
177              }              }
178                            else
179                            {
180                                print $logfh "undef2\t$simfile, $.:\t$id2\t\t\t$sim\n";
181          }          }
182      }      }
183      else      else
184      {      {
185          print STDERR "INVALID FORMAT: $_\n";                          print $logfh "undef2\t$simfile, $.:\t$id2\t\t\t$sim\n";
186                        }
187                    }
188                }
189            }
190            else
191            {
192                print $logfh "INVALID FORMAT\t$simfile, $.:\t$sim\n";
193            }
194      }      }
195  }  }

Legend:
Removed from v.1.1  
changed lines
  Added in v.1.6

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3