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.ArrayList;
  27 import java.util.List;
  28
  29 /**
  30  * Stores mapping between the columns of a protein alignment and a DNA alignment
  31  * and a list of individual codon to amino acid mappings between sequences.
  32  */
  33 public class AlignedCodonFrame
  34 {
  35
  36   /*
  37    * Data bean to hold mappings from one sequence to another
  38    */
  39   private class SequenceToSequenceMapping
  40   {
  41     private SequenceI fromSeq;
  42
  43     private Mapping mapping;
  44
  45     SequenceToSequenceMapping(SequenceI from, Mapping map)
  46     {
  47       this.fromSeq = from;
  48       this.mapping = map;
  49     }
  50
  51     /**
  52      * Readable representation for debugging only, not guaranteed not to change
  53      */
  54     @Override
  55     public String toString()
  56     {
  57       return String.format("From %s %s", fromSeq.getName(),
  58               mapping.toString());
  59     }
  60   }
  61
  62   private List<SequenceToSequenceMapping> mappings;
  63
  64   /**
  65    * Constructor
  66    */
  67   public AlignedCodonFrame()
  68   {
  69     mappings = new ArrayList<SequenceToSequenceMapping>();
  70   }
  71
  72   /**
  73    * Adds a mapping between the dataset sequences for the associated dna and
  74    * protein sequence objects
  75    *
  76    * @param dnaseq
  77    * @param aaseq
  78    * @param map
  79    */
  80   public void addMap(SequenceI dnaseq, SequenceI aaseq, MapList map)
  81   {
  82     // JBPNote DEBUG! THIS !
  83     // dnaseq.transferAnnotation(aaseq, mp);
  84     // aaseq.transferAnnotation(dnaseq, new Mapping(map.getInverse()));
  85
  86     SequenceI fromSeq = (dnaseq.getDatasetSequence() == null) ? dnaseq
  87             : dnaseq.getDatasetSequence();
  88     SequenceI toSeq = (aaseq.getDatasetSequence() == null) ? aaseq : aaseq
  89             .getDatasetSequence();
  90
  91     /*
  92      * if we already hold a mapping between these sequences, just add to it
  93      */
  94     for (SequenceToSequenceMapping ssm : mappings)
  95     {
  96       if (ssm.fromSeq == fromSeq && ssm.mapping.to == toSeq)
  97       {
  98         ssm.mapping.map.addMapList(map);
  99         return;
 100       }
 101     }
 102
 103     /*
 104      * otherwise, add a new sequence mapping
 105      */
 106     Mapping mp = new Mapping(toSeq, map);
 107     mappings.add(new SequenceToSequenceMapping(fromSeq, mp));
 108   }
 109
 110   public SequenceI[] getdnaSeqs()
 111   {
 112     // TODO return a list instead?
 113     // return dnaSeqs;
 114     List<SequenceI> seqs = new ArrayList<SequenceI>();
 115     for (SequenceToSequenceMapping ssm : mappings)
 116     {
 117       seqs.add(ssm.fromSeq);
 118     }
 119     return seqs.toArray(new SequenceI[seqs.size()]);
 120   }
 121
 122   public SequenceI[] getAaSeqs()
 123   {
 124     // TODO not used - remove?
 125     List<SequenceI> seqs = new ArrayList<SequenceI>();
 126     for (SequenceToSequenceMapping ssm : mappings)
 127     {
 128       seqs.add(ssm.mapping.to);
 129     }
 130     return seqs.toArray(new SequenceI[seqs.size()]);
 131   }
 132
 133   public MapList[] getdnaToProt()
 134   {
 135     List<MapList> maps = new ArrayList<MapList>();
 136     for (SequenceToSequenceMapping ssm : mappings)
 137     {
 138       maps.add(ssm.mapping.map);
 139     }
 140     return maps.toArray(new MapList[maps.size()]);
 141   }
 142
 143   public Mapping[] getProtMappings()
 144   {
 145     List<Mapping> maps = new ArrayList<Mapping>();
 146     for (SequenceToSequenceMapping ssm : mappings)
 147     {
 148       maps.add(ssm.mapping);
 149     }
 150     return maps.toArray(new Mapping[maps.size()]);
 151   }
 152
 153   /**
 154    * Returns the first mapping found which is to or from the given sequence, or
 155    * null.
 156    *
 157    * @param seq
 158    * @return
 159    */
 160   public Mapping getMappingForSequence(SequenceI seq)
 161   {
 162     SequenceI seqDs = seq.getDatasetSequence();
 163     seqDs = seqDs != null ? seqDs : seq;
 164
 165     for (SequenceToSequenceMapping ssm : mappings)
 166     {
 167       if (ssm.fromSeq == seqDs || ssm.mapping.to == seqDs)
 168       {
 169         return ssm.mapping;
 170       }
 171     }
 172     return null;
 173   }
 174
 175   /**
 176    * Return the corresponding aligned or dataset aa sequence for given dna
 177    * sequence, null if not found.
 178    *
 179    * @param sequenceRef
 180    * @return
 181    */
 182   public SequenceI getAaForDnaSeq(SequenceI dnaSeqRef)
 183   {
 184     SequenceI dnads = dnaSeqRef.getDatasetSequence();
 185     for (SequenceToSequenceMapping ssm : mappings)
 186     {
 187       if (ssm.fromSeq == dnaSeqRef || ssm.fromSeq == dnads)
 188       {
 189         return ssm.mapping.to;
 190       }
 191     }
 192     return null;
 193   }
 194
 195   /**
 196    *
 197    * @param sequenceRef
 198    * @return null or corresponding aaSeq entry for dnaSeq entry
 199    */
 200   public SequenceI getDnaForAaSeq(SequenceI aaSeqRef)
 201   {
 202     SequenceI aads = aaSeqRef.getDatasetSequence();
 203     for (SequenceToSequenceMapping ssm : mappings)
 204     {
 205       if (ssm.mapping.to == aaSeqRef || ssm.mapping.to == aads)
 206       {
 207         return ssm.fromSeq;
 208       }
 209     }
 210     return null;
 211   }
 212
 213   /**
 214    * test to see if codon frame involves seq in any way
 215    *
 216    * @param seq
 217    *          a nucleotide or protein sequence
 218    * @return true if a mapping exists to or from this sequence to any translated
 219    *         sequence
 220    */
 221   public boolean involvesSequence(SequenceI seq)
 222   {
 223     return getAaForDnaSeq(seq) != null || getDnaForAaSeq(seq) != null;
 224   }
 225
 226   /**
 227    * Add search results for regions in other sequences that translate or are
 228    * translated from a particular position in seq
 229    *
 230    * @param seq
 231    * @param index
 232    *          position in seq
 233    * @param results
 234    *          where highlighted regions go
 235    */
 236   public void markMappedRegion(SequenceI seq, int index,
 237           SearchResults results)
 238   {
 239     int[] codon;
 240     SequenceI ds = seq.getDatasetSequence();
 241     for (SequenceToSequenceMapping ssm : mappings)
 242     {
 243       if (ssm.fromSeq == seq || ssm.fromSeq == ds)
 244       {
 245         codon = ssm.mapping.map.locateInTo(index, index);
 246         if (codon != null)
 247         {
 248           for (int i = 0; i < codon.length; i += 2)
 249           {
 250             results.addResult(ssm.mapping.to, codon[i], codon[i + 1]);
 251           }
 252         }
 253       }
 254       else if (ssm.mapping.to == seq || ssm.mapping.to == ds)
 255       {
 256         {
 257           codon = ssm.mapping.map.locateInFrom(index, index);
 258           if (codon != null)
 259           {
 260             for (int i = 0; i < codon.length; i += 2)
 261             {
 262               results.addResult(ssm.fromSeq, codon[i], codon[i + 1]);
 263             }
 264           }
 265         }
 266       }
 267     }
 268   }
 269
 270   /**
 271    * Returns the DNA codon positions (base 1) for the given position (base 1) in
 272    * a mapped protein sequence, or null if no mapping is found.
 273    *
 274    * Intended for use in aligning cDNA to match aligned protein. Only the first
 275    * mapping found is returned, so not suitable for use if multiple protein
 276    * sequences are mapped to the same cDNA (but aligning cDNA as protein is
 277    * ill-defined for this case anyway).
 278    *
 279    * @param seq
 280    *          the DNA dataset sequence
 281    * @param aaPos
 282    *          residue position (base 1) in a protein sequence
 283    * @return
 284    */
 285   public int[] getDnaPosition(SequenceI seq, int aaPos)
 286   {
 287     /*
 288      * Adapted from markMappedRegion().
 289      */
 290     MapList ml = null;
 291     int i = 0;
 292     for (SequenceToSequenceMapping ssm : mappings)
 293     {
 294       if (ssm.fromSeq == seq)
 295       {
 296         ml = getdnaToProt()[i];
 297         break;
 298       }
 299       i++;
 300     }
 301     return ml == null ? null : ml.locateInFrom(aaPos, aaPos);
 302   }
 303
 304   /**
 305    * Convenience method to return the first aligned sequence in the given
 306    * alignment whose dataset has a mapping with the given (aligned or dataset)
 307    * sequence.
 308    *
 309    * @param seq
 310    *
 311    * @param al
 312    * @return
 313    */
 314   public SequenceI findAlignedSequence(SequenceI seq, AlignmentI al)
 315   {
 316     /*
 317      * Search mapped protein ('to') sequences first.
 318      */
 319     for (SequenceToSequenceMapping ssm : mappings)
 320     {
 321       if (ssm.fromSeq == seq || ssm.fromSeq == seq.getDatasetSequence())
 322       {
 323         for (SequenceI sourceAligned : al.getSequences())
 324         {
 325           if (ssm.mapping.to == sourceAligned.getDatasetSequence()
 326                   || ssm.mapping.to == sourceAligned)
 327           {
 328             return sourceAligned;
 329           }
 330         }
 331       }
 332     }
 333
 334     /*
 335      * Then try mapped dna sequences.
 336      */
 337     for (SequenceToSequenceMapping ssm : mappings)
 338     {
 339       if (ssm.mapping.to == seq
 340               || ssm.mapping.to == seq.getDatasetSequence())
 341       {
 342         for (SequenceI sourceAligned : al.getSequences())
 343         {
 344           if (ssm.fromSeq == sourceAligned.getDatasetSequence())
 345           {
 346             return sourceAligned;
 347           }
 348         }
 349       }
 350     }
 351
 352     return null;
 353   }
 354
 355   /**
 356    * Returns the region in the target sequence's dataset that is mapped to the
 357    * given position (base 1) in the query sequence's dataset. The region is a
 358    * set of start/end position pairs.
 359    *
 360    * @param target
 361    * @param query
 362    * @param queryPos
 363    * @return
 364    */
 365   public int[] getMappedRegion(SequenceI target, SequenceI query,
 366           int queryPos)
 367   {
 368     SequenceI targetDs = target.getDatasetSequence() == null ? target
 369             : target.getDatasetSequence();
 370     SequenceI queryDs = query.getDatasetSequence() == null ? query : query
 371             .getDatasetSequence();
 372     if (targetDs == null || queryDs == null /*|| dnaToProt == null*/)
 373     {
 374       return null;
 375     }
 376     for (SequenceToSequenceMapping ssm : mappings)
 377     {
 378       /*
 379        * try mapping from target to query
 380        */
 381       if (ssm.fromSeq == targetDs && ssm.mapping.to == queryDs)
 382       {
 383         int[] codon = ssm.mapping.map.locateInFrom(queryPos, queryPos);
 384         if (codon != null)
 385         {
 386           return codon;
 387         }
 388       }
 389       /*
 390        * else try mapping from query to target
 391        */
 392       else if (ssm.fromSeq == queryDs && ssm.mapping.to == targetDs)
 393       {
 394         int[] codon = ssm.mapping.map.locateInTo(queryPos, queryPos);
 395         if (codon != null)
 396         {
 397           return codon;
 398         }
 399       }
 400     }
 401     return null;
 402   }
 403
 404   /**
 405    * Returns the mapped DNA codons for the given position in a protein sequence,
 406    * or null if no mapping is found. Returns a list of (e.g.) ['g', 'c', 't']
 407    * codons. There may be more than one codon mapped to the protein if (for
 408    * example), there are mappings to cDNA variants.
 409    *
 410    * @param protein
 411    *          the peptide dataset sequence
 412    * @param aaPos
 413    *          residue position (base 1) in the peptide sequence
 414    * @return
 415    */
 416   public List<char[]> getMappedCodons(SequenceI protein, int aaPos)
 417   {
 418     MapList ml = null;
 419     SequenceI dnaSeq = null;
 420     List<char[]> result = new ArrayList<char[]>();
 421
 422     for (SequenceToSequenceMapping ssm : mappings)
 423     {
 424       if (ssm.mapping.to == protein)
 425       {
 426         ml = ssm.mapping.map;
 427         dnaSeq = ssm.fromSeq;
 428
 429         int[] codonPos = ml.locateInFrom(aaPos, aaPos);
 430         if (codonPos == null)
 431         {
 432           return null;
 433         }
 434
 435         /*
 436          * Read off the mapped nucleotides (converting to position base 0)
 437          */
 438         codonPos = MappingUtils.flattenRanges(codonPos);
 439         char[] dna = dnaSeq.getSequence();
 440         int start = dnaSeq.getStart();
 441         result.add(new char[] { dna[codonPos[0] - start],
 442             dna[codonPos[1] - start], dna[codonPos[2] - start] });
 443       }
 444     }
 445     return result.isEmpty() ? null : result;
 446   }
 447
 448   /**
 449    * Returns any mappings found which are from the given sequence, and to
 450    * distinct sequences.
 451    *
 452    * @param seq
 453    * @return
 454    */
 455   public List<Mapping> getMappingsFromSequence(SequenceI seq)
 456   {
 457     List<Mapping> result = new ArrayList<Mapping>();
 458     List<SequenceI> related = new ArrayList<SequenceI>();
 459     SequenceI seqDs = seq.getDatasetSequence();
 460     seqDs = seqDs != null ? seqDs : seq;
 461
 462     for (SequenceToSequenceMapping ssm : mappings)
 463     {
 464       final Mapping mapping = ssm.mapping;
 465       if (ssm.fromSeq == seqDs)
 466       {
 467         if (!related.contains(mapping.to))
 468         {
 469           result.add(mapping);
 470           related.add(mapping.to);
 471         }
 472       }
 473     }
 474     return result;
 475   }
 476
 477   /**
 478    * Test whether the given sequence is substitutable for one or more dummy
 479    * sequences in this mapping
 480    *
 481    * @param map
 482    * @param seq
 483    * @return
 484    */
 485   public boolean isRealisableWith(SequenceI seq)
 486   {
 487     return realiseWith(seq, false) > 0;
 488   }
 489
 490   /**
 491    * Replace any matchable mapped dummy sequences with the given real one.
 492    * Returns the count of sequence mappings instantiated.
 493    *
 494    * @param seq
 495    * @return
 496    */
 497   public int realiseWith(SequenceI seq)
 498   {
 499     return realiseWith(seq, true);
 500   }
 501
 502   /**
 503    * Returns the number of mapped dummy sequences that could be replaced with
 504    * the given real sequence.
 505    *
 506    * @param seq
 507    *          a dataset sequence
 508    * @param doUpdate
 509    *          if true, performs replacements, else only counts
 510    * @return
 511    */
 512   protected int realiseWith(SequenceI seq, boolean doUpdate)
 513   {
 514     SequenceI ds = seq.getDatasetSequence() != null ? seq
 515             .getDatasetSequence() : seq;
 516     int count = 0;
 517
 518     /*
 519      * check for replaceable DNA ('map from') sequences
 520      */
 521     for (SequenceToSequenceMapping 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 && mapStart <= end)
 590               || (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 (SequenceToSequenceMapping 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 }