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

Annotation of /FigKernelScripts/compute_mw_and_pi.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : mkubal 1.1 # -*- perl -*-
2 :     use FIG;
3 :    
4 : mkubal 1.2 # usage: compute_mw_and_pi genome
5 :     # The results are written in attribute format to /vol/seed-attributes/molecular_weight/$genome_molecular_weight_attributes.txt and
6 :     # /vol/seed-attributes/isoelectric_point/$genome_isoelectric_point_attributes.txt
7 : mkubal 1.1
8 :     $fig = new FIG;
9 :    
10 : mkubal 1.2 my $genome,$name;
11 : mkubal 1.1
12 : mkubal 1.2 $genome = shift(@ARGV);
13 :    
14 :     if(!$genome){
15 :     print "please supply SEED genome id as argumnets\n";
16 : mkubal 1.1 exit;
17 :     }
18 :    
19 : mkubal 1.2 open(MW,">/vol/seed-attributes/molecular_weight/$genome_molecular_weight_attributes.txt");
20 :     open(ISO,">/vol/seed-attributes/isoelectric_point/$genome_isoelectric_point_attributes.txt");
21 : mkubal 1.1
22 :     # moleculare weight or amino acids
23 :     our %aaweights = ('A' => 89.09,
24 :     'C' => 121.16,
25 :     'D' => 133.10,
26 :     'E' => 147.13,
27 :     'F' => 165.19,
28 :     'G' => 75.07,
29 :     'H' => 155.16,
30 :     'I' => 131.17,
31 :     'K' => 146.19,
32 :     'L' => 131.17,
33 :     'M' => 149.21,
34 :     'N' => 132.12,
35 :     'P' => 115.13,
36 :     'Q' => 146.15,
37 :     'R' => 174.20,
38 :     'S' => 105.09,
39 :     'T' => 119.12,
40 :     'V' => 117.15,
41 :     'W' => 204.22,
42 :     'Y' => 181.19,
43 :     '*' => 0,
44 :     'X' => 0);
45 :    
46 :     our $h2oweight = 18.016;
47 :    
48 :     our %aa3letternames = ('A' => 'Ala',
49 :     'B' => 'Asx',
50 :     'C' => 'Cys',
51 :     'D' => 'Asp',
52 :     'E' => 'Glu',
53 :     'F' => 'Phe',
54 :     'G' => 'Gly',
55 :     'H' => 'His',
56 :     'I' => 'Ile',
57 :     'K' => 'Lys',
58 :     'L' => 'Leu',
59 :     'M' => 'Met',
60 :     'N' => 'Asn',
61 :     'P' => 'Pro',
62 :     'Q' => 'Gln',
63 :     'R' => 'Arg',
64 :     'S' => 'Ser',
65 :     'T' => 'Thr',
66 :     'V' => 'Val',
67 :     'W' => 'Trp',
68 :     'Y' => 'Tyr',
69 :     '*' => 'Stop',
70 :     'Z' => 'Glx',
71 :     'X' => '?');
72 :    
73 : mkubal 1.2
74 :     @pegs = $fig->pegs_of($genome);
75 :     foreach my $peg (@pegs){
76 :     my $sequence = $fig->get_translation($peg);
77 :    
78 :     my $mw = &calc_molweight($sequence);
79 :     if($mw =~/(\d+\.\d{2})/){$mw = $1}
80 :     my $pI = &calc_pI($sequence);
81 :     if($pI =~/(\d+\.\d{2})/){$pI = $1}
82 :     print MW "$peg\tmolecular_weight\t$mw\n";
83 :     print ISO "$peg\tisoelectric_point\t$pI\n";
84 : mkubal 1.1 }
85 :    
86 :     ########################################################
87 :     # calculate the molecular weight of a protein sequence #
88 :     ########################################################
89 :     sub calc_molweight {
90 :     my ($aaseq)=@_;
91 :    
92 :     my $molwt = 0;
93 :     my $numaa = length ($aaseq);
94 :    
95 :     my $i = 0;
96 :     my $aa;
97 :     while ($aa = substr($aaseq, $i, 1)) {
98 :     $molwt += $aaweights{$aa};
99 :     $i++;
100 :     }
101 :    
102 :     # Subtract a water for each peptide bond.
103 :     $molwt -= (($numaa - 1) * $h2oweight );
104 :    
105 :     return ($molwt);
106 :     };
107 :    
108 :    
109 :     ###############################
110 :     # calculate isoelectric point #
111 :     ###############################
112 :     sub calc_pI {
113 :     ###################################################################
114 :     # Assume that no electrostatic interactions perturb ionization.
115 :     # Assume that the pI lies between pH1 and pH13.
116 :     #
117 :     # Net Charge = number of positively charges residues
118 :     # minus the number of negatively charged residues
119 :     # plus the number of protonated amino termini
120 :     # minus the number of deprotonated termini
121 :     #
122 :     # For each amino acid capable of charge, the number
123 :     # of protonated residues is determined by the following equation:
124 :     #
125 :     # Np = Nt * [H+] / ([H+] + Kn)
126 :     #
127 :     # where Np is the number of protonated residues, Nt is the number of
128 :     # residues of a specific amino acid, [H+] is the H+ concentration,
129 :     # and Kn is the dissociation constant of the amino acid.
130 :     #
131 :     # Optional Input:
132 :     # picalc (PIVAR,AMINOTERMINI,CARBOXYLTERMINI)
133 :     #
134 :     # where PIVAR is the max. variance in the pI permitted (default=0.000001).
135 :     # AMINOTERMINI is the number of amino termini (default=1).
136 :     # CARBOXYLTERMINI is the number of carboxyl termini (default=1).
137 :     ###################################################################
138 :    
139 :     my ($aaseq)=@_;
140 :    
141 :     my $allowederror = 0.00001;
142 :    
143 :     my $numaterm;
144 :     my $numcterm;
145 :     if (defined($_[1])) { $numaterm = $_[1]; }
146 :     else { $numaterm = 1; }
147 :    
148 :     if (defined($_[2])) { $numcterm = $_[2]; }
149 :     else { $numcterm = 1; }
150 :    
151 :     # Count lysines, arginines, histidines,
152 :     # glutamates, aspartates,
153 :     # cysteines, tyrosines, amino termini, carboxyl termini, ambiguities.
154 :    
155 :     $_ = $aaseq;
156 :     my $numR = tr/R/R/;
157 :     my $numK = tr/K/K/;
158 :     my $numH = tr/H/H/;
159 :     my $numY = tr/Y/Y/;
160 :     my $numC = tr/C/C/;
161 :     my $numE = tr/E/E/;
162 :     my $numD = tr/D/D/;
163 :    
164 :     my $numX = tr/X/X/;
165 :     my $numN = tr/N/N/;
166 :    
167 :     ## A quick search method for the pI will involve a method that requires
168 :     ## a low number of loops through the function. Therefore, an interval
169 :     ## search method will be used to narrow the search range (an interval where
170 :     ## the net charge goes from negative to positive). Then, a binary search
171 :     ## method will be used to search for the 'netcharge=0' point within this
172 :     ## interval. Perhaps a Newton search method would be faster? The desired
173 :     ## error limit is within +/- 0.01 for the pI.
174 :    
175 :     # Set Ka values of charged amino acids:
176 :    
177 :     my $KaR = 10**(-12.48);
178 :     my $KaK = 10**(-10.53);
179 :     my $KaH = 10**(-6);
180 :     my $KaY = 10**(-10.07);
181 :     my $KaC = 10**(-10.28);
182 :     my $KaE = 10**(-4.25);
183 :     my $KaD = 10**(-3.65);
184 :     my $Kaaterm = 10**(-8.56);
185 :     my $Kacterm = 10**(-3.56);
186 :    
187 :     # define helper functions
188 :     # since both methods use a lot of variables defined in this function,
189 :     # both are implemented as coderefs
190 :    
191 :     ##################################################
192 :     # calculate the netcharge #
193 :     # Input: proton concentration #
194 :     # Output: net charge of the amino acid sequence. #
195 :     ##################################################
196 :     my $netcharge = sub {
197 :     my ($in)=@_;
198 :    
199 :     $numR*$in/($in+$KaR) + $numK*$in/($in+$KaK) + $numH*$in/($in+$KaH) - $numY + $numY*$in/($in+$KaY) - $numC + $numC*$in/($in+$KaC) - $numE + $numE*$in/($in+$KaE) - $numD + $numD*$in/($in+$KaD) + $numaterm*$in/($in+$Kaaterm) - $numcterm + $numcterm*$in/($in+$Kacterm);
200 :    
201 :     };
202 :    
203 :     ########################################################################
204 :     # perform binary search #
205 :     # Input: $high and $low boundaries #
206 :     # Exit: when the net charge at $mid is 0 or within the allowed error. #
207 :     ########################################################################
208 :     my $binarysearch;
209 :    
210 :     $binarysearch = sub {
211 :     my ($a,$b, $lastmid)=@_;
212 :    
213 :     my $mid = ($b + $a) / 2;
214 :     my $net = &$netcharge(10**(-$mid));
215 :    
216 :     if ( ($net == 0) || (($mid - $lastmid)**2 <= $allowederror) ) {
217 :     return $mid;
218 :     }
219 :    
220 :     if ($net > 0) {
221 :     &$binarysearch($a, $mid, $mid);
222 :     }
223 :     else {
224 :     &$binarysearch($mid, $b, $mid);
225 :     }
226 :     };
227 :    
228 :    
229 :     my $i;
230 :     for ($i=2; $i<15; $i++) {
231 :     if (&$netcharge(10**(-$i)) <= 0) { last; }
232 :     }
233 :    
234 :     my $high = $i;
235 :     my $low = $i - 1;
236 :     my $pI;
237 :     if (&$netcharge(10**(-$high)) == 0) {
238 :     $pI = $high;
239 :     }
240 :     elsif (&$netcharge(10**(-$low)) == 0) {
241 :     $pI = $low;
242 :     }
243 :     else {
244 :     $pI = &$binarysearch($high, $low, 0);
245 :     }
246 :    
247 :     return ($pI);
248 :     };
249 :    
250 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3