X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fanalysis%2FAAFrequency.java;h=2c1a998a11424efa84f470ab46d8dc8e3ebecb07;hb=a8f483d04205bb8273ee311c12968b7e86d205fa;hp=02569efd6e94d6dcddb74ea7fa72f395b246324c;hpb=99c58ee0ae2a848f982552e53feaf6d5cb9925e5;p=jalview.git diff --git a/src/jalview/analysis/AAFrequency.java b/src/jalview/analysis/AAFrequency.java index 02569ef..2c1a998 100755 --- a/src/jalview/analysis/AAFrequency.java +++ b/src/jalview/analysis/AAFrequency.java @@ -1,412 +1,371 @@ -/* -* Jalview - A Sequence Alignment Editor and Viewer -* Copyright (C) 2005 AM Waterhouse, J Procter, G Barton, M Clamp, S Searle -* -* This program is free software; you can redistribute it and/or -* modify it under the terms of the GNU General Public License -* as published by the Free Software Foundation; either version 2 -* of the License, or (at your option) any later version. -* -* This program is distributed in the hope that it will be useful, -* but WITHOUT ANY WARRANTY; without even the implied warranty of -* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -* GNU General Public License for more details. -* -* You should have received a copy of the GNU General Public License -* along with this program; if not, write to the Free Software -* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA -*/ - -package jalview.analysis; - -import jalview.jbgui.*; -import jalview.datamodel.*; -import jalview.io.*; -import jalview.analysis.*; - -import java.awt.*; -import java.applet.Applet; -import java.util.*; -import java.net.*; -import java.io.*; - -public class AAFrequency { - - // Takes in a vector of sequences and column start and column end - // and returns a vector of size (end-start+1). Each element of the - // vector contains a hashtable with the keys being residues and - // the values being the count of each residue in that column. - // This class is used extensively in calculating alignment colourschemes - // that depend on the amount of conservation in each alignment column. - - - public static Vector calculate(Vector sequences,int start,int end) { - - Vector result = new Vector(); - - for (int i = start;i <= end; i++) - { - - Hashtable residueHash = new Hashtable(); - int maxCount = 0; - String maxResidue = "-"; - int nongap = 0; - for (int j=0; j < sequences.size(); j++) - { - - if (sequences.elementAt(j) instanceof Sequence) - { - Sequence s = (Sequence)sequences.elementAt(j); - - if (s.getSequence().length() > i) - { - - String res = s.getSequence().charAt(i)+""; - - if (!jalview.util.Comparison.isGap(res.charAt(0))) - nongap++; - else - res = "-"; // we always use this for gaps in the property vectors - - if (residueHash.containsKey(res)) - { - - int count = ((Integer)residueHash.get(res)).intValue() ; - count++; - - if (!jalview.util.Comparison.isGap(res.charAt(0)) && count >= maxCount) - { - - if(count>maxCount) - maxResidue = res; - else if(maxResidue.indexOf(res)==-1) - maxResidue += res; - - maxCount = count; - } - - residueHash.put(res,new Integer(count)); - } - else - residueHash.put(res,new Integer(1)); - - - } - else - { - if (residueHash.containsKey("-")) - { - int count = ((Integer)residueHash.get("-")).intValue() ; - count++; - residueHash.put("-",new Integer(count)); - } - else - residueHash.put("-",new Integer(1)); - - } - } - } - - residueHash.put("maxCount",new Integer(maxCount)); - if(maxCount<0) - System.out.println("asasa "+maxCount); - residueHash.put("maxResidue", maxResidue); - residueHash.put("size", new Integer(sequences.size())); - residueHash.put("nongap", new Integer(nongap)); - result.addElement(residueHash); - } - - return result; - } - - public static Vector calculatePID(SequenceI refseq,Vector sequences,int window,int start,int end) { - - Vector result = new Vector(); - - boolean init = true; - - - Vector prev = null; - - for (int i = start;i <= end; i++) { - Vector values = new Vector(); - - result.addElement(values); - - // If start < window/2 then set value to zero. - - if (i< window/2 || i >= refseq.getSequence().length()-window/2) { - for (int j = 0; j < sequences.size(); j++) { - values.addElement(new Integer(0)); - } - } else if (init == true) { - init = false; - - int winstart = i-window/2; - int winend = i+window/2; - - if (window%2 != 0) { - winend++; - } - - for (int j = 0; j < sequences.size(); j++) { - values.addElement(new Integer(0)); - } - - for (int k = winstart; k <= winend; k++) { - String refchar = refseq.getSequence().substring(k,k+1); - if (jalview.util.Comparison.isGap(refchar.charAt(0))) - refchar="-"; - else { - for (int j = 0; j < sequences.size(); j++) { - - Sequence s = (Sequence)sequences.elementAt(j); - - if (s.getSequence().length() > k) { - - String res = s.getSequence().substring(k,k+1); // no gapchar test needed - - if (res.equals(refchar)) { - int val = ((Integer)values.elementAt(j)).intValue(); - val++; - values.setElementAt(new Integer(val),j); - } - } - } - } - - } - - prev = values; - } else { - int winstart = i-window/2; - int winend = i+window/2; - - if (window%2 != 0) { - winend++; - } - // We need to take the previous set of values - // subtract the pid at winstart-1 - // and add the pid at winend; - - String pre_refchar = refseq.getSequence().substring(winstart-1,winstart); - String pos_refchar = "-"; - - if (refseq.getSequence().length() > winend) { - pos_refchar = refseq.getSequence().substring(winend,winend+1); - } - - for (int j = 0; j < sequences.size(); j++) { - // First copy the pid value from i-1 - - int val = ((Integer)prev.elementAt(j)).intValue(); - - Sequence s = (Sequence)sequences.elementAt(j); - - String pre_char = s.getSequence().substring(winstart-1,winstart); - - String pos_char = "-"; - - if (s.getSequence().length() > winend) { - pos_char = s.getSequence().substring(winend,winend+1); - } - - // Now substract 1 if the chars at winstart-1 match - - if (jalview.util.Comparison.isGap(pre_refchar.charAt(0)) == false - && pre_char.equals(pre_refchar)) { - val--; - } - - if (jalview.util.Comparison.isGap(pos_refchar.charAt(0)) == false - && pos_char.equals(pos_refchar)) { - val++; - } - - values.addElement(new Integer(val)); - - - } - prev = values; - } - } - - return result; - } - - public static Hashtable findBlocks(Vector seqs, int start, int end,Vector exc) { - - // start and end are in real (not relative coords); - - // The coords in the hashtable that is returned are in relative coords - // i.e. start from 0 - - Hashtable blocks = new Hashtable(); - - boolean prev = false; - int bstart = -1; - - for (int i = start; i <= end ; i++) { - SequenceI seq = (SequenceI)seqs.elementAt(0); - - char c = seq.getCharAt(i); - - boolean found = true; - - int j = 1; - - while (j < seqs.size() && found == true) { - - SequenceI jseq = (SequenceI)seqs.elementAt(j); - - if (!exc.contains(jseq)) { - - char cc = jseq.getCharAt(i); - - if ( cc != c) { - found = false; - } - } - j++; - } - - - if (prev == false && found == true) { - bstart = i; - } else if (prev == true && found == false && bstart != -1) { - - int blockstart = bstart-start; - int blocklen = i-bstart; - - //System.out.println("Start len " + blockstart + " " + blocklen); - - for (int jj = blockstart; jj < blockstart + blocklen;jj++) { - blocks.put(new Integer(jj),new Integer(blocklen)); - } - - bstart = -1; - } - prev = found; - } - - if (bstart != -1) { - - int blockstart = bstart-start; - int blocklen = end-bstart; - - // System.out.println("Start len " + blockstart + " " + blocklen); - - for (int jj = blockstart; jj < blockstart + blocklen;jj++) { - blocks.put(new Integer(blockstart),new Integer(blocklen)); - } - - } - return blocks; - } - - - - public static Hashtable findKmerCount(SequenceI seq, int start, int end,int window, int step,Vector kmers) { - - int tmpstart = start; - Hashtable vals = new Hashtable(); - - while (tmpstart <= end) { - - String tmpstr = seq.getSequence().substring(tmpstart-window/2,tmpstart+window/2); - - int count = 0; - - //System.out.println("Str " + tmpstr); - - for (int ii = 0; ii < kmers.size(); ii++) { - String kmer = ((SequenceI)kmers.elementAt(ii)).getSequence(); - - int i = -1; - - while (tmpstr.indexOf(kmer,i) != -1) { - i = tmpstr.indexOf(kmer,i); - - i++; - count++; - } - ii++; - } - vals.put(new Integer(tmpstart),new Integer(count)); - tmpstart += step; - } - return vals; - } - - public static Hashtable findBlockStarts(Vector seqs, int start, int end,Vector exc) { - - // start and end are in real (not relative coords); - - // The coords in the hashtable that is returned are in relative coords - // i.e. start from 0 - - Hashtable blocks = new Hashtable(); - - boolean prev = false; - int bstart = -1; - - for (int i = start; i <= end ; i++) { - SequenceI seq = (SequenceI)seqs.elementAt(0); - - char c = seq.getCharAt(i); - - boolean found = true; - - int j = 1; - - while (j < seqs.size() && found == true) { - - SequenceI jseq = (SequenceI)seqs.elementAt(j); - - if (!exc.contains(jseq)) { - - char cc = jseq.getCharAt(i); - - if ( cc != c) { - found = false; - } - } - j++; - } - - - if (prev == false && found == true) { - bstart = i; - } else if (prev == true && found == false && bstart != -1) { - - int blockstart = bstart-start; - int blocklen = i-bstart; - - // System.out.println("Start len " + blockstart + " " + blocklen); - - //for (int jj = blockstart; jj < blockstart + blocklen;jj++) { - blocks.put(new Integer(blockstart),new Integer(blocklen)); - // } - - bstart = -1; - } - prev = found; - } - - if (bstart != -1) { - - int blockstart = bstart-start; - int blocklen = end-bstart; - - // System.out.println("Start len " + blockstart + " " + blocklen); - - //for (int jj = blockstart; jj < blockstart + blocklen;jj++) { - blocks.put(new Integer(blockstart),new Integer(blocklen)); - // } - - } - return blocks; - } - -} - +/* + * Jalview - A Sequence Alignment Editor and Viewer (Version 2.8.2) + * Copyright (C) 2014 The Jalview Authors + * + * This file is part of Jalview. + * + * Jalview is free software: you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. + * + * Jalview is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty + * of MERCHANTABILITY or FITNESS FOR A PARTICULAR + * PURPOSE. See the GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License along with Jalview. If not, see . + * The Jalview Authors are detailed in the 'AUTHORS' file. + */ +package jalview.analysis; + +import java.util.*; + +import jalview.util.Format; +import jalview.datamodel.*; + +/** + * Takes in a vector or array of sequences and column start and column end and + * returns a new Hashtable[] of size maxSeqLength, if Hashtable not supplied. + * This class is used extensively in calculating alignment colourschemes that + * depend on the amount of conservation in each alignment column. + * + * @author $author$ + * @version $Revision$ + */ +public class AAFrequency +{ + // No need to store 1000s of strings which are not + // visible to the user. + public static final String MAXCOUNT = "C"; + + public static final String MAXRESIDUE = "R"; + + public static final String PID_GAPS = "G"; + + public static final String PID_NOGAPS = "N"; + + public static final String PROFILE = "P"; + + public static final Hashtable[] calculate(List list, + int start, int end) + { + return calculate(list, start, end, false); + } + + public static final Hashtable[] calculate(List sequences, + int start, int end, boolean profile) + { + SequenceI[] seqs = new SequenceI[sequences.size()]; + int width = 0; + synchronized (sequences) + { + for (int i = 0; i < sequences.size(); i++) + { + seqs[i] = sequences.get(i); + if (seqs[i].getLength() > width) + { + width = seqs[i].getLength(); + } + } + + Hashtable[] reply = new Hashtable[width]; + + if (end >= width) + { + end = width; + } + + calculate(seqs, start, end, reply, profile); + return reply; + } + } + + public static final void calculate(SequenceI[] sequences, int start, + int end, Hashtable[] result) + { + calculate(sequences, start, end, result, false); + } + + public static final void calculate(SequenceI[] sequences, int start, + int end, Hashtable[] result, boolean profile) + { + Hashtable residueHash; + int maxCount, nongap, i, j, v, jSize = sequences.length; + String maxResidue; + char c='-'; + float percentage; + + int[] values = new int[255]; + + char[] seq; + + for (i = start; i < end; i++) + { + residueHash = new Hashtable(); + maxCount = 0; + maxResidue = ""; + nongap = 0; + values = new int[255]; + + for (j = 0; j < jSize; j++) + { + if (sequences[j] == null) + { + System.err + .println("WARNING: Consensus skipping null sequence - possible race condition."); + continue; + } + seq = sequences[j].getSequence(); + if (seq.length > i) + { + c = seq[i]; + + if (c == '.' || c == ' ') + { + c = '-'; + } + + if (c == '-') + { + values['-']++; + continue; + } + else if ('a' <= c && c <= 'z') + { + c -= 32; // ('a' - 'A'); + } + + nongap++; + values[c]++; + + } + else + { + values['-']++; + } + } + if (jSize==1) + { + maxResidue = String.valueOf(c); + maxCount=1; + } else {for (v = 'A'; v < 'Z'; v++) + { + if (values[v] < 2 || values[v] < maxCount) + { + continue; + } + + if (values[v] > maxCount) + { + maxResidue = String.valueOf((char) v); + } + else if (values[v] == maxCount) + { + maxResidue += String.valueOf((char) v); + } + maxCount = values[v]; + } + } + if (maxResidue.length() == 0) + { + maxResidue = "-"; + } + if (profile) + { + residueHash.put(PROFILE, new int[][] + { values, new int[] + { jSize, nongap } }); + } + residueHash.put(MAXCOUNT, new Integer(maxCount)); + residueHash.put(MAXRESIDUE, maxResidue); + + percentage = ((float) maxCount * 100) / jSize; + residueHash.put(PID_GAPS, new Float(percentage)); + + if (nongap>0) { + // calculate for non-gapped too + percentage = ((float) maxCount * 100) / nongap; + } + residueHash.put(PID_NOGAPS, new Float(percentage)); + + result[i] = residueHash; + } + } + + /** + * Compute all or part of the annotation row from the given consensus + * hashtable + * + * @param consensus + * - pre-allocated annotation row + * @param hconsensus + * @param iStart + * @param width + * @param ignoreGapsInConsensusCalculation + * @param includeAllConsSymbols + * @param nseq + */ + public static void completeConsensus(AlignmentAnnotation consensus, + Hashtable[] hconsensus, int iStart, int width, + boolean ignoreGapsInConsensusCalculation, + boolean includeAllConsSymbols, long nseq) + { + completeConsensus(consensus, hconsensus, iStart, width, + ignoreGapsInConsensusCalculation, includeAllConsSymbols, null, nseq); // new + // char[] + // { 'A', 'C', 'G', 'T', 'U' }); + } + + public static void completeConsensus(AlignmentAnnotation consensus, + Hashtable[] hconsensus, int iStart, int width, + boolean ignoreGapsInConsensusCalculation, + boolean includeAllConsSymbols, char[] alphabet, long nseq) + { + float tval, value; + if (consensus == null || consensus.annotations == null + || consensus.annotations.length < width) + { + // called with a bad alignment annotation row - wait for it to be + // initialised properly + return; + } + String fmtstr="%3.1f"; + int precision=0; + while (nseq>=10) { + precision++; + nseq/=10; + } + final Format fmt; + if (precision>1) + { + //if (precision>2) + { + fmtstr = "%"+(2+precision)+"."+(precision)+"f"; + } + fmt = new Format(fmtstr); + } else { + fmt = null; + } + for (int i = iStart; i < width; i++) + { + Hashtable hci; + if (i >= hconsensus.length || ((hci = hconsensus[i]) == null)) + { + // happens if sequences calculated over were shorter than alignment + // width + consensus.annotations[i] = null; + continue; + } + value = 0; + Float fv; + if (ignoreGapsInConsensusCalculation) + { + fv = (Float) hci.get(AAFrequency.PID_NOGAPS); + } + else + { + fv = (Float) hci.get(AAFrequency.PID_GAPS); + } + if (fv == null) + { + consensus.annotations[i] = null; + // data has changed below us .. give up and + continue; + } + value = fv.floatValue(); + String maxRes = hci.get(AAFrequency.MAXRESIDUE).toString(); + String mouseOver = hci.get(AAFrequency.MAXRESIDUE) + " "; + if (maxRes.length() > 1) + { + mouseOver = "[" + maxRes + "] "; + maxRes = "+"; + } + int[][] profile = (int[][]) hci.get(AAFrequency.PROFILE); + if (profile != null && includeAllConsSymbols) + { + mouseOver = ""; + if (alphabet != null) + { + for (int c = 0; c < alphabet.length; c++) + { + tval = profile[0][alphabet[c]] * 100f + / profile[1][ignoreGapsInConsensusCalculation ? 1 : 0]; + mouseOver += ((c == 0) ? "" : "; ") + alphabet[c] + " " + + ((fmt!=null) ? fmt.form(tval) : ((int) tval)) + "%"; + } + } + else + { + Object[] ca = new Object[profile[0].length]; + float[] vl = new float[profile[0].length]; + for (int c = 0; c < ca.length; c++) + { + ca[c] = new char[] + { (char) c }; + vl[c] = profile[0][c]; + } + ; + jalview.util.QuickSort.sort(vl, ca); + for (int p = 0, c = ca.length - 1; profile[0][((char[]) ca[c])[0]] > 0; c--) + { + if (((char[]) ca[c])[0] != '-') + { + tval = profile[0][((char[]) ca[c])[0]] + * 100f + / profile[1][ignoreGapsInConsensusCalculation ? 1 : 0]; + mouseOver += ((p == 0) ? "" : "; ") + ((char[]) ca[c])[0] + + " " + ((fmt!=null) ? fmt.form(tval) : ((int) tval)) + "%"; + p++; + + } + } + + } + } + else + { + mouseOver += ((fmt!=null) ? fmt.form(value) : ((int) value)) + "%"; + } + consensus.annotations[i] = new Annotation(maxRes, mouseOver, ' ', + value); + } + } + + /** + * get the sorted profile for the given position of the consensus + * + * @param hconsensus + * @return + */ + public static int[] extractProfile(Hashtable hconsensus, + boolean ignoreGapsInConsensusCalculation) + { + int[] rtnval = new int[64]; + int[][] profile = (int[][]) hconsensus.get(AAFrequency.PROFILE); + if (profile == null) + return null; + Object[] ca = new Object[profile[0].length]; + float[] vl = new float[profile[0].length]; + for (int c = 0; c < ca.length; c++) + { + ca[c] = new char[] + { (char) c }; + vl[c] = profile[0][c]; + } + ; + jalview.util.QuickSort.sort(vl, ca); + rtnval[0] = 2; + rtnval[1] = 0; + for (int c = ca.length - 1; profile[0][((char[]) ca[c])[0]] > 0; c--) + { + if (((char[]) ca[c])[0] != '-') + { + rtnval[rtnval[0]++] = ((char[]) ca[c])[0]; + rtnval[rtnval[0]] = (int) (profile[0][((char[]) ca[c])[0]] * 100f / profile[1][ignoreGapsInConsensusCalculation ? 1 + : 0]); + rtnval[1] += rtnval[rtnval[0]++]; + } + } + return rtnval; + } +}