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

Annotation of /FigKernelScripts/get_families_final.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (view) (download) (as text)

1 : overbeek 1.1 use strict;
2 :     use Data::Dumper;
3 :     use Getopt::Long;
4 :     use SeedEnv;
5 :     use gjoseqlib;
6 :     use gjostat;
7 :    
8 :     my $usage = "usage: get_families_final -f families -s Seq\n";
9 :     my $families;
10 :     my $seqD;
11 :     my $rc = GetOptions('f=s' => \$families,
12 :     's=s' => \$seqD);
13 :    
14 :     if ((! $rc) || (! $families) || (! $seqD))
15 :     {
16 :     print STDERR $usage; exit ;
17 :     }
18 :    
19 :     open(IN,"cat $families.good $families.bad.fixed $families.missed |")
20 :     || die "could not read families";
21 :    
22 :     my %genomes;
23 :     my %needed_pegs;
24 :    
25 :     while (defined($_ = <IN>))
26 :     {
27 :     chomp;
28 :     my(undef,undef,$peg) = split(/\t/,$_);
29 :     my $g = &SeedUtils::genome_of($peg);
30 :     $genomes{$g} = 1;
31 :     $needed_pegs{$peg} = 1;
32 :     }
33 :     close(IN);
34 :     my %to_seq;
35 :     foreach my $g (keys(%genomes))
36 :     {
37 :     my @tuples = grep { $needed_pegs{$_->[0]} } &gjoseqlib::read_fasta("$seqD/$g");
38 :     foreach my $tuple (@tuples)
39 :     {
40 :     my($peg,undef,$seq) = @$tuple;
41 :     $to_seq{$peg} = $seq;
42 :     }
43 :     }
44 :    
45 :     open(FAMS,"cat $families.good $families.bad.fixed $families.missed |")
46 :     || die "could not open families";
47 :    
48 :     my @sets;
49 :     my $last = <FAMS>;
50 :     while ($last && ($last =~ /^(\S[^\t]+\S)\t(\d+)\t(\S+)/))
51 :     {
52 :     my $fam = $1;
53 :     my $subfam = $2;
54 :     my @set;
55 :     while ($last && ($last =~ /^(\S[^\t]+\S)\t(\d+)\t(\S+)/) && ($fam eq $1) && ($subfam == $2))
56 :     {
57 :     push(@set,[$1,$2,$3]);
58 :     $last = <FAMS>;
59 :     }
60 :     push(@sets,\@set);
61 :     }
62 :    
63 :     my $famN = 1;
64 :     foreach my $set (sort { (@$b <=> @$a) or
65 :     ($a->[0] cmp $b->[0]) or
66 :     ($a->[1] <=> $b->[1])
67 :     } @sets)
68 :     {
69 :     my @lens;
70 :     foreach my $tuple (@$set)
71 :     {
72 :     my $ln = length($to_seq{$tuple->[2]});
73 :     push(@lens,$ln);
74 :     }
75 :     my($mean,$std_dev) = &gjostat::mean_stddev(@lens);
76 :    
77 :     foreach my $tuple (@$set)
78 :     {
79 :     my($fam,$subfam,$peg) = @$tuple;
80 :     my $ln = length($to_seq{$peg});
81 :     if (! $ln) { print STDERR &Dumper($peg,$to_seq{$peg}); die "HERE"; }
82 :     my $z_sc = ($ln - $mean) / ($std_dev + 0.000000001);
83 :     print join("\t",($famN,$fam,$subfam,$peg,$ln,sprintf("%0.3f",$mean),sprintf("%0.3f",$std_dev),sprintf("%0.3f",$z_sc))),"\n";
84 :     }
85 :     $famN++;
86 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3