[Bio] / FigKernelPackages / raelib.pm Repository:
ViewVC logotype

Annotation of /FigKernelPackages/raelib.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : redwards 1.1 # -*- perl -*-
2 :    
3 :     =pod
4 :    
5 :     =head1
6 :    
7 :     Some routines and things that Rob uses. Please feel free to use at will and incorporate into
8 :     your own code or move them into FIG.pm or elsewhere.
9 :    
10 :     =cut
11 :    
12 :     package raelib;
13 :     use strict;
14 :     use FIG;
15 :     my $fig=new FIG;
16 :    
17 : redwards 1.4 =head2 features_on_contig
18 :    
19 :     Returns a reference to an array containing all the features on a contig in a genome.
20 :    
21 :     use:
22 :    
23 :     my $arrayref=$rae->features_on_contig($genome, $contig);
24 :    
25 :     or
26 :    
27 :     foreach my $peg (@{$rae->features_on_contig($genome, $contig)}) {
28 :     ... blah blah ...
29 :     }
30 :    
31 :     returns undef if contig is not a part of genome or there is nothing to return, otherwise returns a list of pegs
32 :    
33 :     v. experimental and guaranteed not to work!
34 :    
35 :     =cut
36 :    
37 :     sub features_on_contig {
38 :     my ($self, $genome, $contig)=@_;
39 :     # were this in FIG.pm you'd use this line:
40 :     #my $rdbH = $self->db_handle;
41 :    
42 :     my $rdbH = $fig->db_handle;
43 :     my $relational_db_response=$rdbH->SQL('SELECT id FROM features WHERE (genome = \'' . $genome . '\' AND location ~* \'' . $contig . '\')');
44 :     # this is complicated. A reference to an array of references to arrays, and we only want the first element.
45 :     # simplify.
46 :     my @results;
47 :     foreach my $res (@$relational_db_response) {push @results, $res->[0]}
48 :     return \@results;
49 :     }
50 :    
51 :    
52 :    
53 :    
54 :    
55 :    
56 :    
57 : redwards 1.1 =head2 pirsfcorrespondance
58 :    
59 :     Generate the pirsf->fig id correspondance. This is only done once and the correspondance file is written.
60 :     This is so that we can easily go back and forth.
61 :    
62 :     The correspondance has PIR ID \t FIG ID\n, and is probably based on ftp://ftp.pir.georgetown.edu/pir_databases/pirsf/data/pirsfinfo.dat
63 :    
64 :     =cut
65 :    
66 :     sub pirsfcorrespondance {
67 :     my ($self, $from, $to)=@_;
68 :     die "File $from does not exist as called in $0" unless (-e $from);
69 :     open (IN, $from) || die "Can't open $from";
70 :     open (OUT, ">$to") || die "Can't write tot $to";
71 :     while (<IN>) {
72 :     if (/^>/) {print OUT; next}
73 :     chomp;
74 :     my $done;
75 : redwards 1.2 foreach my $peg ($fig->by_alias("uni|$_")) {
76 : redwards 1.1 print OUT $_, "\t", $peg, "\n";
77 :     $done=1;
78 :     }
79 :     unless ($done) {print OUT $_, "\t\n"}
80 :     }
81 :     close IN;
82 :     close OUT;
83 :     }
84 :    
85 :    
86 :     =head2 ss_by_id
87 :    
88 :     Generate a list of subsystems that a peg occurs in. This is a ; separated list.
89 :     This is a wrapper that removes roles and ignores essential things
90 :    
91 :     =cut
92 :    
93 :     sub ss_by_id {
94 :     my ($self, $peg)=@_;
95 :     my $ssout;
96 :     foreach my $ss (sort $fig->subsystems_for_peg($peg))
97 :     {
98 :     next if ($$ss[0] =~ /essential/i); # Ignore the Essential B-subtilis subsystems
99 :     $ssout.=$$ss[0]."; ";
100 :     }
101 :     $ssout =~ s/; $//;
102 :     return $ssout;
103 :     }
104 :    
105 : redwards 1.3 =head2 ss_by_homol
106 :    
107 :     Generate a list of subsystems that homologs of a peg occur in. This is a ; separated list.
108 :     This is also a wrapper around sims and ss, but makes everything unified
109 :    
110 :     =cut
111 :    
112 :     sub ss_by_homol {
113 :     my ($self, $peg)=@_;
114 :     return unless ($peg);
115 :     my ($maxN, $maxP)=(50, 1e-20);
116 :    
117 :     # find the sims
118 :     my @sims=$fig->sims($peg, $maxN, $maxP, 'fig');
119 :    
120 :     # we are only going to keep the best hit for each peg
121 :     # in a subsystem
122 :     my $best_ss_score; my $best_ss_id;
123 :     foreach my $sim (@sims)
124 :     {
125 :     my $simpeg=$$sim[1];
126 :     my $simscore=$$sim[10];
127 :     my @subsys=$fig->subsystems_for_peg($simpeg);
128 :     foreach my $ss (@subsys)
129 :     {
130 :     if (! defined $best_ss_score->{$$ss[0]}) {$best_ss_score->{$$ss[0]}=$simscore; $best_ss_id->{$$ss[0]}=$simpeg}
131 :     elsif ($best_ss_score->{$$ss[0]} > $simscore)
132 :     {
133 :     $best_ss_score->{$$ss[0]}=$simscore;
134 :     $best_ss_id->{$$ss[0]}=$simpeg;
135 :     }
136 :     }
137 :     }
138 :    
139 :     my $ssoutput=join "", (map {"$_ (".$best_ss_id->{$_}."), "} keys %$best_ss_id);
140 :    
141 :     $ssoutput =~ s/, $//;
142 :     return $ssoutput;
143 :     }
144 :    
145 :     =head2 tagvalue
146 :    
147 :     This will just check for tag value pairs and return either an array of values or a single ; separated list (if called as a scalar)
148 :    
149 :     e.g. $values=raelib->tagvalue($peg, "PIRSF"); print join "\n", @$values;
150 :    
151 :     Returns an empty array if no tag/value appropriate.
152 :    
153 :     Just because I use this a lot I don't want to waste rewriting it.
154 :    
155 :     =cut
156 :    
157 :     sub tagvalue {
158 :     my ($self, $peg, $tag)=@_;
159 :     my @return;
160 :     my @attr=$fig->feature_attributes($peg);
161 :     foreach my $attr (@attr) {
162 :     my ($gottag, $val, $link)=@$attr;
163 :     push @return, $val if ($gottag eq $tag);
164 :     }
165 :     return wantarray ? @return : join "; ", @return;
166 :     }
167 : redwards 1.1
168 :    
169 :    
170 :    
171 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3