src/jalview/util/MappingUtils.java

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