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