src/jalview/datamodel/AlignedCodonFrame.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.datamodel;
  22
  23 import jalview.util.MapList;
  24 import jalview.util.MappingUtils;
  25
  26 import java.util.AbstractList;
  27 import java.util.ArrayList;
  28 import java.util.List;
  29
  30 /**
  31  * Stores mapping between the columns of a protein alignment and a DNA alignment
  32  * and a list of individual codon to amino acid mappings between sequences.
  33  */
  34 public class AlignedCodonFrame
  35 {
  36
  37   private List<SequenceMapping> mappings;
  38
  39   /**
  40    * Constructor
  41    */
  42   public AlignedCodonFrame()
  43   {
  44     mappings = new ArrayList<>();
  45   }
  46
  47   /**
  48    * Adds a mapping between the dataset sequences for the associated dna and
  49    * protein sequence objects
  50    *
  51    * @param dnaseq
  52    * @param aaseq
  53    * @param map
  54    */
  55   public void addMap(SequenceI dnaseq, SequenceI aaseq, MapList map)
  56   {
  57     addMap(dnaseq, aaseq, map, null);
  58   }
  59
  60   /**
  61    * Adds a mapping between the dataset sequences for the associated dna and
  62    * protein sequence objects
  63    *
  64    * @param dnaseq
  65    * @param aaseq
  66    * @param map
  67    * @param mapFromId
  68    */
  69   public void addMap(SequenceI dnaseq, SequenceI aaseq, MapList map,
  70           String mapFromId)
  71   {
  72     // JBPNote DEBUG! THIS !
  73     // dnaseq.transferAnnotation(aaseq, mp);
  74     // aaseq.transferAnnotation(dnaseq, new Mapping(map.getInverse()));
  75
  76     SequenceI fromSeq = (dnaseq.getDatasetSequence() == null) ? dnaseq
  77             : dnaseq.getDatasetSequence();
  78     SequenceI toSeq = (aaseq.getDatasetSequence() == null) ? aaseq
  79             : aaseq.getDatasetSequence();
  80
  81     /*
  82      * if we already hold a mapping between these sequences, just add to it
  83      * note that 'adding' a duplicate map does nothing; this protects against
  84      * creating duplicate mappings in AlignedCodonFrame
  85      */
  86     for (SequenceMapping ssm : mappings)
  87     {
  88       if (ssm.fromSeq == fromSeq && ssm.mapping.to == toSeq)
  89       {
  90         ssm.mapping.map.addMapList(map);
  91         return;
  92       }
  93     }
  94
  95     /*
  96      * otherwise, add a new sequence mapping
  97      */
  98     Mapping mp = new Mapping(toSeq, map);
  99     mp.setMappedFromId(mapFromId);
 100     mappings.add(new SequenceMapping(fromSeq, mp));
 101   }
 102
 103   public SequenceI[] getdnaSeqs()
 104   {
 105     // TODO return a list instead?
 106     // return dnaSeqs;
 107     List<SequenceI> seqs = new ArrayList<>();
 108     for (SequenceMapping ssm : mappings)
 109     {
 110       seqs.add(ssm.fromSeq);
 111     }
 112     return seqs.toArray(new SequenceI[seqs.size()]);
 113   }
 114
 115   public SequenceI[] getAaSeqs()
 116   {
 117     // TODO not used - remove?
 118     List<SequenceI> seqs = new ArrayList<>();
 119     for (SequenceMapping ssm : mappings)
 120     {
 121       seqs.add(ssm.mapping.to);
 122     }
 123     return seqs.toArray(new SequenceI[seqs.size()]);
 124   }
 125
 126   public MapList[] getdnaToProt()
 127   {
 128     List<MapList> maps = new ArrayList<>();
 129     for (SequenceMapping ssm : mappings)
 130     {
 131       maps.add(ssm.mapping.map);
 132     }
 133     return maps.toArray(new MapList[maps.size()]);
 134   }
 135
 136   public Mapping[] getProtMappings()
 137   {
 138     List<Mapping> maps = new ArrayList<>();
 139     for (SequenceMapping ssm : mappings)
 140     {
 141       maps.add(ssm.mapping);
 142     }
 143     return maps.toArray(new Mapping[maps.size()]);
 144   }
 145
 146   /**
 147    * Returns the first mapping found which is to or from the given sequence, or
 148    * null.
 149    *
 150    * @param seq
 151    * @return
 152    */
 153   public Mapping getMappingForSequence(SequenceI seq, boolean cdsOnly)
 154   {
 155     SequenceI seqDs = seq.getDatasetSequence();
 156     seqDs = seqDs != null ? seqDs : seq;
 157
 158     for (SequenceMapping ssm : mappings)
 159     {
 160       if (ssm.fromSeq == seqDs || ssm.mapping.to == seqDs)
 161       {
 162         if (!cdsOnly || ssm.fromSeq.getName().startsWith("CDS")
 163                 || ssm.mapping.to.getName().startsWith("CDS"))
 164         {
 165           return ssm.mapping;
 166         }
 167       }
 168     }
 169     return null;
 170   }
 171
 172   /**
 173    * Return the corresponding aligned or dataset aa sequence for given dna
 174    * sequence, null if not found.
 175    *
 176    * @param sequenceRef
 177    * @return
 178    */
 179   public SequenceI getAaForDnaSeq(SequenceI dnaSeqRef)
 180   {
 181     SequenceI dnads = dnaSeqRef.getDatasetSequence();
 182     for (SequenceMapping ssm : mappings)
 183     {
 184       if (ssm.fromSeq == dnaSeqRef || ssm.fromSeq == dnads)
 185       {
 186         return ssm.mapping.to;
 187       }
 188     }
 189     return null;
 190   }
 191
 192   /**
 193    *
 194    * @param sequenceRef
 195    * @return null or corresponding aaSeq entry for dnaSeq entry
 196    */
 197   public SequenceI getDnaForAaSeq(SequenceI aaSeqRef)
 198   {
 199     SequenceI aads = aaSeqRef.getDatasetSequence();
 200     for (SequenceMapping ssm : mappings)
 201     {
 202       if (ssm.mapping.to == aaSeqRef || ssm.mapping.to == aads)
 203       {
 204         return ssm.fromSeq;
 205       }
 206     }
 207     return null;
 208   }
 209
 210   /**
 211    * test to see if codon frame involves seq in any way
 212    *
 213    * @param seq
 214    *          a nucleotide or protein sequence
 215    * @return true if a mapping exists to or from this sequence to any translated
 216    *         sequence
 217    */
 218   public boolean involvesSequence(SequenceI seq)
 219   {
 220     return getAaForDnaSeq(seq) != null || getDnaForAaSeq(seq) != null;
 221   }
 222
 223   /**
 224    * Add search results for regions in other sequences that translate or are
 225    * translated from a particular position in seq
 226    *
 227    * @param seq
 228    * @param index
 229    *          position in seq
 230    * @param results
 231    *          where highlighted regions go
 232    */
 233   public void markMappedRegion(SequenceI seq, int index,
 234           SearchResultsI results)
 235   {
 236     int[] codon;
 237     SequenceI ds = seq.getDatasetSequence();
 238     for (SequenceMapping ssm : mappings)
 239     {
 240       if (ssm.fromSeq == seq || ssm.fromSeq == ds)
 241       {
 242         codon = ssm.mapping.map.locateInTo(index, index);
 243         if (codon != null)
 244         {
 245           for (int i = 0; i < codon.length; i += 2)
 246           {
 247             results.addResult(ssm.mapping.to, codon[i], codon[i + 1]);
 248           }
 249         }
 250       }
 251       else if (ssm.mapping.to == seq || ssm.mapping.to == ds)
 252       {
 253         {
 254           codon = ssm.mapping.map.locateInFrom(index, index);
 255           if (codon != null)
 256           {
 257             for (int i = 0; i < codon.length; i += 2)
 258             {
 259               results.addResult(ssm.fromSeq, codon[i], codon[i + 1]);
 260             }
 261           }
 262         }
 263       }
 264     }
 265   }
 266
 267   /**
 268    * Returns the DNA codon positions (base 1) for the given position (base 1) in
 269    * a mapped protein sequence, or null if no mapping is found.
 270    *
 271    * Intended for use in aligning cDNA to match aligned protein. Only the first
 272    * mapping found is returned, so not suitable for use if multiple protein
 273    * sequences are mapped to the same cDNA (but aligning cDNA as protein is
 274    * ill-defined for this case anyway).
 275    *
 276    * @param seq
 277    *          the DNA dataset sequence
 278    * @param aaPos
 279    *          residue position (base 1) in a protein sequence
 280    * @return
 281    */
 282   public int[] getDnaPosition(SequenceI seq, int aaPos)
 283   {
 284     /*
 285      * Adapted from markMappedRegion().
 286      */
 287     MapList ml = null;
 288     int i = 0;
 289     for (SequenceMapping ssm : mappings)
 290     {
 291       if (ssm.fromSeq == seq)
 292       {
 293         ml = getdnaToProt()[i];
 294         break;
 295       }
 296       i++;
 297     }
 298     return ml == null ? null : ml.locateInFrom(aaPos, aaPos);
 299   }
 300
 301   /**
 302    * Convenience method to return the first aligned sequence in the given
 303    * alignment whose dataset has a mapping with the given (aligned or dataset)
 304    * sequence.
 305    *
 306    * @param seq
 307    *
 308    * @param al
 309    * @return
 310    */
 311   public SequenceI findAlignedSequence(SequenceI seq, AlignmentI al)
 312   {
 313     /*
 314      * Search mapped protein ('to') sequences first.
 315      */
 316     for (SequenceMapping ssm : mappings)
 317     {
 318       if (ssm.fromSeq == seq || ssm.fromSeq == seq.getDatasetSequence())
 319       {
 320         for (SequenceI sourceAligned : al.getSequences())
 321         {
 322           if (ssm.mapping.to == sourceAligned.getDatasetSequence()
 323                   || ssm.mapping.to == sourceAligned)
 324           {
 325             return sourceAligned;
 326           }
 327         }
 328       }
 329     }
 330
 331     /*
 332      * Then try mapped dna sequences.
 333      */
 334     for (SequenceMapping ssm : mappings)
 335     {
 336       if (ssm.mapping.to == seq
 337               || ssm.mapping.to == seq.getDatasetSequence())
 338       {
 339         for (SequenceI sourceAligned : al.getSequences())
 340         {
 341           if (ssm.fromSeq == sourceAligned.getDatasetSequence())
 342           {
 343             return sourceAligned;
 344           }
 345         }
 346       }
 347     }
 348
 349     return null;
 350   }
 351
 352   /**
 353    * Returns the region in the target sequence's dataset that is mapped to the
 354    * given position (base 1) in the query sequence's dataset. The region is a
 355    * set of start/end position pairs.
 356    *
 357    * @param target
 358    * @param query
 359    * @param queryPos
 360    * @return
 361    */
 362   public int[] getMappedRegion(SequenceI target, SequenceI query,
 363           int queryPos)
 364   {
 365     SequenceI targetDs = target.getDatasetSequence() == null ? target
 366             : target.getDatasetSequence();
 367     SequenceI queryDs = query.getDatasetSequence() == null ? query
 368             : query.getDatasetSequence();
 369     if (targetDs == null || queryDs == null /*|| dnaToProt == null*/)
 370     {
 371       return null;
 372     }
 373     for (SequenceMapping ssm : mappings)
 374     {
 375       /*
 376        * try mapping from target to query
 377        */
 378       if (ssm.fromSeq == targetDs && ssm.mapping.to == queryDs)
 379       {
 380         int[] codon = ssm.mapping.map.locateInFrom(queryPos, queryPos);
 381         if (codon != null)
 382         {
 383           return codon;
 384         }
 385       }
 386       /*
 387        * else try mapping from query to target
 388        */
 389       else if (ssm.fromSeq == queryDs && ssm.mapping.to == targetDs)
 390       {
 391         int[] codon = ssm.mapping.map.locateInTo(queryPos, queryPos);
 392         if (codon != null)
 393         {
 394           return codon;
 395         }
 396       }
 397     }
 398     return null;
 399   }
 400
 401   /**
 402    * Returns the mapped DNA codons for the given position in a protein sequence,
 403    * or null if no mapping is found. Returns a list of (e.g.) ['g', 'c', 't']
 404    * codons. There may be more than one codon mapped to the protein if (for
 405    * example), there are mappings to cDNA variants.
 406    *
 407    * @param protein
 408    *          the peptide dataset sequence
 409    * @param aaPos
 410    *          residue position (base 1) in the peptide sequence
 411    * @return
 412    */
 413   public List<char[]> getMappedCodons(SequenceI protein, int aaPos)
 414   {
 415     MapList ml = null;
 416     SequenceI dnaSeq = null;
 417     List<char[]> result = new ArrayList<>();
 418
 419     for (SequenceMapping ssm : mappings)
 420     {
 421       if (ssm.mapping.to == protein
 422               && ssm.mapping.getMap().getFromRatio() == 3)
 423       {
 424         ml = ssm.mapping.map;
 425         dnaSeq = ssm.fromSeq;
 426
 427         int[] codonPos = ml.locateInFrom(aaPos, aaPos);
 428         if (codonPos == null)
 429         {
 430           return null;
 431         }
 432
 433         /*
 434          * Read off the mapped nucleotides (converting to position base 0)
 435          */
 436         codonPos = MappingUtils.flattenRanges(codonPos);
 437         int start = dnaSeq.getStart();
 438         char c1 = dnaSeq.getCharAt(codonPos[0] - start);
 439         char c2 = dnaSeq.getCharAt(codonPos[1] - start);
 440         char c3 = dnaSeq.getCharAt(codonPos[2] - start);
 441         result.add(new char[] { c1, c2, c3 });
 442       }
 443     }
 444     return result.isEmpty() ? null : result;
 445   }
 446
 447   /**
 448    * Returns any mappings found which are from the given sequence, and to
 449    * distinct sequences.
 450    *
 451    * @param seq
 452    * @return
 453    */
 454   public List<Mapping> getMappingsFromSequence(SequenceI seq)
 455   {
 456     List<Mapping> result = new ArrayList<>();
 457     List<SequenceI> related = new ArrayList<>();
 458     SequenceI seqDs = seq.getDatasetSequence();
 459     seqDs = seqDs != null ? seqDs : seq;
 460
 461     for (SequenceMapping ssm : mappings)
 462     {
 463       final Mapping mapping = ssm.mapping;
 464       if (ssm.fromSeq == seqDs)
 465       {
 466         if (!related.contains(mapping.to))
 467         {
 468           result.add(mapping);
 469           related.add(mapping.to);
 470         }
 471       }
 472     }
 473     return result;
 474   }
 475
 476   /**
 477    * Test whether the given sequence is substitutable for one or more dummy
 478    * sequences in this mapping
 479    *
 480    * @param map
 481    * @param seq
 482    * @return
 483    */
 484   public boolean isRealisableWith(SequenceI seq)
 485   {
 486     return realiseWith(seq, false) > 0;
 487   }
 488
 489   /**
 490    * Replace any matchable mapped dummy sequences with the given real one.
 491    * Returns the count of sequence mappings instantiated.
 492    *
 493    * @param seq
 494    * @return
 495    */
 496   public int realiseWith(SequenceI seq)
 497   {
 498     return realiseWith(seq, true);
 499   }
 500
 501   /**
 502    * Returns the number of mapped dummy sequences that could be replaced with
 503    * the given real sequence.
 504    *
 505    * @param seq
 506    *          a dataset sequence
 507    * @param doUpdate
 508    *          if true, performs replacements, else only counts
 509    * @return
 510    */
 511   protected int realiseWith(SequenceI seq, boolean doUpdate)
 512   {
 513     SequenceI ds = seq.getDatasetSequence() != null
 514             ? seq.getDatasetSequence()
 515             : seq;
 516     int count = 0;
 517
 518     /*
 519      * check for replaceable DNA ('map from') sequences
 520      */
 521     for (SequenceMapping ssm : mappings)
 522     {
 523       SequenceI dna = ssm.fromSeq;
 524       if (dna instanceof SequenceDummy
 525               && dna.getName().equals(ds.getName()))
 526       {
 527         Mapping mapping = ssm.mapping;
 528         int mapStart = mapping.getMap().getFromLowest();
 529         int mapEnd = mapping.getMap().getFromHighest();
 530         boolean mappable = couldRealiseSequence(dna, ds, mapStart, mapEnd);
 531         if (mappable)
 532         {
 533           count++;
 534           if (doUpdate)
 535           {
 536             // TODO: new method ? ds.realise(dna);
 537             // might want to copy database refs as well
 538             ds.setSequenceFeatures(dna.getSequenceFeatures());
 539             // dnaSeqs[i] = ds;
 540             ssm.fromSeq = ds;
 541             System.out.println("Realised mapped sequence " + ds.getName());
 542           }
 543         }
 544       }
 545
 546       /*
 547        * check for replaceable protein ('map to') sequences
 548        */
 549       Mapping mapping = ssm.mapping;
 550       SequenceI prot = mapping.getTo();
 551       int mapStart = mapping.getMap().getToLowest();
 552       int mapEnd = mapping.getMap().getToHighest();
 553       boolean mappable = couldRealiseSequence(prot, ds, mapStart, mapEnd);
 554       if (mappable)
 555       {
 556         count++;
 557         if (doUpdate)
 558         {
 559           // TODO: new method ? ds.realise(dna);
 560           // might want to copy database refs as well
 561           ds.setSequenceFeatures(dna.getSequenceFeatures());
 562           ssm.mapping.setTo(ds);
 563         }
 564       }
 565     }
 566     return count;
 567   }
 568
 569   /**
 570    * Helper method to test whether a 'real' sequence could replace a 'dummy'
 571    * sequence in the map. The criteria are that they have the same name, and
 572    * that the mapped region overlaps the candidate sequence.
 573    *
 574    * @param existing
 575    * @param replacement
 576    * @param mapStart
 577    * @param mapEnd
 578    * @return
 579    */
 580   protected static boolean couldRealiseSequence(SequenceI existing,
 581           SequenceI replacement, int mapStart, int mapEnd)
 582   {
 583     if (existing instanceof SequenceDummy
 584             && !(replacement instanceof SequenceDummy)
 585             && existing.getName().equals(replacement.getName()))
 586     {
 587       int start = replacement.getStart();
 588       int end = replacement.getEnd();
 589       boolean mappingOverlapsSequence = (mapStart >= start
 590               && mapStart <= end) || (mapEnd >= start && mapEnd <= end);
 591       if (mappingOverlapsSequence)
 592       {
 593         return true;
 594       }
 595     }
 596     return false;
 597   }
 598
 599   /**
 600    * Change any mapping to the given sequence to be to its dataset sequence
 601    * instead. For use when mappings are created before their referenced
 602    * sequences are instantiated, for example when parsing GFF data.
 603    *
 604    * @param seq
 605    */
 606   public void updateToDataset(SequenceI seq)
 607   {
 608     if (seq == null || seq.getDatasetSequence() == null)
 609     {
 610       return;
 611     }
 612     SequenceI ds = seq.getDatasetSequence();
 613
 614     for (SequenceMapping ssm : mappings)
 615     /*
 616      * 'from' sequences
 617      */
 618     {
 619       if (ssm.fromSeq == seq)
 620       {
 621         ssm.fromSeq = ds;
 622       }
 623
 624       /*
 625        * 'to' sequences
 626        */
 627       if (ssm.mapping.to == seq)
 628       {
 629         ssm.mapping.to = ds;
 630       }
 631     }
 632   }
 633
 634   /**
 635    * Answers true if this object contains no mappings
 636    *
 637    * @return
 638    */
 639   public boolean isEmpty()
 640   {
 641     return mappings.isEmpty();
 642   }
 643
 644   /**
 645    * Method for debug / inspection purposes only, may change in future
 646    */
 647   @Override
 648   public String toString()
 649   {
 650     return mappings == null ? "null" : mappings.toString();
 651   }
 652
 653   /**
 654    * Returns the first mapping found that is between 'fromSeq' and 'toSeq', or
 655    * null if none found
 656    *
 657    * @param fromSeq
 658    *          aligned or dataset sequence
 659    * @param toSeq
 660    *          aligned or dataset sequence
 661    * @return
 662    */
 663   public Mapping getMappingBetween(SequenceI fromSeq, SequenceI toSeq)
 664   {
 665     SequenceI dssFrom = fromSeq.getDatasetSequence() == null ? fromSeq
 666             : fromSeq.getDatasetSequence();
 667     SequenceI dssTo = toSeq.getDatasetSequence() == null ? toSeq
 668             : toSeq.getDatasetSequence();
 669
 670     for (SequenceMapping mapping : mappings)
 671     {
 672       SequenceI from = mapping.fromSeq;
 673       SequenceI to = mapping.mapping.to;
 674       if ((from == dssFrom && to == dssTo)
 675               || (from == dssTo && to == dssFrom))
 676       {
 677         return mapping.mapping;
 678       }
 679     }
 680     return null;
 681   }
 682
 683   /**
 684    * Returns a hashcode derived from the list of sequence mappings
 685    *
 686    * @see SequenceMapping#hashCode()
 687    * @see AbstractList#hashCode()
 688    */
 689   @Override
 690   public int hashCode()
 691   {
 692     return this.mappings.hashCode();
 693   }
 694
 695   /**
 696    * Two AlignedCodonFrame objects are equal if they hold the same ordered list
 697    * of mappings
 698    *
 699    * @see SequenceToSequenceMapping#
 700    */
 701   @Override
 702   public boolean equals(Object obj)
 703   {
 704     if (!(obj instanceof AlignedCodonFrame))
 705     {
 706       return false;
 707     }
 708     return this.mappings.equals(((AlignedCodonFrame) obj).mappings);
 709   }
 710
 711   public List<SequenceMapping> getMappings()
 712   {
 713     return mappings;
 714   }
 715 }