src/jalview/util/MappingUtils.java

   1 package jalview.util;
   2
   3 import jalview.analysis.AlignmentSorter;
   4 import jalview.api.AlignViewportI;
   5 import jalview.commands.CommandI;
   6 import jalview.commands.EditCommand;
   7 import jalview.commands.EditCommand.Action;
   8 import jalview.commands.EditCommand.Edit;
   9 import jalview.commands.OrderCommand;
  10 import jalview.datamodel.AlignedCodonFrame;
  11 import jalview.datamodel.AlignmentI;
  12 import jalview.datamodel.AlignmentOrder;
  13 import jalview.datamodel.ColumnSelection;
  14 import jalview.datamodel.SearchResults;
  15 import jalview.datamodel.SearchResults.Match;
  16 import jalview.datamodel.Sequence;
  17 import jalview.datamodel.SequenceGroup;
  18 import jalview.datamodel.SequenceI;
  19
  20 import java.util.ArrayList;
  21 import java.util.HashMap;
  22 import java.util.Iterator;
  23 import java.util.List;
  24 import java.util.Map;
  25 import java.util.Set;
  26
  27 /**
  28  * Helper methods for manipulations involving sequence mappings.
  29  *
  30  * @author gmcarstairs
  31  *
  32  */
  33 public final class MappingUtils
  34 {
  35
  36   /**
  37    * Helper method to map a CUT or PASTE command.
  38    *
  39    * @param edit
  40    *          the original command
  41    * @param undo
  42    *          if true, the command is to be undone
  43    * @param targetSeqs
  44    *          the mapped sequences to apply the mapped command to
  45    * @param result
  46    *          the mapped EditCommand to add to
  47    * @param mappings
  48    */
  49   protected static void mapCutOrPaste(Edit edit, boolean undo,
  50           List<SequenceI> targetSeqs, EditCommand result,
  51           Set<AlignedCodonFrame> mappings)
  52   {
  53     Action action = edit.getAction();
  54     if (undo)
  55     {
  56       action = action.getUndoAction();
  57     }
  58     // TODO write this
  59     System.err.println("MappingUtils.mapCutOrPaste not yet implemented");
  60   }
  61
  62   /**
  63    * Returns a new EditCommand representing the given command as mapped to the
  64    * given sequences. If there is no mapping, returns null.
  65    *
  66    * @param command
  67    * @param undo
  68    * @param mapTo
  69    * @param gapChar
  70    * @param mappings
  71    * @return
  72    */
  73   public static EditCommand mapEditCommand(EditCommand command,
  74           boolean undo, final AlignmentI mapTo, char gapChar,
  75           Set<AlignedCodonFrame> mappings)
  76   {
  77     /*
  78      * For now, only support mapping from protein edits to cDna
  79      */
  80     if (!mapTo.isNucleotide())
  81     {
  82       return null;
  83     }
  84
  85     /*
  86      * Cache a copy of the target sequences so we can mimic successive edits on
  87      * them. This lets us compute mappings for all edits in the set.
  88      */
  89     Map<SequenceI, SequenceI> targetCopies = new HashMap<SequenceI, SequenceI>();
  90     for (SequenceI seq : mapTo.getSequences())
  91     {
  92       SequenceI ds = seq.getDatasetSequence();
  93       if (ds != null)
  94       {
  95         final SequenceI copy = new Sequence(seq);
  96         copy.setDatasetSequence(ds);
  97         targetCopies.put(ds, copy);
  98       }
  99     }
 100
 101     /*
 102      * Compute 'source' sequences as they were before applying edits:
 103      */
 104     Map<SequenceI, SequenceI> originalSequences = command.priorState(undo);
 105
 106     EditCommand result = new EditCommand();
 107     Iterator<Edit> edits = command.getEditIterator(!undo);
 108     while (edits.hasNext())
 109     {
 110       Edit edit = edits.next();
 111       if (edit.getAction() == Action.CUT
 112               || edit.getAction() == Action.PASTE)
 113       {
 114         mapCutOrPaste(edit, undo, mapTo.getSequences(), result, mappings);
 115       }
 116       else if (edit.getAction() == Action.INSERT_GAP
 117               || edit.getAction() == Action.DELETE_GAP)
 118       {
 119         mapInsertOrDelete(edit, undo, originalSequences,
 120                 mapTo.getSequences(), targetCopies, gapChar, result,
 121                 mappings);
 122       }
 123     }
 124     return result.getSize() > 0 ? result : null;
 125   }
 126
 127   /**
 128    * Helper method to map an edit command to insert or delete gaps.
 129    *
 130    * @param edit
 131    *          the original command
 132    * @param undo
 133    *          if true, the action is to undo the command
 134    * @param originalSequences
 135    *          the sequences the command acted on
 136    * @param targetSeqs
 137    * @param targetCopies
 138    * @param gapChar
 139    * @param result
 140    *          the new EditCommand to add mapped commands to
 141    * @param mappings
 142    */
 143   protected static void mapInsertOrDelete(Edit edit, boolean undo,
 144           Map<SequenceI, SequenceI> originalSequences,
 145           final List<SequenceI> targetSeqs,
 146           Map<SequenceI, SequenceI> targetCopies, char gapChar,
 147           EditCommand result, Set<AlignedCodonFrame> mappings)
 148   {
 149     Action action = edit.getAction();
 150
 151     /*
 152      * Invert sense of action if an Undo.
 153      */
 154     if (undo)
 155     {
 156       action = action.getUndoAction();
 157     }
 158     final int count = edit.getNumber();
 159     final int editPos = edit.getPosition();
 160     for (SequenceI seq : edit.getSequences())
 161     {
 162       /*
 163        * Get residue position at (or to right of) edit location. Note we use our
 164        * 'copy' of the sequence before editing for this.
 165        */
 166       SequenceI ds = seq.getDatasetSequence();
 167       if (ds == null)
 168       {
 169         continue;
 170       }
 171       final SequenceI actedOn = originalSequences.get(ds);
 172       final int seqpos = actedOn.findPosition(editPos);
 173
 174       /*
 175        * Determine all mappings from this position to mapped sequences.
 176        */
 177       SearchResults sr = buildSearchResults(seq, seqpos, mappings);
 178
 179       if (!sr.isEmpty())
 180       {
 181         for (SequenceI targetSeq : targetSeqs)
 182         {
 183           ds = targetSeq.getDatasetSequence();
 184           if (ds == null)
 185           {
 186             continue;
 187           }
 188           SequenceI copyTarget = targetCopies.get(ds);
 189           final int[] match = sr.getResults(copyTarget, 0,
 190                   copyTarget.getLength());
 191           if (match != null)
 192           {
 193             final int ratio = 3; // TODO: compute this - how?
 194             final int mappedCount = count * ratio;
 195
 196             /*
 197              * Shift Delete start position left, as it acts on positions to its
 198              * right.
 199              */
 200             int mappedEditPos = action == Action.DELETE_GAP ? match[0]
 201                     - mappedCount : match[0];
 202             Edit e = result.new Edit(action, new SequenceI[]
 203             { targetSeq }, mappedEditPos, mappedCount, gapChar);
 204             result.addEdit(e);
 205
 206             /*
 207              * and 'apply' the edit to our copy of its target sequence
 208              */
 209             if (action == Action.INSERT_GAP)
 210             {
 211               copyTarget.setSequence(new String(StringUtils.insertCharAt(
 212                       copyTarget.getSequence(), mappedEditPos, mappedCount,
 213                       gapChar)));
 214             }
 215             else if (action == Action.DELETE_GAP)
 216             {
 217               copyTarget.setSequence(new String(StringUtils.deleteChars(
 218                       copyTarget.getSequence(), mappedEditPos,
 219                       mappedEditPos + mappedCount)));
 220             }
 221           }
 222         }
 223       }
 224       /*
 225        * and 'apply' the edit to our copy of its source sequence
 226        */
 227       if (action == Action.INSERT_GAP)
 228       {
 229         actedOn.setSequence(new String(StringUtils.insertCharAt(
 230                 actedOn.getSequence(), editPos, count, gapChar)));
 231       }
 232       else if (action == Action.DELETE_GAP)
 233       {
 234         actedOn.setSequence(new String(StringUtils.deleteChars(
 235                 actedOn.getSequence(), editPos, editPos + count)));
 236       }
 237     }
 238   }
 239
 240   /**
 241    * Returns a SearchResults object describing the mapped region corresponding
 242    * to the specified sequence position.
 243    *
 244    * @param seq
 245    * @param index
 246    * @param seqmappings
 247    * @return
 248    */
 249   public static SearchResults buildSearchResults(SequenceI seq, int index,
 250           Set<AlignedCodonFrame> seqmappings)
 251   {
 252     SearchResults results;
 253     results = new SearchResults();
 254     if (index >= seq.getStart() && index <= seq.getEnd())
 255     {
 256       for (AlignedCodonFrame acf : seqmappings)
 257       {
 258         acf.markMappedRegion(seq, index, results);
 259       }
 260     }
 261     return results;
 262   }
 263
 264   /**
 265    * Returns a (possibly empty) SequenceGroup containing any sequences in the
 266    * mapped viewport corresponding to the given group in the source viewport.
 267    *
 268    * @param sg
 269    * @param mapFrom
 270    * @param mapTo
 271    * @return
 272    */
 273   public static SequenceGroup mapSequenceGroup(SequenceGroup sg,
 274           AlignViewportI mapFrom, AlignViewportI mapTo)
 275   {
 276     /*
 277      * Note the SequenceGroup holds aligned sequences, the mappings hold dataset
 278      * sequences.
 279      */
 280     boolean targetIsNucleotide = mapTo.isNucleotide();
 281     AlignViewportI protein = targetIsNucleotide ? mapFrom : mapTo;
 282     Set<AlignedCodonFrame> codonFrames = protein.getAlignment()
 283             .getCodonFrames();
 284
 285     /*
 286      * Copy group name, name colours, but not sequences or sequence colour
 287      * scheme
 288      */
 289     SequenceGroup mappedGroup = new SequenceGroup(sg);
 290     mappedGroup.cs = mapTo.getGlobalColourScheme();
 291     mappedGroup.clear();
 292     // TODO set width of mapped group
 293
 294     for (SequenceI selected : sg.getSequences())
 295     {
 296       for (AlignedCodonFrame acf : codonFrames)
 297       {
 298         SequenceI mappedSequence = targetIsNucleotide ? acf
 299                 .getDnaForAaSeq(selected) : acf.getAaForDnaSeq(selected);
 300         if (mappedSequence != null)
 301         {
 302           for (SequenceI seq : mapTo.getAlignment().getSequences())
 303           {
 304             if (seq.getDatasetSequence() == mappedSequence)
 305             {
 306               mappedGroup.addSequence(seq, false);
 307               break;
 308             }
 309           }
 310         }
 311       }
 312     }
 313     return mappedGroup;
 314   }
 315
 316   /**
 317    * Returns an OrderCommand equivalent to the given one, but acting on mapped
 318    * sequences as described by the mappings, or null if no mapping can be made.
 319    *
 320    * @param command
 321    *          the original order command
 322    * @param undo
 323    *          if true, the action is to undo the sort
 324    * @param mapTo
 325    *          the alignment we are mapping to
 326    * @param mappings
 327    *          the mappings available
 328    * @return
 329    */
 330   public static CommandI mapOrderCommand(OrderCommand command,
 331           boolean undo, AlignmentI mapTo, Set<AlignedCodonFrame> mappings)
 332   {
 333     SequenceI[] sortOrder = command.getSequenceOrder(undo);
 334     List<SequenceI> mappedOrder = new ArrayList<SequenceI>();
 335     int j = 0;
 336     for (SequenceI seq : sortOrder)
 337     {
 338       for (AlignedCodonFrame acf : mappings)
 339       {
 340         /*
 341          * Try protein-to-Dna, failing that try dna-to-protein
 342          */
 343         SequenceI mappedSeq = acf.getDnaForAaSeq(seq);
 344         if (mappedSeq == null)
 345         {
 346           mappedSeq = acf.getAaForDnaSeq(seq);
 347         }
 348         if (mappedSeq != null)
 349         {
 350           for (SequenceI seq2 : mapTo.getSequences())
 351           {
 352             if (seq2.getDatasetSequence() == mappedSeq)
 353             {
 354               mappedOrder.add(seq2);
 355               j++;
 356               break;
 357             }
 358           }
 359         }
 360       }
 361     }
 362
 363     /*
 364      * Return null if no mappings made.
 365      */
 366     if (j == 0)
 367     {
 368       return null;
 369     }
 370
 371     /*
 372      * Add any unmapped sequences on the end of the sort in their original
 373      * ordering.
 374      */
 375     if (j < mapTo.getHeight())
 376     {
 377       for (SequenceI seq : mapTo.getSequences())
 378       {
 379         if (!mappedOrder.contains(seq))
 380         {
 381           mappedOrder.add(seq);
 382         }
 383       }
 384     }
 385
 386     /*
 387      * Have to sort the sequences before constructing the OrderCommand - which
 388      * then resorts them?!?
 389      */
 390     final SequenceI[] mappedOrderArray = mappedOrder
 391             .toArray(new SequenceI[mappedOrder.size()]);
 392     SequenceI[] oldOrder = mapTo.getSequencesArray();
 393     AlignmentSorter.sortBy(mapTo, new AlignmentOrder(mappedOrderArray));
 394     final OrderCommand result = new OrderCommand(command.getDescription(),
 395             oldOrder, mapTo);
 396     return result;
 397   }
 398
 399   /**
 400    * Returns a ColumnSelection in the 'mapTo' view which corresponds to the
 401    * given selection in the 'mapFrom' view. We assume one is nucleotide, the
 402    * other is protein (and holds the mappings from codons to protein residues).
 403    *
 404    * @param colsel
 405    * @param mapFrom
 406    * @param mapTo
 407    * @return
 408    */
 409   public static ColumnSelection mapColumnSelection(ColumnSelection colsel,
 410           AlignViewportI mapFrom, AlignViewportI mapTo)
 411   {
 412     boolean targetIsNucleotide = mapTo.isNucleotide();
 413     AlignViewportI protein = targetIsNucleotide ? mapFrom : mapTo;
 414     Set<AlignedCodonFrame> codonFrames = protein.getAlignment()
 415             .getCodonFrames();
 416     ColumnSelection mappedColumns = new ColumnSelection();
 417     char fromGapChar = mapFrom.getAlignment().getGapCharacter();
 418
 419     // FIXME allow for hidden columns
 420
 421     /*
 422      * For each mapped column, find the range of columns that residues in that
 423      * column map to.
 424      */
 425     for (Object obj : colsel.getSelected())
 426     {
 427       int col = ((Integer) obj).intValue();
 428       int mappedToMin = Integer.MAX_VALUE;
 429       int mappedToMax = Integer.MIN_VALUE;
 430
 431       /*
 432        * For each sequence in the 'from' alignment
 433        */
 434       for (SequenceI fromSeq : mapFrom.getAlignment().getSequences())
 435       {
 436         /*
 437          * Ignore gaps (unmapped anyway)
 438          */
 439         if (fromSeq.getCharAt(col) == fromGapChar)
 440         {
 441           continue;
 442         }
 443
 444         /*
 445          * Get the residue position and find the mapped position.
 446          */
 447         int residuePos = fromSeq.findPosition(col);
 448         SearchResults sr = buildSearchResults(fromSeq, residuePos,
 449                 codonFrames);
 450         for (Match m : sr.getResults())
 451         {
 452           int mappedStartResidue = m.getStart();
 453           int mappedEndResidue = m.getEnd();
 454           SequenceI mappedSeq = m.getSequence();
 455
 456           /*
 457            * Locate the aligned sequence whose dataset is mappedSeq. TODO a
 458            * datamodel that can do this efficiently.
 459            */
 460           for (SequenceI toSeq : mapTo.getAlignment().getSequences())
 461           {
 462             if (toSeq.getDatasetSequence() == mappedSeq)
 463             {
 464               int mappedStartCol = toSeq.findIndex(mappedStartResidue);
 465               int mappedEndCol = toSeq.findIndex(mappedEndResidue);
 466               mappedToMin = Math.min(mappedToMin, mappedStartCol);
 467               mappedToMax = Math.max(mappedToMax, mappedEndCol);
 468               // System.out.println(fromSeq.getName() + " mapped to cols "
 469               // + mappedStartCol + ":" + mappedEndCol);
 470               break;
 471               // note: remove break if we ever want to map one to many sequences
 472             }
 473           }
 474         }
 475       }
 476       /*
 477        * Add the range of mapped columns to the mapped selection (converting
 478        * base 1 to base 0). Note that this may include intron-only regions which
 479        * lie between the start and end ranges of the selection.
 480        */
 481       for (int i = mappedToMin; i <= mappedToMax; i++)
 482       {
 483         mappedColumns.addElement(i - 1);
 484       }
 485     }
 486     return mappedColumns;
 487   }
 488
 489   /**
 490    * Returns the mapped codon for a given aligned sequence column position (base
 491    * 0).
 492    *
 493    * @param seq
 494    *          an aligned peptide sequence
 495    * @param col
 496    *          an aligned column position (base 0)
 497    * @param mappings
 498    *          a set of codon mappings
 499    * @return the bases of the mapped codon in the cDNA dataset sequence, or null
 500    *         if not found
 501    */
 502   public static char[] findCodonFor(SequenceI seq, int col,
 503           Set<AlignedCodonFrame> mappings)
 504   {
 505     int dsPos = seq.findPosition(col);
 506     for (AlignedCodonFrame mapping : mappings)
 507     {
 508       if (mapping.involvesSequence(seq))
 509       {
 510         return mapping.getMappedCodon(seq.getDatasetSequence(), dsPos);
 511       }
 512     }
 513     return null;
 514   }
 515
 516   /**
 517    * Converts a series of [start, end] ranges into an array of individual
 518    * positions.
 519    *
 520    * @param ranges
 521    * @return
 522    */
 523   public static int[] flattenRanges(int[] ranges)
 524   {
 525     /*
 526      * Count how many positions altogether
 527      */
 528     int count = 0;
 529     for (int i = 0; i < ranges.length - 1; i += 2)
 530     {
 531       count += ranges[i + 1] - ranges[i] + 1;
 532     }
 533
 534     int[] result = new int[count];
 535     int k = 0;
 536     for (int i = 0; i < ranges.length - 1; i += 2)
 537     {
 538       for (int j = ranges[i]; j <= ranges[i + 1]; j++)
 539       {
 540         result[k++] = j;
 541       }
 542     }
 543     return result;
 544   }
 545 }