*/
package jalview.analysis;
+import jalview.analysis.scoremodels.ScoreMatrix;
+import jalview.analysis.scoremodels.ScoreModels;
import jalview.datamodel.AlignmentAnnotation;
import jalview.datamodel.Annotation;
import jalview.datamodel.ResidueCount;
+import jalview.datamodel.ResidueCount.SymbolCounts;
import jalview.datamodel.Sequence;
import jalview.datamodel.SequenceI;
-import jalview.datamodel.ResidueCount.SymbolCounts;
import jalview.schemes.ResidueProperties;
import jalview.util.Comparison;
import java.util.List;
import java.util.Map;
import java.util.Map.Entry;
+import java.util.SortedMap;
import java.util.TreeMap;
import java.util.Vector;
/**
* Calculates conservation values for a given set of sequences
- *
- * @author $author$
- * @version $Revision$
*/
public class Conservation
{
+ /*
+ * need to have a minimum of 3% of sequences with a residue
+ * for it to be included in the conservation calculation
+ */
+ private static final int THRESHOLD_PERCENT = 3;
+
private static final int TOUPPERCASE = 'a' - 'A';
SequenceI[] sequences;
private String[] consSymbs;
/**
- * Creates a new Conservation object.
+ * Constructor using default threshold of 3%
*
* @param name
* Name of conservation
- * @param threshold
- * to count the residues in residueHash(). commonly used value is 3
* @param sequences
* sequences to be used in calculation
* @param start
* @param end
* end residue position
*/
+ public Conservation(String name, List<SequenceI> sequences, int start,
+ int end)
+ {
+ this(name, THRESHOLD_PERCENT, sequences, start, end);
+ }
+
+ /**
+ * Constructor
+ *
+ * @param name
+ * Name of conservation
+ * @param threshold
+ * percentage of sequences at or below which property conservation is
+ * ignored
+ * @param sequences
+ * sequences to be used in calculation
+ * @param start
+ * start column position
+ * @param end
+ * end column position
+ */
public Conservation(String name, int threshold,
List<SequenceI> sequences, int start, int end)
{
}
/**
- * Translate sequence i into a numerical representation and store it in the
- * i'th position of the seqNums array.
+ * Translate sequence i into score matrix indices and store it in the i'th
+ * position of the seqNums array.
*
* @param i
+ * @param sm
*/
- private void calcSeqNum(int i)
+ private void calcSeqNum(int i, ScoreMatrix sm)
{
- String sq = null; // for dumb jbuilder not-inited exception warning
- int[] sqnum = null;
-
+ int gapIndex = sm.getGapIndex();
int sSize = sequences.length;
if ((i > -1) && (i < sSize))
{
- sq = sequences[i].getSequenceAsString();
+ String sq = sequences[i].getSequenceAsString();
if (seqNums.size() <= i)
{
seqNums.addElement(new int[sq.length() + 1]);
}
+ /*
+ * the first entry in the array is the sequence's hashcode,
+ * following entries are matrix indices of sequence characters
+ */
if (sq.hashCode() != seqNums.elementAt(i)[0])
{
int j;
maxLength = len;
}
- sqnum = new int[len + 1]; // better to always make a new array -
+ int[] sqnum = new int[len + 1]; // better to always make a new array -
// sequence can change its length
sqnum[0] = sq.hashCode();
for (j = 1; j <= len; j++)
{
- sqnum[j] = jalview.schemes.ResidueProperties.aaIndex[sq
- .charAt(j - 1)];
+ // sqnum[j] = ResidueProperties.aaIndex[sq.charAt(j - 1)];
+ sqnum[j] = sm.getMatrixIndex(sq.charAt(j - 1));
+ if (sqnum[j] == -1)
+ {
+ sqnum[j] = gapIndex;
+ }
}
seqNums.setElementAt(sqnum, i);
{
ResidueCount values = countResidues(column);
- // TODO is threshold a percentage or count value?
+ /*
+ * percentage count at or below which we ignore residues
+ */
int thresh = (threshold * height) / 100;
/*
* or not conserved (-1)
* Using TreeMap means properties are displayed in alphabetical order
*/
- Map<String, Integer> resultHash = new TreeMap<String, Integer>();
+ SortedMap<String, Integer> resultHash = new TreeMap<String, Integer>();
SymbolCounts symbolCounts = values.getSymbolCounts();
char[] symbols = symbolCounts.symbols;
int[] counts = symbolCounts.values;
*
* @return Conservation sequence
*/
- public Sequence getConsSequence()
+ public SequenceI getConsSequence()
{
return consSequence;
}
/**
* DOCUMENT ME!
+ *
+ * @param sm
*/
- private void percentIdentity2()
+ private void percentIdentity(ScoreMatrix sm)
{
seqNums = new Vector<int[]>();
- // calcSeqNum(s);
int i = 0, iSize = sequences.length;
// Do we need to calculate this again?
for (i = 0; i < iSize; i++)
{
- calcSeqNum(i);
+ calcSeqNum(i, sm);
}
+ int gapIndex = sm.getGapIndex();
+
if ((cons2 == null) || seqNumsChanged)
{
+ // FIXME remove magic number 24 without changing calc
+ // sm.getSize() returns 25 so doesn't quite do it...
cons2 = new int[maxLength][24];
- // Initialize the array
- for (int j = 0; j < 24; j++)
- {
- for (i = 0; i < maxLength; i++)
- {
- cons2[i][j] = 0;
- }
- }
-
- int[] sqnum;
int j = 0;
while (j < sequences.length)
{
- sqnum = seqNums.elementAt(j);
+ int[] sqnum = seqNums.elementAt(j);
for (i = 1; i < sqnum.length; i++)
{
for (i = sqnum.length - 1; i < maxLength; i++)
{
- cons2[i][23]++; // gap count
+ cons2[i][gapIndex]++; // gap count
}
-
j++;
}
-
- // unnecessary ?
-
- /*
- * for (int i=start; i <= end; i++) { int max = -1000; int maxi = -1; int
- * maxj = -1;
- *
- * for (int j=0;j<24;j++) { if (cons2[i][j] > max) { max = cons2[i][j];
- * maxi = i; maxj = j; } } }
- */
}
}
{
quality = new Vector<Double>();
- double max = -10000;
- int[][] BLOSUM62 = ResidueProperties.getBLOSUM62();
+ double max = -Double.MAX_VALUE;
+ ScoreMatrix blosum62 = ScoreModels.getInstance().getBlosum62();
+ float[][] blosumScores = blosum62.getMatrix();
// Loop over columns // JBPNote Profiling info
// long ts = System.currentTimeMillis();
// long te = System.currentTimeMillis();
- percentIdentity2();
+ percentIdentity(blosum62);
int size = seqNums.size();
int[] lengths = new int[size];
lengths[l] = seqNums.elementAt(l).length - 1;
}
+ // todo ? remove '*' (unused?) from score matrix and
+ // use getSize() here instead of getSize() - 1 ??
+ final int symbolCount = blosum62.getSize() - 1; // 24;
+ int gapIndex = blosum62.getGapIndex();
+
for (j = startRes; j <= endRes; j++)
{
bigtot = 0;
// First Xr = depends on column only
- x = new double[24];
+ x = new double[symbolCount];
- for (ii = 0; ii < 24; ii++)
+ for (ii = 0; ii < symbolCount; ii++)
{
x[ii] = 0;
- for (i2 = 0; i2 < 24; i2++)
+ for (i2 = 0; i2 < symbolCount; i2++)
{
- x[ii] += (((double) cons2[j][i2] * BLOSUM62[ii][i2]) + 4);
+ x[ii] += (((double) cons2[j][i2] * blosumScores[ii][i2]) + 4);
}
x[ii] /= size;
for (k = 0; k < size; k++)
{
tot = 0;
- xx = new double[24];
- seqNum = (j < lengths[k]) ? seqNums.elementAt(k)[j + 1] : 23; // Sequence,
- // or gap
- // at the
- // end
+ xx = new double[symbolCount];
+ seqNum = (j < lengths[k]) ? seqNums.elementAt(k)[j + 1] : gapIndex;
+ // Sequence, or gap at the end
// This is a loop over r
- for (i = 0; i < 23; i++)
+ for (i = 0; i < symbolCount - 1; i++)
{
sr = 0;
- sr = (double) BLOSUM62[i][seqNum] + 4;
+ sr = (double) blosumScores[i][seqNum] + 4;
// Calculate X with another loop over residues
// System.out.println("Xi " + i + " " + x[i] + " " + sr);
// Need to normalize by gaps
}
- double newmax = -10000;
+ double newmax = -Double.MAX_VALUE;
for (j = startRes; j <= endRes; j++)
{
tmp = quality.elementAt(j).doubleValue();
- tmp = ((max - tmp) * (size - cons2[j][23])) / size;
+ // tmp = ((max - tmp) * (size - cons2[j][23])) / size;
+ tmp = ((max - tmp) * (size - cons2[j][gapIndex])) / size;
// System.out.println(tmp+ " " + j);
quality.setElementAt(new Double(tmp), j);
*
* @param name
* - name of conservation
- * @param threshold
- * - minimum number of conserved residues needed to indicate
- * conservation (typically 3)
* @param seqs
* @param start
* first column in calculation window
* @return Conservation object ready for use in visualization
*/
public static Conservation calculateConservation(String name,
- int threshold, List<SequenceI> seqs, int start, int end,
- boolean positiveOnly, int maxPercentGaps, boolean calcQuality)
+ List<SequenceI> seqs, int start, int end, boolean positiveOnly,
+ int maxPercentGaps, boolean calcQuality)
{
- Conservation cons = new Conservation(name, threshold, seqs, start, end);
+ Conservation cons = new Conservation(name, seqs, start, end);
cons.calculate();
cons.verdict(positiveOnly, maxPercentGaps);
git://source.jalview.org / jalview.git / blobdiff