001/*
002 *                    BioJava development code
003 *
004 * This code may be freely distributed and modified under the
005 * terms of the GNU Lesser General Public Licence.  This should
006 * be distributed with the code.  If you do not have a copy,
007 * see:
008 *
009 *      http://www.gnu.org/copyleft/lesser.html
010 *
011 * Copyright for this code is held jointly by the individual
012 * authors.  These should be listed in @author doc comments.
013 *
014 * For more information on the BioJava project and its aims,
015 * or to join the biojava-l mailing list, visit the home page
016 * at:
017 *
018 *      http://www.biojava.org/
019 *
020 */
021package org.biojava.nbio.aaproperties;
022
023import org.biojava.nbio.aaproperties.xml.AminoAcidCompositionTable;
024import org.biojava.nbio.aaproperties.xml.ElementTable;
025import org.biojava.nbio.aaproperties.xml.MyValidationEventHandler;
026import org.biojava.nbio.core.sequence.ProteinSequence;
027import org.biojava.nbio.core.sequence.compound.AminoAcidCompound;
028import org.biojava.nbio.core.sequence.compound.AminoAcidCompoundSet;
029import org.slf4j.Logger;
030import org.slf4j.LoggerFactory;
031
032import javax.xml.bind.JAXBContext;
033import javax.xml.bind.JAXBException;
034import javax.xml.bind.Unmarshaller;
035import java.io.File;
036import java.io.FileInputStream;
037import java.io.FileNotFoundException;
038import java.util.HashMap;
039import java.util.Map;
040
041/**
042 * This class contains the actual implementation of IPeptideProperties and is wrapped around by PeptideProperties for ease of use.
043 *
044 * @author kohchuanhock
045 * @version 2011.08.22
046 * @since 3.0.2
047 * @see IPeptideProperties
048 * @see PeptideProperties
049 */
050public class PeptidePropertiesImpl implements IPeptideProperties{
051
052        private final static Logger logger = LoggerFactory.getLogger(PeptidePropertiesImpl.class);
053
054        /**
055         * @return the molecular weight of water
056         */
057        private double getWaterMoleculeWeight(){
058                final double hydrogenMW = 1.0079;
059                final double hydroxideMW = 17.0073;
060                //H     1.0079  OH      17.0073
061                return hydrogenMW + hydroxideMW;
062        }
063
064        private char[] getSequence(String sequence, boolean ignoreCase){
065                if(ignoreCase){
066                        return sequence.toUpperCase().toCharArray();
067                }else{
068                        return sequence.toCharArray();
069                }
070        }
071
072        @Override
073        public double getMolecularWeight(ProteinSequence sequence) {
074                double value = 0.0;
075                AminoAcidCompoundSet aaSet = new AminoAcidCompoundSet();
076                char[] seq = getSequence(sequence.toString(), true);//ignore case
077                for(char aa:seq){
078                        AminoAcidCompound c = aaSet.getCompoundForString(String.valueOf(aa));
079                        if(Constraints.aa2MolecularWeight.containsKey(c)){
080                                value += Constraints.aa2MolecularWeight.get(c);
081                        }
082                }
083                if(value == 0)
084                        return value;
085                else
086                        return value + getWaterMoleculeWeight();
087        }
088
089        @Override
090        public double getMolecularWeight(ProteinSequence sequence, File aminoAcidCompositionFile) throws JAXBException, FileNotFoundException {
091                File elementMassFile = new File("./src/main/resources/ElementMass.xml");
092                if(!elementMassFile.exists()){
093                        throw new FileNotFoundException("Cannot locate ElementMass.xml. " +
094                                        "Please use getMolecularWeight(ProteinSequence, File, File) to specify ElementMass.xml location.");
095                }
096                return getMolecularWeightBasedOnXML(sequence, obtainAminoAcidCompositionTable(elementMassFile, aminoAcidCompositionFile));
097        }
098
099        @Override
100        public double getMolecularWeight(ProteinSequence sequence, File elementMassFile, File aminoAcidCompositionFile)
101                        throws JAXBException, FileNotFoundException{
102                return getMolecularWeightBasedOnXML(sequence, obtainAminoAcidCompositionTable(elementMassFile, aminoAcidCompositionFile));
103        }
104
105        @Override
106        public double getMolecularWeightBasedOnXML(ProteinSequence sequence, AminoAcidCompositionTable aminoAcidCompositionTable){
107                double value = 0.0;
108                char[] seq = sequence.toString().toCharArray();
109                for(char aa:seq){
110                        Double weight = aminoAcidCompositionTable.getMolecularWeight(aa);
111                        if(weight != null){
112                                value += weight;
113                        }
114                }
115                if(value == 0.0)
116                        return value;
117                else
118                        return value + getWaterMoleculeWeight();
119        }
120
121        @Override
122        public AminoAcidCompositionTable obtainAminoAcidCompositionTable(File aminoAcidCompositionFile)
123                throws JAXBException, FileNotFoundException{
124                File elementMassFile = new File("./src/main/resources/ElementMass.xml");
125                if(!elementMassFile.exists()){
126                        throw new FileNotFoundException("Cannot locate ElementMass.xml. " +
127                                        "Please use getMolecularWeight(ProteinSequence, File, File) to specify ElementMass.xml location.");
128                }
129                return obtainAminoAcidCompositionTable(elementMassFile, aminoAcidCompositionFile);
130        }
131
132        @Override
133        public AminoAcidCompositionTable obtainAminoAcidCompositionTable(File elementMassFile, File aminoAcidCompositionFile)
134                throws JAXBException, FileNotFoundException{
135                //Parse elementMassFile
136                ElementTable iTable = new ElementTable();
137                // Get a JAXB Context for the object we created above
138                JAXBContext jc = JAXBContext.newInstance(iTable.getClass());
139                Unmarshaller u = jc.createUnmarshaller();
140                u.setEventHandler(new MyValidationEventHandler());
141                iTable = (ElementTable)u.unmarshal(new FileInputStream(elementMassFile));
142                iTable.populateMaps();
143
144                //Parse aminoAcidCompositionFile
145                AminoAcidCompositionTable aTable = new AminoAcidCompositionTable();
146                // Get a JAXB Context for the object we created above
147                JAXBContext jc2 = JAXBContext.newInstance(aTable.getClass());
148                Unmarshaller u2 = jc2.createUnmarshaller();
149                u2.setEventHandler(new MyValidationEventHandler());
150                aTable = (AminoAcidCompositionTable)u2.unmarshal(new FileInputStream(aminoAcidCompositionFile));
151                aTable.computeMolecularWeight(iTable);
152                return aTable;
153        }
154
155        @Override
156        public double getExtinctionCoefficient(ProteinSequence sequence, boolean assumeCysReduced) {
157                //Tyr => Y
158                //Trp => W
159                //Cys => C
160                //E(Prot) = Numb(Tyr)*Ext(Tyr) + Numb(Trp)*Ext(Trp) + Numb(Cystine)*Ext(Cystine)
161                //where (for proteins in water measured at 280 nm): Ext(Tyr) = 1490, Ext(Trp) = 5500, Ext(Cystine) = 125;
162                AminoAcidCompoundSet aaSet = new AminoAcidCompoundSet();
163                Map<AminoAcidCompound, Integer> extinctAA2Count = this.getExtinctAACount(sequence);
164
165                double eProt;
166                if(!assumeCysReduced){
167                        eProt = extinctAA2Count.get(aaSet.getCompoundForString("Y")) *
168                                Constraints.aa2ExtinctionCoefficient.get(aaSet.getCompoundForString("Y")) +
169                                extinctAA2Count.get(aaSet.getCompoundForString("W")) *
170                                Constraints.aa2ExtinctionCoefficient.get(aaSet.getCompoundForString("W")) +
171                                extinctAA2Count.get(aaSet.getCompoundForString("C")) *
172                                Constraints.aa2ExtinctionCoefficient.get(aaSet.getCompoundForString("C"));
173                }else
174                        eProt = extinctAA2Count.get(aaSet.getCompoundForString("Y")) *
175                                Constraints.aa2ExtinctionCoefficient.get(aaSet.getCompoundForString("Y")) +
176                                extinctAA2Count.get(aaSet.getCompoundForString("W")) *
177                                Constraints.aa2ExtinctionCoefficient.get(aaSet.getCompoundForString("W"));
178
179                return eProt;
180        }
181
182        @Override
183        public double getAbsorbance(ProteinSequence sequence, boolean assumeCysReduced){
184                //Absorb(Prot) = E(Prot) / Molecular_weight
185                double mw = this.getMolecularWeight(sequence);
186                double eProt = this.getExtinctionCoefficient(sequence, assumeCysReduced);
187                if (mw == 0.0) {
188                        logger.warn("Molecular weight is 0.0, can't divide by 0: setting absorbance to 0.0");
189                        return 0.0;
190                }
191                return eProt / mw;
192        }
193
194        private Map<AminoAcidCompound, Integer> getExtinctAACount(ProteinSequence sequence){
195                //Cys => C, Tyr => Y, Trp => W
196                int numW = 0;
197                int smallW = 0;
198                double numC = 0;
199                double smallC = 0;
200                int numY = 0;
201                int smallY = 0;
202                for(char aa:sequence.getSequenceAsString().toCharArray()){
203                        switch(aa){
204                        case 'W': numW++; break;
205                        case 'w': smallW++; break;
206                        case 'C': numC += 0.5; break;
207                        case 'c': smallC += 0.5; break;
208                        case 'Y': numY++; break;
209                        case 'y': smallY++; break;
210                        }
211                }
212                AminoAcidCompoundSet aaSet = new AminoAcidCompoundSet();
213                Map<AminoAcidCompound, Integer> extinctAA2Count = new HashMap<AminoAcidCompound, Integer>();
214                //Ignore Case is always true
215                extinctAA2Count.put(aaSet.getCompoundForString("W"), numW + smallW);
216                extinctAA2Count.put(aaSet.getCompoundForString("C"), (int) (numC + smallC));
217                extinctAA2Count.put(aaSet.getCompoundForString("Y"), numY + smallY);
218                return extinctAA2Count;
219        }
220
221        @Override
222        public double getInstabilityIndex(ProteinSequence sequence) {
223                double sum = 0.0;
224                String s = sequence.getSequenceAsString().toUpperCase();
225                for(int i = 0; i < sequence.getLength() - 1; i++){
226                        String dipeptide = s.substring(i, i+2);
227                        if(Constraints.diAA2Instability.containsKey(dipeptide)){
228                                sum += Constraints.diAA2Instability.get(dipeptide);
229                        }
230                }
231                int denominator = s.length() - Utils.getNumberOfInvalidChar(s, null, true);
232
233                if (denominator==0) {
234                        logger.warn("Valid length of sequence is 0, can't divide by 0 to calculate instability index: setting instability index value to 0.0");
235                        return 0.0;
236                }
237                return sum * 10.0 / denominator;
238        }
239
240        @Override
241        public double getApliphaticIndex(ProteinSequence sequence) {
242//              Aliphatic index = X(Ala) + a * X(Val) + b * ( X(Ile) + X(Leu) )
243//              where X(Ala), X(Val), X(Ile), and X(Leu) are mole percent (100 X mole fraction)
244//              of alanine, valine, isoleucine, and leucine.
245//              The coefficients a and b are the relative volume of valine side chain (a = 2.9)
246//              and of Leu/Ile side chains (b = 3.9) to the side chain of alanine.
247//              Ala => A, Val => V, Ile => I, Leu => L
248                AminoAcidCompoundSet aaSet = new AminoAcidCompoundSet();
249                Map<AminoAcidCompound, Double> aa2Composition = getAAComposition(sequence);
250                final double a = 2.9;
251                final double b = 3.9;
252                double xAla = aa2Composition.get(aaSet.getCompoundForString("A"));
253                double xVal = aa2Composition.get(aaSet.getCompoundForString("V"));
254                double xIle = aa2Composition.get(aaSet.getCompoundForString("I"));
255                double xLeu = aa2Composition.get(aaSet.getCompoundForString("L"));
256                return (xAla + (a * xVal) + (b * (xIle + xLeu))) * 100;
257        }
258
259        @Override
260        public double getAvgHydropathy(ProteinSequence sequence) {
261                int validLength = 0;
262                double total = 0.0;
263                AminoAcidCompoundSet aaSet = new AminoAcidCompoundSet();
264                char[] seq = this.getSequence(sequence.toString(), true);
265                for(char aa:seq){
266                        AminoAcidCompound c = aaSet.getCompoundForString(String.valueOf(aa));
267                        if(Constraints.aa2Hydrophathicity.containsKey(c)){
268                                total += Constraints.aa2Hydrophathicity.get(c);
269                                validLength++;
270                        }
271                }
272                if (validLength==0) {
273                        logger.warn("Valid length of sequence is 0, can't divide by 0 to calculate average hydropathy: setting average hydropathy to 0");
274                        return 0.0;
275                }
276
277                return total / validLength;
278        }
279
280        @Override
281        public double getIsoelectricPoint(ProteinSequence sequence, boolean useExpasyValues) {
282                if(useExpasyValues){
283                        return this.getIsoelectricPointExpasy(sequence.toString().toUpperCase());
284                }else{
285                        return this.getIsoelectricPointInnovagen(sequence);
286                }
287        }
288
289        private double getIsoelectricPointInnovagen(ProteinSequence sequence){
290                double currentPH = 7.0;
291                double changeSize = 7.0;
292                String sequenceString = sequence.toString();
293                char nTerminalChar = sequenceString.charAt(0);
294                char cTerminalChar = sequenceString.charAt(sequenceString.length() - 1);
295
296                Map<AminoAcidCompound, Integer> chargedAA2Count = this.getChargedAACount(sequence);
297                double margin;
298                final double difference = 0.0001;
299
300                while(true){
301                        margin = this.getNetChargeInnovagen(chargedAA2Count, currentPH, nTerminalChar, cTerminalChar);
302                        //Within allowed difference
303                        if(margin <= difference && margin >= -difference) break;
304                        changeSize /= 2.0;
305                        if(margin > 0){
306                                currentPH += changeSize;
307                        }else{
308                                currentPH -= changeSize;
309                        }
310                }
311                return currentPH;
312        }
313
314        /*
315         *  Pseudo code obtained from email correspondance with ExPASy Helpdesk, Gregoire Rossier
316         */
317        //
318        // Table of pk values :
319        // Note: the current algorithm does not use the last two columns.
320        // Each row corresponds to an amino acid starting with Ala. J, O and U are
321        // inexistant, but here only in order to have the complete alphabet.
322        //
323        // Ct Nt Sm Sc Sn
324        //
325        private final double[][] cPk = {
326                        {3.55, 7.59, 0.0},  // A
327                        {3.55, 7.50, 0.0},  // B
328                        {3.55, 7.50, 9.00}, // C
329//                      {4.55, 7.50, 4.05}, // D
330//                      {4.75, 7.70, 4.45}, // E
331                        {3.55, 7.50, 4.05}, // D
332                        {3.55, 7.70, 4.45}, // E
333                        {3.55, 7.50, 0}, // F
334                        {3.55, 7.50, 0}, // G
335                        {3.55, 7.50, 5.98}, // H
336                        {3.55, 7.50, 0.0}, // I
337                        {0.0, 0.0, 0.0}, // J
338                        {3.55, 7.50, 10.00}, // K
339                        {3.55, 7.50, 0.0}, // L
340                        {3.55, 7.00, 0.0},// M
341                        {3.55, 7.50, 0.0},// N
342                        {0.00, 0.00, 0.0},// O
343                        {3.55, 8.36, 0.0},// P
344                        {3.55, 7.50, 0.0}, // Q
345                        {3.55, 7.50, 12.0},// R
346                        {3.55, 6.93, 0.0},// S
347                        {3.55, 6.82, 0.0}, // T
348                        {0.00, 0.00, 0.0}, // U
349                        {3.55, 7.44, 0.0},// V
350                        {3.55, 7.50, 0.0},// W
351                        {3.55, 7.50, 0.0},// X
352                        {3.55, 7.50, 10.00},// Y
353                        {3.55, 7.50, 0.0}}; // Z
354
355        private final double PH_MIN = 0.0; /* minimum pH value */
356        private final double PH_MAX = 14.0; /* maximum pH value */
357        private final double MAXLOOP = 2000.0; /* maximum number of iterations */
358        private final double EPSI = 0.0001; /* desired precision */
359
360        private double exp10(double pka){
361                return Math.pow(10, pka);
362        }
363
364        private double getIsoelectricPointExpasy(String sequence){
365                //
366                // Compute the amino-acid composition.
367                //
368                int[] comp = new int[26];
369                for(int i = 0; i < sequence.length(); i++){
370                        int index = sequence.charAt(i) - 'A';
371                        if(index < 0 || index >= 26) continue;
372                        comp[index]++;
373                }
374                //
375                // Look up N-terminal and C-terminal residue.
376                //
377                int nTermResidue = -1;
378                int index = 0;
379                while((nTermResidue < 0 || nTermResidue >= 26) && index < 25){
380                        nTermResidue = sequence.charAt(index++) - 'A';
381                }
382
383                int cTermResidue = -1;
384                index = 1;
385                while((cTermResidue < 0 || cTermResidue >= 26) && index < 25){
386                        cTermResidue = sequence.charAt(sequence.length() - index++) - 'A';
387                }
388
389                double phMin = PH_MIN;
390                double phMax = PH_MAX;
391
392                double phMid = 0.0;
393                double charge = 1.0;
394                for (int i = 0; i < MAXLOOP && (phMax - phMin) > EPSI; i++){
395                        phMid = phMin + (phMax - phMin) / 2.0;
396
397                        charge = getNetChargeExpasy(comp, nTermResidue, cTermResidue, phMid);
398
399                        if (charge > 0.0) phMin = phMid;
400                        else phMax = phMid;
401                }
402                return phMid;
403        }
404
405        @Override
406        public double getIsoelectricPoint(ProteinSequence sequence){
407                return getIsoelectricPoint(sequence, true);
408        }
409
410        @Override
411        public double getNetCharge(ProteinSequence sequence) {
412                return getNetCharge(sequence, true);
413        }
414
415        @Override
416        public double getNetCharge(ProteinSequence sequence, boolean useExpasyValues){
417                return getNetCharge(sequence, true, 7.0);
418        }
419
420        @Override
421        public double getNetCharge(ProteinSequence sequence, boolean useExpasyValues, double pHPoint){
422                if(useExpasyValues){
423                        return getNetChargeExpasy(sequence.toString().toUpperCase(), pHPoint);
424                }else{
425                        return getNetChargeInnovagen(sequence, pHPoint);
426                }
427        }
428
429        private double getNetChargeExpasy(String sequence, double pHPoint){
430                //
431                // Compute the amino-acid composition.
432                //
433                int[] comp = new int[26];
434                for(int i = 0; i < sequence.length(); i++){
435                        int index = sequence.charAt(i) - 'A';
436                        if(index < 0 || index >= 26) continue;
437                        comp[index]++;
438                }
439                //
440                // Look up N-terminal and C-terminal residue.
441                //
442                int nTermResidue = sequence.charAt(0) - 'A';
443                int cTermResidue = sequence.charAt(sequence.length() - 1) - 'A';
444                return getNetChargeExpasy(comp, nTermResidue, cTermResidue, pHPoint);
445        }
446
447        private double getNetChargeExpasy(int[] comp, int nTermResidue, int cTermResidue, double ph){
448                double cter = 0.0;
449                if(cTermResidue >= 0 && cTermResidue < 26) cter = exp10(-cPk[cTermResidue][0]) / (exp10(-cPk[cTermResidue][0]) + exp10(-ph));
450                double nter = 0.0;
451                if(nTermResidue >= 0 && nTermResidue < 26) nter = exp10(-ph) / (exp10(-cPk[nTermResidue][1]) + exp10(-ph));
452
453                double carg = comp['R' - 'A'] * exp10(-ph) / (exp10(-cPk['R' - 'A'][2]) + exp10(-ph));
454                double chis = comp['H' - 'A'] * exp10(-ph) / (exp10(-cPk['H' - 'A'][2]) + exp10(-ph));
455                double clys = comp['K' - 'A'] * exp10(-ph) / (exp10(-cPk['K' - 'A'][2]) + exp10(-ph));
456
457                double casp = comp['D' - 'A'] * exp10(-cPk['D' - 'A'][2]) / (exp10(-cPk['D' - 'A'][2]) + exp10(-ph));
458                double cglu = comp['E' - 'A'] * exp10(-cPk['E' - 'A'][2]) / (exp10(-cPk['E' - 'A'][2]) + exp10(-ph));
459
460                double ccys = comp['C' - 'A'] * exp10(-cPk['C' - 'A'][2]) / (exp10(-cPk['C' - 'A'][2]) + exp10(-ph));
461                double ctyr = comp['Y' - 'A'] * exp10(-cPk['Y' - 'A'][2]) / (exp10(-cPk['Y' - 'A'][2]) + exp10(-ph));
462
463                return (carg + clys + chis + nter) - (casp + cglu + ctyr + ccys + cter);
464        }
465
466        private double getNetChargeInnovagen(ProteinSequence sequence, double pHPoint) {
467                Map<AminoAcidCompound, Integer> chargedAA2Count = this.getChargedAACount(sequence);
468                String sequenceString = sequence.getSequenceAsString();
469                return getNetChargeInnovagen(chargedAA2Count, pHPoint, sequenceString.charAt(0), sequenceString.charAt(sequenceString.length() - 1));
470        }
471
472        private double getNetChargeInnovagen(Map<AminoAcidCompound, Integer> chargedAA2Count, double ph, char nTerminalChar, char cTerminalChar){
473                //Constraints.aa2PKa is aleady reinitialized in getChargedAACount hence no need to do it again
474
475                //Lys => K, Arg => R, His => H
476                //Asp => D, Glu => E, Cys => C, Tyr => Y
477                AminoAcidCompoundSet aaSet = new AminoAcidCompoundSet();
478
479                double nTerminalCharge = 0.0;
480                AminoAcidCompound nTermCompound = aaSet.getCompoundForString(String.valueOf(nTerminalChar));
481                if(Constraints.aa2NTerminalPka.containsKey(nTermCompound)){
482                        nTerminalCharge = this.getPosCharge(Constraints.aa2NTerminalPka.get(nTermCompound), ph);
483                }
484
485                double cTerminalCharge = 0.0;
486                AminoAcidCompound cTermCompound = aaSet.getCompoundForString(String.valueOf(cTerminalChar));
487                if(Constraints.aa2CTerminalPka.containsKey(cTermCompound)){
488                        cTerminalCharge = this.getNegCharge(Constraints.aa2CTerminalPka.get(cTermCompound), ph);
489                }
490
491                double kCharge = chargedAA2Count.get(aaSet.getCompoundForString("K")) * this.getPosCharge(Constraints.aa2PKa.get(aaSet.getCompoundForString("K")), ph);
492                double rCharge = chargedAA2Count.get(aaSet.getCompoundForString("R")) * this.getPosCharge(Constraints.aa2PKa.get(aaSet.getCompoundForString("R")), ph);
493                double hCharge = chargedAA2Count.get(aaSet.getCompoundForString("H")) * this.getPosCharge(Constraints.aa2PKa.get(aaSet.getCompoundForString("H")), ph);
494                double dCharge = chargedAA2Count.get(aaSet.getCompoundForString("D")) * this.getNegCharge(Constraints.aa2PKa.get(aaSet.getCompoundForString("D")), ph);
495                double eCharge = chargedAA2Count.get(aaSet.getCompoundForString("E")) * this.getNegCharge(Constraints.aa2PKa.get(aaSet.getCompoundForString("E")), ph);
496                double cCharge = chargedAA2Count.get(aaSet.getCompoundForString("C")) * this.getNegCharge(Constraints.aa2PKa.get(aaSet.getCompoundForString("C")), ph);
497                double yCharge = chargedAA2Count.get(aaSet.getCompoundForString("Y")) * this.getNegCharge(Constraints.aa2PKa.get(aaSet.getCompoundForString("Y")), ph);
498//              if((kCharge + rCharge + hCharge) == 0.0 && (dCharge + eCharge + cCharge + yCharge) == 0.0){
499//                      return 0.0;
500//              }
501                return (nTerminalCharge + kCharge + rCharge + hCharge) - (dCharge + eCharge + cCharge + yCharge + cTerminalCharge);
502        }
503
504        private double getPosCharge(double pka, double ph){
505                return Math.pow(10, pka) / (Math.pow(10, pka) + Math.pow(10, ph));
506        }
507
508        private double getNegCharge(double pka, double ph){
509                return Math.pow(10, ph) / (Math.pow(10, pka) + Math.pow(10, ph));
510        }
511
512        private Map<AminoAcidCompound, Integer> getChargedAACount(ProteinSequence sequence){
513                //Lys => K, Arg => R, His => H
514                //Asp => D, Glu => E, Cys => C, Tyr => Y
515                int numK = 0;
516                int numR = 0;
517                int numH = 0;
518                int numD = 0;
519                int numE = 0;
520                int numC = 0;
521                int numY = 0;
522                char[] seq = this.getSequence(sequence.getSequenceAsString(), true);
523                for(char aa:seq){
524                        switch(aa){
525                        case 'K': numK++; break;
526                        case 'R': numR++; break;
527                        case 'H': numH++; break;
528                        case 'D': numD++; break;
529                        case 'E': numE++; break;
530                        case 'C': numC++; break;
531                        case 'Y': numY++; break;
532                        }
533                }
534                AminoAcidCompoundSet aaSet = new AminoAcidCompoundSet();
535                Map<AminoAcidCompound, Integer> chargedAA2Count = new HashMap<AminoAcidCompound, Integer>();
536                chargedAA2Count.put(aaSet.getCompoundForString("K"), numK);
537                chargedAA2Count.put(aaSet.getCompoundForString("R"), numR);
538                chargedAA2Count.put(aaSet.getCompoundForString("H"), numH);
539                chargedAA2Count.put(aaSet.getCompoundForString("D"), numD);
540                chargedAA2Count.put(aaSet.getCompoundForString("E"), numE);
541                chargedAA2Count.put(aaSet.getCompoundForString("C"), numC);
542                chargedAA2Count.put(aaSet.getCompoundForString("Y"), numY);
543                return chargedAA2Count;
544        }
545
546        @Override
547        public double getEnrichment(ProteinSequence sequence, AminoAcidCompound aminoAcidCode) {
548                double counter = 0.0;
549                char[] seq = this.getSequence(sequence.getSequenceAsString(), true);
550                for(char aa:seq){
551                        if(aminoAcidCode.getShortName().equals(String.valueOf(aa))){
552                                counter++;
553                        }
554                }
555                return counter/sequence.getLength();
556        }
557
558        @Override
559        public Map<AminoAcidCompound, Double> getAAComposition(ProteinSequence sequence) {
560                int validLength = 0;
561                Map<AminoAcidCompound, Double> aa2Composition = new HashMap<AminoAcidCompound, Double>();
562                AminoAcidCompoundSet aaSet = new AminoAcidCompoundSet();
563                for(AminoAcidCompound aa:aaSet.getAllCompounds()){
564                        aa2Composition.put(aa, 0.0);
565                }
566                char[] seq = this.getSequence(sequence.toString(), true);
567                for(char aa:seq){
568                        if(PeptideProperties.standardAASet.contains(aa)){
569                                AminoAcidCompound compound = aaSet.getCompoundForString(String.valueOf(aa));
570                                aa2Composition.put(compound, aa2Composition.get(compound) + 1.0);
571                                validLength++;
572                        }
573                }
574                if(validLength > 0){
575                        for(AminoAcidCompound aa:aaSet.getAllCompounds()){
576                                aa2Composition.put(aa, aa2Composition.get(aa) / validLength);
577                        }
578                }else{
579                        for(AminoAcidCompound aa:aaSet.getAllCompounds()){
580                                aa2Composition.put(aa, 0.0);
581                        }
582                }
583                return aa2Composition;
584        }
585}