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.HiddenColumns;
  35 import jalview.datamodel.SearchResultMatchI;
  36 import jalview.datamodel.SearchResults;
  37 import jalview.datamodel.SearchResultsI;
  38 import jalview.datamodel.Sequence;
  39 import jalview.datamodel.SequenceGroup;
  40 import jalview.datamodel.SequenceI;
  41
  42 import java.util.ArrayList;
  43 import java.util.Arrays;
  44 import java.util.HashMap;
  45 import java.util.Iterator;
  46 import java.util.List;
  47 import java.util.Map;
  48
  49 /**
  50  * Helper methods for manipulations involving sequence mappings.
  51  *
  52  * @author gmcarstairs
  53  *
  54  */
  55 public final class MappingUtils
  56 {
  57
  58   /**
  59    * Helper method to map a CUT or PASTE command.
  60    *
  61    * @param edit
  62    *          the original command
  63    * @param undo
  64    *          if true, the command is to be undone
  65    * @param targetSeqs
  66    *          the mapped sequences to apply the mapped command to
  67    * @param result
  68    *          the mapped EditCommand to add to
  69    * @param mappings
  70    */
  71   protected static void mapCutOrPaste(Edit edit, boolean undo,
  72           List<SequenceI> targetSeqs, EditCommand result,
  73           List<AlignedCodonFrame> mappings)
  74   {
  75     Action action = edit.getAction();
  76     if (undo)
  77     {
  78       action = action.getUndoAction();
  79     }
  80     // TODO write this
  81     System.err.println("MappingUtils.mapCutOrPaste not yet implemented");
  82   }
  83
  84   /**
  85    * Returns a new EditCommand representing the given command as mapped to the
  86    * given sequences. If there is no mapping, returns null.
  87    *
  88    * @param command
  89    * @param undo
  90    * @param mapTo
  91    * @param gapChar
  92    * @param mappings
  93    * @return
  94    */
  95   public static EditCommand mapEditCommand(EditCommand command,
  96           boolean undo, final AlignmentI mapTo, char gapChar,
  97           List<AlignedCodonFrame> mappings)
  98   {
  99     /*
 100      * For now, only support mapping from protein edits to cDna
 101      */
 102     if (!mapTo.isNucleotide())
 103     {
 104       return null;
 105     }
 106
 107     /*
 108      * Cache a copy of the target sequences so we can mimic successive edits on
 109      * them. This lets us compute mappings for all edits in the set.
 110      */
 111     Map<SequenceI, SequenceI> targetCopies = new HashMap<>();
 112     for (SequenceI seq : mapTo.getSequences())
 113     {
 114       SequenceI ds = seq.getDatasetSequence();
 115       if (ds != null)
 116       {
 117         final SequenceI copy = new Sequence(seq);
 118         copy.setDatasetSequence(ds);
 119         targetCopies.put(ds, copy);
 120       }
 121     }
 122
 123     /*
 124      * Compute 'source' sequences as they were before applying edits:
 125      */
 126     Map<SequenceI, SequenceI> originalSequences = command.priorState(undo);
 127
 128     EditCommand result = new EditCommand();
 129     Iterator<Edit> edits = command.getEditIterator(!undo);
 130     while (edits.hasNext())
 131     {
 132       Edit edit = edits.next();
 133       if (edit.getAction() == Action.CUT
 134               || edit.getAction() == Action.PASTE)
 135       {
 136         mapCutOrPaste(edit, undo, mapTo.getSequences(), result, mappings);
 137       }
 138       else if (edit.getAction() == Action.INSERT_GAP
 139               || edit.getAction() == Action.DELETE_GAP)
 140       {
 141         mapInsertOrDelete(edit, undo, originalSequences,
 142                 mapTo.getSequences(), targetCopies, gapChar, result,
 143                 mappings);
 144       }
 145     }
 146     return result.getSize() > 0 ? result : null;
 147   }
 148
 149   /**
 150    * Helper method to map an edit command to insert or delete gaps.
 151    *
 152    * @param edit
 153    *          the original command
 154    * @param undo
 155    *          if true, the action is to undo the command
 156    * @param originalSequences
 157    *          the sequences the command acted on
 158    * @param targetSeqs
 159    * @param targetCopies
 160    * @param gapChar
 161    * @param result
 162    *          the new EditCommand to add mapped commands to
 163    * @param mappings
 164    */
 165   protected static void mapInsertOrDelete(Edit edit, boolean undo,
 166           Map<SequenceI, SequenceI> originalSequences,
 167           final List<SequenceI> targetSeqs,
 168           Map<SequenceI, SequenceI> targetCopies, char gapChar,
 169           EditCommand result, List<AlignedCodonFrame> mappings)
 170   {
 171     Action action = edit.getAction();
 172
 173     /*
 174      * Invert sense of action if an Undo.
 175      */
 176     if (undo)
 177     {
 178       action = action.getUndoAction();
 179     }
 180     final int count = edit.getNumber();
 181     final int editPos = edit.getPosition();
 182     for (SequenceI seq : edit.getSequences())
 183     {
 184       /*
 185        * Get residue position at (or to right of) edit location. Note we use our
 186        * 'copy' of the sequence before editing for this.
 187        */
 188       SequenceI ds = seq.getDatasetSequence();
 189       if (ds == null)
 190       {
 191         continue;
 192       }
 193       final SequenceI actedOn = originalSequences.get(ds);
 194       final int seqpos = actedOn.findPosition(editPos);
 195
 196       /*
 197        * Determine all mappings from this position to mapped sequences.
 198        */
 199       SearchResultsI sr = buildSearchResults(seq, seqpos, mappings);
 200
 201       if (!sr.isEmpty())
 202       {
 203         for (SequenceI targetSeq : targetSeqs)
 204         {
 205           ds = targetSeq.getDatasetSequence();
 206           if (ds == null)
 207           {
 208             continue;
 209           }
 210           SequenceI copyTarget = targetCopies.get(ds);
 211           final int[] match = sr.getResults(copyTarget, 0,
 212                   copyTarget.getLength());
 213           if (match != null)
 214           {
 215             final int ratio = 3; // TODO: compute this - how?
 216             final int mappedCount = count * ratio;
 217
 218             /*
 219              * Shift Delete start position left, as it acts on positions to its
 220              * right.
 221              */
 222             int mappedEditPos = action == Action.DELETE_GAP
 223                     ? match[0] - mappedCount
 224                     : match[0];
 225             Edit e = result.new Edit(action, new SequenceI[] { targetSeq },
 226                     mappedEditPos, mappedCount, gapChar);
 227             result.addEdit(e);
 228
 229             /*
 230              * and 'apply' the edit to our copy of its target sequence
 231              */
 232             if (action == Action.INSERT_GAP)
 233             {
 234               copyTarget.setSequence(new String(
 235                       StringUtils.insertCharAt(copyTarget.getSequence(),
 236                               mappedEditPos, mappedCount, gapChar)));
 237             }
 238             else if (action == Action.DELETE_GAP)
 239             {
 240               copyTarget.setSequence(new String(
 241                       StringUtils.deleteChars(copyTarget.getSequence(),
 242                               mappedEditPos, mappedEditPos + mappedCount)));
 243             }
 244           }
 245         }
 246       }
 247       /*
 248        * and 'apply' the edit to our copy of its source sequence
 249        */
 250       if (action == Action.INSERT_GAP)
 251       {
 252         actedOn.setSequence(new String(StringUtils.insertCharAt(
 253                 actedOn.getSequence(), editPos, count, gapChar)));
 254       }
 255       else if (action == Action.DELETE_GAP)
 256       {
 257         actedOn.setSequence(new String(StringUtils.deleteChars(
 258                 actedOn.getSequence(), editPos, editPos + count)));
 259       }
 260     }
 261   }
 262
 263   /**
 264    * Returns a SearchResults object describing the mapped region corresponding
 265    * to the specified sequence position.
 266    *
 267    * @param seq
 268    * @param index
 269    * @param seqmappings
 270    * @return
 271    */
 272   public static SearchResultsI buildSearchResults(SequenceI seq, int index,
 273           List<AlignedCodonFrame> seqmappings)
 274   {
 275     SearchResultsI results = new SearchResults();
 276     addSearchResults(results, seq, index, seqmappings);
 277     return results;
 278   }
 279
 280   /**
 281    * Adds entries to a SearchResults object describing the mapped region
 282    * corresponding to the specified sequence position.
 283    *
 284    * @param results
 285    * @param seq
 286    * @param index
 287    * @param seqmappings
 288    */
 289   public static void addSearchResults(SearchResultsI results, SequenceI seq,
 290           int index, List<AlignedCodonFrame> seqmappings)
 291   {
 292     if (index >= seq.getStart() && index <= seq.getEnd())
 293     {
 294       for (AlignedCodonFrame acf : seqmappings)
 295       {
 296         acf.markMappedRegion(seq, index, results);
 297       }
 298     }
 299   }
 300
 301   /**
 302    * Returns a (possibly empty) SequenceGroup containing any sequences in the
 303    * mapped viewport corresponding to the given group in the source viewport.
 304    *
 305    * @param sg
 306    * @param mapFrom
 307    * @param mapTo
 308    * @return
 309    */
 310   public static SequenceGroup mapSequenceGroup(final SequenceGroup sg,
 311           final AlignViewportI mapFrom, final AlignViewportI mapTo)
 312   {
 313     /*
 314      * Note the SequenceGroup holds aligned sequences, the mappings hold dataset
 315      * sequences.
 316      */
 317     boolean targetIsNucleotide = mapTo.isNucleotide();
 318     AlignViewportI protein = targetIsNucleotide ? mapFrom : mapTo;
 319     List<AlignedCodonFrame> codonFrames = protein.getAlignment()
 320             .getCodonFrames();
 321     /*
 322      * Copy group name, colours etc, but not sequences or sequence colour scheme
 323      */
 324     SequenceGroup mappedGroup = new SequenceGroup(sg);
 325     mappedGroup.setColourScheme(mapTo.getGlobalColourScheme());
 326     mappedGroup.clear();
 327
 328     int minStartCol = -1;
 329     int maxEndCol = -1;
 330     final int selectionStartRes = sg.getStartRes();
 331     final int selectionEndRes = sg.getEndRes();
 332     for (SequenceI selected : sg.getSequences())
 333     {
 334       /*
 335        * Find the widest range of non-gapped positions in the selection range
 336        */
 337       int firstUngappedPos = selectionStartRes;
 338       while (firstUngappedPos <= selectionEndRes
 339               && Comparison.isGap(selected.getCharAt(firstUngappedPos)))
 340       {
 341         firstUngappedPos++;
 342       }
 343
 344       /*
 345        * If this sequence is only gaps in the selected range, skip it
 346        */
 347       if (firstUngappedPos > selectionEndRes)
 348       {
 349         continue;
 350       }
 351
 352       int lastUngappedPos = selectionEndRes;
 353       while (lastUngappedPos >= selectionStartRes
 354               && Comparison.isGap(selected.getCharAt(lastUngappedPos)))
 355       {
 356         lastUngappedPos--;
 357       }
 358
 359       /*
 360        * Find the selected start/end residue positions in sequence
 361        */
 362       int startResiduePos = selected.findPosition(firstUngappedPos);
 363       int endResiduePos = selected.findPosition(lastUngappedPos);
 364
 365       for (AlignedCodonFrame acf : codonFrames)
 366       {
 367         SequenceI mappedSequence = targetIsNucleotide
 368                 ? acf.getDnaForAaSeq(selected)
 369                 : acf.getAaForDnaSeq(selected);
 370         if (mappedSequence != null)
 371         {
 372           for (SequenceI seq : mapTo.getAlignment().getSequences())
 373           {
 374             int mappedStartResidue = 0;
 375             int mappedEndResidue = 0;
 376             if (seq.getDatasetSequence() == mappedSequence)
 377             {
 378               /*
 379                * Found a sequence mapping. Locate the start/end mapped residues.
 380                */
 381               List<AlignedCodonFrame> mapping = Arrays
 382                       .asList(new AlignedCodonFrame[]
 383                       { acf });
 384               SearchResultsI sr = buildSearchResults(selected,
 385                       startResiduePos, mapping);
 386               for (SearchResultMatchI m : sr.getResults())
 387               {
 388                 mappedStartResidue = m.getStart();
 389                 mappedEndResidue = m.getEnd();
 390               }
 391               sr = buildSearchResults(selected, endResiduePos, mapping);
 392               for (SearchResultMatchI m : sr.getResults())
 393               {
 394                 mappedStartResidue = Math.min(mappedStartResidue,
 395                         m.getStart());
 396                 mappedEndResidue = Math.max(mappedEndResidue, m.getEnd());
 397               }
 398
 399               /*
 400                * Find the mapped aligned columns, save the range. Note findIndex
 401                * returns a base 1 position, SequenceGroup uses base 0
 402                */
 403               int mappedStartCol = seq.findIndex(mappedStartResidue) - 1;
 404               minStartCol = minStartCol == -1 ? mappedStartCol
 405                       : Math.min(minStartCol, mappedStartCol);
 406               int mappedEndCol = seq.findIndex(mappedEndResidue) - 1;
 407               maxEndCol = maxEndCol == -1 ? mappedEndCol
 408                       : Math.max(maxEndCol, mappedEndCol);
 409               mappedGroup.addSequence(seq, false);
 410               break;
 411             }
 412           }
 413         }
 414       }
 415     }
 416     mappedGroup.setStartRes(minStartCol < 0 ? 0 : minStartCol);
 417     mappedGroup.setEndRes(maxEndCol < 0 ? 0 : maxEndCol);
 418     return mappedGroup;
 419   }
 420
 421   /**
 422    * Returns an OrderCommand equivalent to the given one, but acting on mapped
 423    * sequences as described by the mappings, or null if no mapping can be made.
 424    *
 425    * @param command
 426    *          the original order command
 427    * @param undo
 428    *          if true, the action is to undo the sort
 429    * @param mapTo
 430    *          the alignment we are mapping to
 431    * @param mappings
 432    *          the mappings available
 433    * @return
 434    */
 435   public static CommandI mapOrderCommand(OrderCommand command, boolean undo,
 436           AlignmentI mapTo, List<AlignedCodonFrame> mappings)
 437   {
 438     SequenceI[] sortOrder = command.getSequenceOrder(undo);
 439     List<SequenceI> mappedOrder = new ArrayList<>();
 440     int j = 0;
 441
 442     /*
 443      * Assumption: we are only interested in a cDNA/protein mapping; refactor in
 444      * future if we want to support sorting (c)dna as (c)dna or protein as
 445      * protein
 446      */
 447     boolean mappingToNucleotide = mapTo.isNucleotide();
 448     for (SequenceI seq : sortOrder)
 449     {
 450       for (AlignedCodonFrame acf : mappings)
 451       {
 452         SequenceI mappedSeq = mappingToNucleotide ? acf.getDnaForAaSeq(seq)
 453                 : acf.getAaForDnaSeq(seq);
 454         if (mappedSeq != null)
 455         {
 456           for (SequenceI seq2 : mapTo.getSequences())
 457           {
 458             if (seq2.getDatasetSequence() == mappedSeq)
 459             {
 460               mappedOrder.add(seq2);
 461               j++;
 462               break;
 463             }
 464           }
 465         }
 466       }
 467     }
 468
 469     /*
 470      * Return null if no mappings made.
 471      */
 472     if (j == 0)
 473     {
 474       return null;
 475     }
 476
 477     /*
 478      * Add any unmapped sequences on the end of the sort in their original
 479      * ordering.
 480      */
 481     if (j < mapTo.getHeight())
 482     {
 483       for (SequenceI seq : mapTo.getSequences())
 484       {
 485         if (!mappedOrder.contains(seq))
 486         {
 487           mappedOrder.add(seq);
 488         }
 489       }
 490     }
 491
 492     /*
 493      * Have to sort the sequences before constructing the OrderCommand - which
 494      * then resorts them?!?
 495      */
 496     final SequenceI[] mappedOrderArray = mappedOrder
 497             .toArray(new SequenceI[mappedOrder.size()]);
 498     SequenceI[] oldOrder = mapTo.getSequencesArray();
 499     AlignmentSorter.sortBy(mapTo, new AlignmentOrder(mappedOrderArray));
 500     final OrderCommand result = new OrderCommand(command.getDescription(),
 501             oldOrder, mapTo);
 502     return result;
 503   }
 504
 505   /**
 506    * Returns a ColumnSelection in the 'mapTo' view which corresponds to the
 507    * given selection in the 'mapFrom' view. We assume one is nucleotide, the
 508    * other is protein (and holds the mappings from codons to protein residues).
 509    *
 510    * @param colsel
 511    * @param mapFrom
 512    * @param mapTo
 513    * @return
 514    */
 515   public static void mapColumnSelection(ColumnSelection colsel,
 516           HiddenColumns hiddencols, AlignViewportI mapFrom,
 517           AlignViewportI mapTo, ColumnSelection newColSel,
 518           HiddenColumns newHidden)
 519   {
 520     boolean targetIsNucleotide = mapTo.isNucleotide();
 521     AlignViewportI protein = targetIsNucleotide ? mapFrom : mapTo;
 522     List<AlignedCodonFrame> codonFrames = protein.getAlignment()
 523             .getCodonFrames();
 524
 525     if (colsel == null)
 526     {
 527       return; // mappedColumns;
 528     }
 529
 530     char fromGapChar = mapFrom.getAlignment().getGapCharacter();
 531
 532     /*
 533      * For each mapped column, find the range of columns that residues in that
 534      * column map to.
 535      */
 536     List<SequenceI> fromSequences = mapFrom.getAlignment().getSequences();
 537     List<SequenceI> toSequences = mapTo.getAlignment().getSequences();
 538
 539     for (Integer sel : colsel.getSelected())
 540     {
 541       mapColumn(sel.intValue(), codonFrames, newColSel, fromSequences,
 542               toSequences, fromGapChar);
 543     }
 544
 545     for (int[] hidden : hiddencols.getHiddenColumnsCopy())
 546     {
 547       mapHiddenColumns(hidden, codonFrames, newHidden, fromSequences,
 548               toSequences, fromGapChar);
 549     }
 550     return; // mappedColumns;
 551   }
 552
 553   /**
 554    * Helper method that maps a [start, end] hidden column range to its mapped
 555    * equivalent
 556    *
 557    * @param hidden
 558    * @param mappings
 559    * @param mappedColumns
 560    * @param fromSequences
 561    * @param toSequences
 562    * @param fromGapChar
 563    */
 564   protected static void mapHiddenColumns(int[] hidden,
 565           List<AlignedCodonFrame> mappings, HiddenColumns mappedColumns,
 566           List<SequenceI> fromSequences, List<SequenceI> toSequences,
 567           char fromGapChar)
 568   {
 569     for (int col = hidden[0]; col <= hidden[1]; col++)
 570     {
 571       int[] mappedTo = findMappedColumns(col, mappings, fromSequences,
 572               toSequences, fromGapChar);
 573
 574       /*
 575        * Add the range of hidden columns to the mapped selection (converting
 576        * base 1 to base 0).
 577        */
 578       if (mappedTo != null)
 579       {
 580         mappedColumns.hideColumns(mappedTo[0] - 1, mappedTo[1] - 1);
 581       }
 582     }
 583   }
 584
 585   /**
 586    * Helper method to map one column selection
 587    *
 588    * @param col
 589    *          the column number (base 0)
 590    * @param mappings
 591    *          the sequence mappings
 592    * @param mappedColumns
 593    *          the mapped column selections to add to
 594    * @param fromSequences
 595    * @param toSequences
 596    * @param fromGapChar
 597    */
 598   protected static void mapColumn(int col, List<AlignedCodonFrame> mappings,
 599           ColumnSelection mappedColumns, List<SequenceI> fromSequences,
 600           List<SequenceI> toSequences, char fromGapChar)
 601   {
 602     int[] mappedTo = findMappedColumns(col, mappings, fromSequences,
 603             toSequences, fromGapChar);
 604
 605     /*
 606      * Add the range of mapped columns to the mapped selection (converting
 607      * base 1 to base 0). Note that this may include intron-only regions which
 608      * lie between the start and end ranges of the selection.
 609      */
 610     if (mappedTo != null)
 611     {
 612       for (int i = mappedTo[0]; i <= mappedTo[1]; i++)
 613       {
 614         mappedColumns.addElement(i - 1);
 615       }
 616     }
 617   }
 618
 619   /**
 620    * Helper method to find the range of columns mapped to from one column.
 621    * Returns the maximal range of columns mapped to from all sequences in the
 622    * source column, or null if no mappings were found.
 623    *
 624    * @param col
 625    * @param mappings
 626    * @param fromSequences
 627    * @param toSequences
 628    * @param fromGapChar
 629    * @return
 630    */
 631   protected static int[] findMappedColumns(int col,
 632           List<AlignedCodonFrame> mappings, List<SequenceI> fromSequences,
 633           List<SequenceI> toSequences, char fromGapChar)
 634   {
 635     int[] mappedTo = new int[] { Integer.MAX_VALUE, Integer.MIN_VALUE };
 636     boolean found = false;
 637
 638     /*
 639      * For each sequence in the 'from' alignment
 640      */
 641     for (SequenceI fromSeq : fromSequences)
 642     {
 643       /*
 644        * Ignore gaps (unmapped anyway)
 645        */
 646       if (fromSeq.getCharAt(col) == fromGapChar)
 647       {
 648         continue;
 649       }
 650
 651       /*
 652        * Get the residue position and find the mapped position.
 653        */
 654       int residuePos = fromSeq.findPosition(col);
 655       SearchResultsI sr = buildSearchResults(fromSeq, residuePos, mappings);
 656       for (SearchResultMatchI m : sr.getResults())
 657       {
 658         int mappedStartResidue = m.getStart();
 659         int mappedEndResidue = m.getEnd();
 660         SequenceI mappedSeq = m.getSequence();
 661
 662         /*
 663          * Locate the aligned sequence whose dataset is mappedSeq. TODO a
 664          * datamodel that can do this efficiently.
 665          */
 666         for (SequenceI toSeq : toSequences)
 667         {
 668           if (toSeq.getDatasetSequence() == mappedSeq)
 669           {
 670             int mappedStartCol = toSeq.findIndex(mappedStartResidue);
 671             int mappedEndCol = toSeq.findIndex(mappedEndResidue);
 672             mappedTo[0] = Math.min(mappedTo[0], mappedStartCol);
 673             mappedTo[1] = Math.max(mappedTo[1], mappedEndCol);
 674             found = true;
 675             break;
 676             // note: remove break if we ever want to map one to many sequences
 677           }
 678         }
 679       }
 680     }
 681     return found ? mappedTo : null;
 682   }
 683
 684   /**
 685    * Returns the mapped codon or codons for a given aligned sequence column
 686    * position (base 0).
 687    *
 688    * @param seq
 689    *          an aligned peptide sequence
 690    * @param col
 691    *          an aligned column position (base 0)
 692    * @param mappings
 693    *          a set of codon mappings
 694    * @return the bases of the mapped codon(s) in the cDNA dataset sequence(s),
 695    *         or an empty list if none found
 696    */
 697   public static List<char[]> findCodonsFor(SequenceI seq, int col,
 698           List<AlignedCodonFrame> mappings)
 699   {
 700     List<char[]> result = new ArrayList<>();
 701     int dsPos = seq.findPosition(col);
 702     for (AlignedCodonFrame mapping : mappings)
 703     {
 704       if (mapping.involvesSequence(seq))
 705       {
 706         List<char[]> codons = mapping
 707                 .getMappedCodons(seq.getDatasetSequence(), dsPos);
 708         if (codons != null)
 709         {
 710           result.addAll(codons);
 711         }
 712       }
 713     }
 714     return result;
 715   }
 716
 717   /**
 718    * Converts a series of [start, end] range pairs into an array of individual
 719    * positions. This also caters for 'reverse strand' (start > end) cases.
 720    *
 721    * @param ranges
 722    * @return
 723    */
 724   public static int[] flattenRanges(int[] ranges)
 725   {
 726     /*
 727      * Count how many positions altogether
 728      */
 729     int count = 0;
 730     for (int i = 0; i < ranges.length - 1; i += 2)
 731     {
 732       count += Math.abs(ranges[i + 1] - ranges[i]) + 1;
 733     }
 734
 735     int[] result = new int[count];
 736     int k = 0;
 737     for (int i = 0; i < ranges.length - 1; i += 2)
 738     {
 739       int from = ranges[i];
 740       final int to = ranges[i + 1];
 741       int step = from <= to ? 1 : -1;
 742       do
 743       {
 744         result[k++] = from;
 745         from += step;
 746       } while (from != to + step);
 747     }
 748     return result;
 749   }
 750
 751   /**
 752    * Returns a list of any mappings that are from or to the given (aligned or
 753    * dataset) sequence.
 754    *
 755    * @param sequence
 756    * @param mappings
 757    * @return
 758    */
 759   public static List<AlignedCodonFrame> findMappingsForSequence(
 760           SequenceI sequence, List<AlignedCodonFrame> mappings)
 761   {
 762     return findMappingsForSequenceAndOthers(sequence, mappings, null);
 763   }
 764
 765   /**
 766    * Returns a list of any mappings that are from or to the given (aligned or
 767    * dataset) sequence, optionally limited to mappings involving one of a given
 768    * list of sequences.
 769    *
 770    * @param sequence
 771    * @param mappings
 772    * @param filterList
 773    * @return
 774    */
 775   public static List<AlignedCodonFrame> findMappingsForSequenceAndOthers(
 776           SequenceI sequence, List<AlignedCodonFrame> mappings,
 777           List<SequenceI> filterList)
 778   {
 779     List<AlignedCodonFrame> result = new ArrayList<>();
 780     if (sequence == null || mappings == null)
 781     {
 782       return result;
 783     }
 784     for (AlignedCodonFrame mapping : mappings)
 785     {
 786       if (mapping.involvesSequence(sequence))
 787       {
 788         if (filterList != null)
 789         {
 790           for (SequenceI otherseq : filterList)
 791           {
 792             SequenceI otherDataset = otherseq.getDatasetSequence();
 793             if (otherseq == sequence
 794                     || otherseq == sequence.getDatasetSequence()
 795                     || (otherDataset != null && (otherDataset == sequence
 796                             || otherDataset == sequence
 797                                     .getDatasetSequence())))
 798             {
 799               // skip sequences in subset which directly relate to sequence
 800               continue;
 801             }
 802             if (mapping.involvesSequence(otherseq))
 803             {
 804               // selected a mapping contained in subselect alignment
 805               result.add(mapping);
 806               break;
 807             }
 808           }
 809         }
 810         else
 811         {
 812           result.add(mapping);
 813         }
 814       }
 815     }
 816     return result;
 817   }
 818
 819   /**
 820    * Returns the total length of the supplied ranges, which may be as single
 821    * [start, end] or multiple [start, end, start, end ...]
 822    *
 823    * @param ranges
 824    * @return
 825    */
 826   public static int getLength(List<int[]> ranges)
 827   {
 828     if (ranges == null)
 829     {
 830       return 0;
 831     }
 832     int length = 0;
 833     for (int[] range : ranges)
 834     {
 835       if (range.length % 2 != 0)
 836       {
 837         System.err.println(
 838                 "Error unbalance start/end ranges: " + ranges.toString());
 839         return 0;
 840       }
 841       for (int i = 0; i < range.length - 1; i += 2)
 842       {
 843         length += Math.abs(range[i + 1] - range[i]) + 1;
 844       }
 845     }
 846     return length;
 847   }
 848
 849   /**
 850    * Answers true if any range includes the given value
 851    *
 852    * @param ranges
 853    * @param value
 854    * @return
 855    */
 856   public static boolean contains(List<int[]> ranges, int value)
 857   {
 858     if (ranges == null)
 859     {
 860       return false;
 861     }
 862     for (int[] range : ranges)
 863     {
 864       if (range[1] >= range[0] && value >= range[0] && value <= range[1])
 865       {
 866         /*
 867          * value within ascending range
 868          */
 869         return true;
 870       }
 871       if (range[1] < range[0] && value <= range[0] && value >= range[1])
 872       {
 873         /*
 874          * value within descending range
 875          */
 876         return true;
 877       }
 878     }
 879     return false;
 880   }
 881
 882   /**
 883    * Removes a specified number of positions from the start of a ranges list.
 884    * For example, could be used to adjust cds ranges to allow for an incomplete
 885    * start codon. Subranges are removed completely, or their start positions
 886    * adjusted, until the required number of positions has been removed from the
 887    * range. Reverse strand ranges are supported. The input array is not
 888    * modified.
 889    *
 890    * @param removeCount
 891    * @param ranges
 892    *          an array of [start, end, start, end...] positions
 893    * @return a new array with the first removeCount positions removed
 894    */
 895   public static int[] removeStartPositions(int removeCount,
 896           final int[] ranges)
 897   {
 898     if (removeCount <= 0)
 899     {
 900       return ranges;
 901     }
 902
 903     int[] copy = Arrays.copyOf(ranges, ranges.length);
 904     int sxpos = -1;
 905     int cdspos = 0;
 906     for (int x = 0; x < copy.length && sxpos == -1; x += 2)
 907     {
 908       cdspos += Math.abs(copy[x + 1] - copy[x]) + 1;
 909       if (removeCount < cdspos)
 910       {
 911         /*
 912          * we have removed enough, time to finish
 913          */
 914         sxpos = x;
 915
 916         /*
 917          * increment start of first exon, or decrement if reverse strand
 918          */
 919         if (copy[x] <= copy[x + 1])
 920         {
 921           copy[x] = copy[x + 1] - cdspos + removeCount + 1;
 922         }
 923         else
 924         {
 925           copy[x] = copy[x + 1] + cdspos - removeCount - 1;
 926         }
 927         break;
 928       }
 929     }
 930
 931     if (sxpos > 0)
 932     {
 933       /*
 934        * we dropped at least one entire sub-range - compact the array
 935        */
 936       int[] nxon = new int[copy.length - sxpos];
 937       System.arraycopy(copy, sxpos, nxon, 0, copy.length - sxpos);
 938       return nxon;
 939     }
 940     return copy;
 941   }
 942
 943   /**
 944    * Answers true if range's start-end positions include those of queryRange,
 945    * where either range might be in reverse direction, else false
 946    *
 947    * @param range
 948    *          a start-end range
 949    * @param queryRange
 950    *          a candidate subrange of range (start2-end2)
 951    * @return
 952    */
 953   public static boolean rangeContains(int[] range, int[] queryRange)
 954   {
 955     if (range == null || queryRange == null || range.length != 2
 956             || queryRange.length != 2)
 957     {
 958       /*
 959        * invalid arguments
 960        */
 961       return false;
 962     }
 963
 964     int min = Math.min(range[0], range[1]);
 965     int max = Math.max(range[0], range[1]);
 966
 967     return (min <= queryRange[0] && max >= queryRange[0]
 968             && min <= queryRange[1] && max >= queryRange[1]);
 969   }
 970
 971   /**
 972    * Removes the specified number of positions from the given ranges. Provided
 973    * to allow a stop codon to be stripped from a CDS sequence so that it matches
 974    * the peptide translation length.
 975    *
 976    * @param positions
 977    * @param ranges
 978    *          a list of (single) [start, end] ranges
 979    * @return
 980    */
 981   public static void removeEndPositions(int positions,
 982           List<int[]> ranges)
 983   {
 984     int toRemove = positions;
 985     Iterator<int[]> it = new ReverseListIterator<>(ranges);
 986     while (toRemove > 0)
 987     {
 988       int[] endRange = it.next();
 989       if (endRange.length != 2)
 990       {
 991         /*
 992          * not coded for [start1, end1, start2, end2, ...]
 993          */
 994         System.err
 995                 .println("MappingUtils.removeEndPositions doesn't handle multiple  ranges");
 996         return;
 997       }
 998
 999       int length = endRange[1] - endRange[0] + 1;
1000       if (length <= 0)
1001       {
1002         /*
1003          * not coded for a reverse strand range (end < start)
1004          */
1005         System.err
1006                 .println("MappingUtils.removeEndPositions doesn't handle reverse strand");
1007         return;
1008       }
1009       if (length > toRemove)
1010       {
1011         endRange[1] -= toRemove;
1012         toRemove = 0;
1013       }
1014       else
1015       {
1016         toRemove -= length;
1017         it.remove();
1018       }
1019     }
1020   }
1021 }