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   private static final String NAME = "Sequence Feature Similarity";
  45
  46   private String description;
  47
  48   FeatureRenderer fr;
  49
  50   /**
  51    * Constructor
  52    */
  53   public FeatureDistanceModel()
  54   {
  55   }
  56
  57   @Override
  58   public boolean configureFromAlignmentView(AlignmentViewPanel view)
  59
  60   {
  61     fr = view.cloneFeatureRenderer();
  62     return true;
  63   }
  64
  65   /**
  66    * Calculates a distance measure [i][j] between each pair of sequences as the
  67    * average number of features they have but do not share. That is, find the
  68    * features each sequence pair has at each column, ignore feature types they
  69    * have in common, and count the rest. The totals are normalised by the number
  70    * of columns processed.
  71    * <p>
  72    * The parameters argument provides settings for treatment of gap-residue
  73    * aligned positions, and whether the score is over the longer or shorter of
  74    * each pair of sequences
  75    *
  76    * @param seqData
  77    * @param params
  78    */
  79   @Override
  80   public MatrixI findDistances(AlignmentView seqData,
  81           SimilarityParamsI params)
  82   {
  83     SeqCigar[] seqs = seqData.getSequences();
  84     int noseqs = seqs.length;
  85     int cpwidth = 0;// = seqData.getWidth();
  86     double[][] distances = new double[noseqs][noseqs];
  87     List<String> dft = null;
  88     if (fr != null)
  89     {
  90       dft = fr.getDisplayedFeatureTypes();
  91     }
  92     if (dft == null || dft.isEmpty())
  93     {
  94       return new Matrix(distances);
  95     }
  96
  97     // need to get real position for view position
  98     int[] viscont = seqData.getVisibleContigs();
  99
 100     /*
 101      * scan each column, compute and add to each distance[i, j]
 102      * the number of feature types that seqi and seqj do not share
 103      */
 104     for (int vc = 0; vc < viscont.length; vc += 2)
 105     {
 106       for (int cpos = viscont[vc]; cpos <= viscont[vc + 1]; cpos++)
 107       {
 108         cpwidth++;
 109
 110         /*
 111          * first record feature types in this column for each sequence
 112          */
 113         Map<SeqCigar, Set<String>> sfap = findFeatureTypesAtColumn(
 114                 seqs, cpos);
 115
 116         /*
 117          * count feature types on either i'th or j'th sequence but not both
 118          * and add this 'distance' measure to the total for [i, j] for j > i
 119          */
 120         for (int i = 0; i < (noseqs - 1); i++)
 121         {
 122           for (int j = i + 1; j < noseqs; j++)
 123           {
 124             SeqCigar sc1 = seqs[i];
 125             SeqCigar sc2 = seqs[j];
 126             Set<String> set1 = sfap.get(sc1);
 127             Set<String> set2 = sfap.get(sc2);
 128             boolean gap1 = set1 == null;
 129             boolean gap2 = set2 == null;
 130
 131             /*
 132              * gap-gap always scores zero
 133              * residue-residue is always scored
 134              * include gap-residue score if params say to do so
 135              */
 136             if ((!gap1 && !gap2) || params.includeGaps())
 137             {
 138               int seqDistance = SetUtils.countDisjunction(set1, set2);
 139               distances[i][j] += seqDistance;
 140             }
 141           }
 142         }
 143       }
 144     }
 145
 146     /*
 147      * normalise the distance scores (summed over columns) by the
 148      * number of visible columns used in the calculation
 149      * and fill in the bottom half of the matrix
 150      */
 151     // TODO JAL-2424 cpwidth may be out by 1 - affects scores but not tree shape
 152     for (int i = 0; i < noseqs; i++)
 153     {
 154       for (int j = i + 1; j < noseqs; j++)
 155       {
 156         distances[i][j] /= cpwidth;
 157         distances[j][i] = distances[i][j];
 158       }
 159     }
 160     return new Matrix(distances);
 161   }
 162
 163   /**
 164    * Builds and returns a map containing a (possibly empty) list (one per
 165    * SeqCigar) of visible feature types at the given column position. The map
 166    * has no entry for sequences which are gapped at the column position.
 167    *
 168    * @param seqs
 169    * @param columnPosition
 170    * @return
 171    */
 172   protected Map<SeqCigar, Set<String>> findFeatureTypesAtColumn(
 173           SeqCigar[] seqs, int columnPosition)
 174   {
 175     Map<SeqCigar, Set<String>> sfap = new HashMap<SeqCigar, Set<String>>();
 176     for (SeqCigar seq : seqs)
 177     {
 178       int spos = seq.findPosition(columnPosition);
 179       if (spos != -1)
 180       {
 181         Set<String> types = new HashSet<String>();
 182         List<SequenceFeature> sfs = fr.findFeaturesAtRes(seq.getRefSeq(),
 183                 spos);
 184         for (SequenceFeature sf : sfs)
 185         {
 186           types.add(sf.getType());
 187         }
 188         sfap.put(seq, types);
 189       }
 190     }
 191     return sfap;
 192   }
 193
 194   @Override
 195   public String getName()
 196   {
 197     return NAME;
 198   }
 199
 200   @Override
 201   public String getDescription()
 202   {
 203     return description;
 204   }
 205
 206   @Override
 207   public boolean isDNA()
 208   {
 209     return true;
 210   }
 211
 212   @Override
 213   public boolean isProtein()
 214   {
 215     return true;
 216   }
 217
 218   @Override
 219   public String toString()
 220   {
 221     return "Score between sequences based on hamming distance between binary vectors marking features displayed at each column";
 222   }
 223 }