src/jalview/analysis/scoremodels/FeatureDistanceModel.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
   3  * Copyright (C) $$Year-Rel$$ The Jalview Authors
   4  *
   5  * This file is part of Jalview.
   6  *
   7  * Jalview is free software: you can redistribute it and/or
   8  * modify it under the terms of the GNU General Public License
   9  * as published by the Free Software Foundation, either version 3
  10  * of the License, or (at your option) any later version.
  11  *
  12  * Jalview is distributed in the hope that it will be useful, but
  13  * WITHOUT ANY WARRANTY; without even the implied warranty
  14  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  15  * PURPOSE.  See the GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
  19  * The Jalview Authors are detailed in the 'AUTHORS' file.
  20  */
  21 package jalview.analysis.scoremodels;
  22
  23 import jalview.api.AlignmentViewPanel;
  24 import jalview.api.FeatureRenderer;
  25 import jalview.api.analysis.DistanceScoreModelI;
  26 import jalview.api.analysis.SimilarityParamsI;
  27 import jalview.api.analysis.ViewBasedAnalysisI;
  28 import jalview.datamodel.AlignmentView;
  29 import jalview.datamodel.SeqCigar;
  30 import jalview.datamodel.SequenceFeature;
  31 import jalview.math.Matrix;
  32 import jalview.math.MatrixI;
  33 import jalview.util.SetUtils;
  34
  35 import java.util.HashMap;
  36 import java.util.HashSet;
  37 import java.util.List;
  38 import java.util.Map;
  39 import java.util.Set;
  40
  41 public class FeatureDistanceModel implements DistanceScoreModelI,
  42         ViewBasedAnalysisI
  43 {
  44   FeatureRenderer fr;
  45
  46   @Override
  47   public boolean configureFromAlignmentView(AlignmentViewPanel view)
  48
  49   {
  50     fr = view.cloneFeatureRenderer();
  51     return true;
  52   }
  53
  54   /**
  55    * Calculates a distance measure [i][j] between each pair of sequences as the
  56    * average number of features they have but do not share. That is, find the
  57    * features each sequence pair has at each column, ignore feature types they
  58    * have in common, and count the rest. The totals are normalised by the number
  59    * of columns processed.
  60    * <p>
  61    * The parameters argument provides settings for treatment of gap-residue
  62    * aligned positions, and whether the score is over the longer or shorter of
  63    * each pair of sequences
  64    *
  65    * @param seqData
  66    * @param params
  67    */
  68   @Override
  69   public MatrixI findDistances(AlignmentView seqData,
  70           SimilarityParamsI params)
  71   {
  72     SeqCigar[] seqs = seqData.getSequences();
  73     int noseqs = seqs.length;
  74     int cpwidth = 0;// = seqData.getWidth();
  75     double[][] distances = new double[noseqs][noseqs];
  76     List<String> dft = null;
  77     if (fr != null)
  78     {
  79       dft = fr.getDisplayedFeatureTypes();
  80     }
  81     if (dft == null || dft.isEmpty())
  82     {
  83       return new Matrix(distances);
  84     }
  85
  86     // need to get real position for view position
  87     int[] viscont = seqData.getVisibleContigs();
  88
  89     /*
  90      * scan each column, compute and add to each distance[i, j]
  91      * the number of feature types that seqi and seqj do not share
  92      */
  93     for (int vc = 0; vc < viscont.length; vc += 2)
  94     {
  95       for (int cpos = viscont[vc]; cpos <= viscont[vc + 1]; cpos++)
  96       {
  97         cpwidth++;
  98
  99         /*
 100          * first record feature types in this column for each sequence
 101          */
 102         Map<SeqCigar, Set<String>> sfap = findFeatureTypesAtColumn(
 103                 seqs, cpos);
 104
 105         /*
 106          * count feature types on either i'th or j'th sequence but not both
 107          * and add this 'distance' measure to the total for [i, j] for j > i
 108          */
 109         for (int i = 0; i < (noseqs - 1); i++)
 110         {
 111           for (int j = i + 1; j < noseqs; j++)
 112           {
 113             SeqCigar sc1 = seqs[i];
 114             SeqCigar sc2 = seqs[j];
 115             Set<String> set1 = sfap.get(sc1);
 116             Set<String> set2 = sfap.get(sc2);
 117             boolean gap1 = set1 == null;
 118             boolean gap2 = set2 == null;
 119
 120             /*
 121              * gap-gap always scores zero
 122              * residue-residue is always scored
 123              * include gap-residue score if params say to do so
 124              */
 125             if ((!gap1 && !gap2) || params.includeGaps())
 126             {
 127               int seqDistance = SetUtils.countDisjunction(set1, set2);
 128               distances[i][j] += seqDistance;
 129             }
 130           }
 131         }
 132       }
 133     }
 134
 135     /*
 136      * normalise the distance scores (summed over columns) by the
 137      * number of visible columns used in the calculation
 138      * and fill in the bottom half of the matrix
 139      */
 140     // TODO JAL-2424 cpwidth may be out by 1 - affects scores but not tree shape
 141     for (int i = 0; i < noseqs; i++)
 142     {
 143       for (int j = i + 1; j < noseqs; j++)
 144       {
 145         distances[i][j] /= cpwidth;
 146         distances[j][i] = distances[i][j];
 147       }
 148     }
 149     return new Matrix(distances);
 150   }
 151
 152   /**
 153    * Builds and returns a map containing a (possibly empty) list (one per
 154    * SeqCigar) of visible feature types at the given column position. The map
 155    * has no entry for sequences which are gapped at the column position.
 156    *
 157    * @param seqs
 158    * @param columnPosition
 159    * @return
 160    */
 161   protected Map<SeqCigar, Set<String>> findFeatureTypesAtColumn(
 162           SeqCigar[] seqs, int columnPosition)
 163   {
 164     Map<SeqCigar, Set<String>> sfap = new HashMap<SeqCigar, Set<String>>();
 165     for (SeqCigar seq : seqs)
 166     {
 167       int spos = seq.findPosition(columnPosition);
 168       if (spos != -1)
 169       {
 170         Set<String> types = new HashSet<String>();
 171         List<SequenceFeature> sfs = fr.findFeaturesAtRes(seq.getRefSeq(),
 172                 spos);
 173         for (SequenceFeature sf : sfs)
 174         {
 175           types.add(sf.getType());
 176         }
 177         sfap.put(seq, types);
 178       }
 179     }
 180     return sfap;
 181   }
 182
 183   @Override
 184   public String getName()
 185   {
 186     return "Sequence Feature Similarity";
 187   }
 188
 189   @Override
 190   public boolean isDNA()
 191   {
 192     return true;
 193   }
 194
 195   @Override
 196   public boolean isProtein()
 197   {
 198     return true;
 199   }
 200
 201   @Override
 202   public String toString()
 203   {
 204     return "Score between sequences based on hamming distance between binary vectors marking features displayed at each column";
 205   }
 206 }