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