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