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

Annotation of /FigKernelScripts/pg_get_codon_usage_stats.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 gjoseqlib;
4 :     use Getopt::Long;
5 :     use SeedEnv;
6 :    
7 :     my $usage = "usage: pg_get_codon_usage -d Data\n";
8 :     my $dataD;
9 :     my $rc = GetOptions('d=s' => \$dataD,);
10 :    
11 :     if ((! $rc) || (! -d $dataD)) { print STDERR $usage; exit }
12 :    
13 :     if (! -d "$dataD/CodonUsage")
14 :     {
15 :     mkdir("$dataD/CodonUsage",0777);
16 :     }
17 :     my $funcs = &load_funcs($dataD);
18 :     foreach $_ (`cut -f3,4 $dataD/genomes.with.job.and.genomeID`)
19 :     {
20 :     chop;
21 :     my($job,$genome) = split(/\t/,$_);
22 :     print STDERR "job=$job genome=$genome\n";
23 :     if (! -s "$dataD/CodonUsage/$genome")
24 :     {
25 :     mkdir("$dataD/CodonUsage/$genome",0777);
26 :     my $contigs = &load_contigs($job,$genome);
27 :     my $tmp = "tmp.$$.codon.usage";
28 :     open(TMP,">$tmp") || die "could not open $tmp";
29 :     foreach $_ (`cat /vol/rast-prod/jobs/$job/rp/$genome/Features/peg/tbl`)
30 :     {
31 :     if ($_ =~ /^(\S+)\t(\S+)/)
32 :     {
33 :     my($peg,$loc) = ($1,$2);
34 :     if ($loc =~ /^(\S+)_(\d+)_(\d+)$/)
35 :     {
36 :     my($contig,$beg,$end) = ($1,$2,$3);
37 :     my $dna;
38 :     my $contig_dna = $contigs->{$contig};
39 :     if ($contig_dna && ($beg > 0) && ($end > 0) && ($beg < length($contig_dna)) && ($end < length($contig_dna)))
40 :     {
41 :     if ($beg < $end)
42 :     {
43 :     $dna = substr($contig_dna,$beg-1,($end+1-$beg));
44 :     }
45 :     else
46 :     {
47 :     $dna = &SeedUtils::rev_comp(substr($contig_dna,$end-1,($beg+1-$end)));
48 :     }
49 :     my $f = $funcs->{$peg};
50 :     print TMP ">$peg $f\n$dna\n";
51 :     }
52 :     }
53 :     }
54 :     }
55 :     close(TMP);
56 :     if (-s $tmp)
57 :     {
58 :     &SeedUtils::run("match_to_axes -f $tmp > $dataD/CodonUsage/$genome/$genome.html");
59 :     }
60 :     unlink($tmp);
61 :     }
62 :     }
63 :    
64 :     sub load_contigs {
65 :     my($job,$genome) = @_;
66 :    
67 :     my $file = "/vol/rast-prod/jobs/$job/rp/$genome/contigs";
68 :     my @seqs = &gjoseqlib::read_fasta($file);
69 :     my %contigs = map { ($_->[0] => $_->[2]) } @seqs;
70 :     return \%contigs;
71 :     }
72 :    
73 :     sub load_funcs {
74 :     my($dataD) = @_;
75 :    
76 :     my $to_func = {};
77 :    
78 :     foreach my $job (`cut -f 3 $dataD/genomes.with.job`)
79 :     {
80 :     chop $job;
81 :     foreach $_ (`cat /vol/rast-prod/jobs/$job/rp/*/proposed*functions`)
82 :     {
83 :     if ($_ =~ /^(\S+)\t(\S.*\S)/)
84 :     {
85 :     $to_func->{$1} = $2;
86 :     }
87 :     }
88 :     }
89 :     return $to_func;
90 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3