[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.2 - (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 : overbeek 1.2 while ($last && ($last =~ /^(\S[^\t]*)\t(\d+)\t(\S+)/))
51 : overbeek 1.1 {
52 :     my $fam = $1;
53 :     my $subfam = $2;
54 :     my @set;
55 : overbeek 1.2 while ($last && ($last =~ /^(\S[^\t]*\S)\t(\d+)\t(\S+)/) && ($fam eq $1) && ($subfam == $2))
56 : overbeek 1.1 {
57 :     push(@set,[$1,$2,$3]);
58 :     $last = <FAMS>;
59 :     }
60 :     push(@sets,\@set);
61 :     }
62 : overbeek 1.2 if ($last) { die "POORLY FORMATTED FAMILY: CHECK $last" }
63 : overbeek 1.1 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 : overbeek 1.2 my($mean,$std_dev) = &mean_stddev(@lens);
76 : overbeek 1.1
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 :     }
87 : overbeek 1.2
88 :     sub mean_stddev
89 :     {
90 :     my $n = @_;
91 :     return $n ? ( shift, undef) : ( undef, undef ) if $n < 2;
92 :     my ( $sum, $sum2 ) = ( 0, 0 );
93 :     foreach ( @_ ) { $sum += $_; $sum2 += $_ * $_ }
94 :     my $x = $sum / $n;
95 :     my $to_check = ( $sum2 - ($n * $x * $x )) / ( $n - 1 );
96 :     ( $x, ($to_check > 0) ? sqrt($to_check) : 0);
97 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3