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