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 dataset sequence.
 307    *
 308    * @param seq
 309    *
 310    * @param al
 311    * @return
 312    */
 313   public SequenceI findAlignedSequence(SequenceI seq, AlignmentI al)
 314   {
 315     /*
 316      * Search mapped protein ('to') sequences first.
 317      */
 318     for (SequenceToSequenceMapping ssm : mappings)
 319     {
 320       if (ssm.fromSeq == seq)
 321       {
 322         for (SequenceI sourceAligned : al.getSequences())
 323         {
 324           if (ssm.mapping.to == sourceAligned.getDatasetSequence())
 325           {
 326             return sourceAligned;
 327           }
 328         }
 329       }
 330     }
 331
 332     /*
 333      * Then try mapped dna sequences.
 334      */
 335     for (SequenceToSequenceMapping ssm : mappings)
 336     {
 337       if (ssm.mapping.to == seq)
 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 : query
 368             .getDatasetSequence();
 369     if (targetDs == null || queryDs == null /*|| dnaToProt == null*/)
 370     {
 371       return null;
 372     }
 373     for (SequenceToSequenceMapping 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<char[]>();
 418
 419     for (SequenceToSequenceMapping ssm : mappings)
 420     {
 421       if (ssm.mapping.to == protein)
 422       {
 423         ml = ssm.mapping.map;
 424         dnaSeq = ssm.fromSeq;
 425
 426         int[] codonPos = ml.locateInFrom(aaPos, aaPos);
 427         if (codonPos == null)
 428         {
 429           return null;
 430         }
 431
 432         /*
 433          * Read off the mapped nucleotides (converting to position base 0)
 434          */
 435         codonPos = MappingUtils.flattenRanges(codonPos);
 436         char[] dna = dnaSeq.getSequence();
 437         int start = dnaSeq.getStart();
 438         result.add(new char[] { dna[codonPos[0] - start],
 439             dna[codonPos[1] - start], dna[codonPos[2] - start] });
 440       }
 441     }
 442     return result.isEmpty() ? null : result;
 443   }
 444
 445   /**
 446    * Returns any mappings found which are to (or from) the given sequence, and
 447    * to distinct sequences.
 448    *
 449    * @param seq
 450    * @return
 451    */
 452   public List<Mapping> getMappingsForSequence(SequenceI seq)
 453   {
 454     List<Mapping> result = new ArrayList<Mapping>();
 455     List<SequenceI> related = new ArrayList<SequenceI>();
 456     SequenceI seqDs = seq.getDatasetSequence();
 457     seqDs = seqDs != null ? seqDs : seq;
 458
 459     for (SequenceToSequenceMapping ssm : mappings)
 460     {
 461       final Mapping mapping = ssm.mapping;
 462       if (ssm.fromSeq == seqDs || mapping.to == seqDs)
 463       {
 464         if (!related.contains(mapping.to))
 465         {
 466           result.add(mapping);
 467           related.add(mapping.to);
 468         }
 469       }
 470     }
 471     return result;
 472   }
 473
 474   /**
 475    * Test whether the given sequence is substitutable for one or more dummy
 476    * sequences in this mapping
 477    *
 478    * @param map
 479    * @param seq
 480    * @return
 481    */
 482   public boolean isRealisableWith(SequenceI seq)
 483   {
 484     return realiseWith(seq, false) > 0;
 485   }
 486
 487   /**
 488    * Replace any matchable mapped dummy sequences with the given real one.
 489    * Returns the count of sequence mappings instantiated.
 490    *
 491    * @param seq
 492    * @return
 493    */
 494   public int realiseWith(SequenceI seq)
 495   {
 496     return realiseWith(seq, true);
 497   }
 498
 499   /**
 500    * Returns the number of mapped dummy sequences that could be replaced with
 501    * the given real sequence.
 502    *
 503    * @param seq
 504    *          a dataset sequence
 505    * @param doUpdate
 506    *          if true, performs replacements, else only counts
 507    * @return
 508    */
 509   protected int realiseWith(SequenceI seq, boolean doUpdate)
 510   {
 511     SequenceI ds = seq.getDatasetSequence() != null ? seq
 512             .getDatasetSequence() : seq;
 513     int count = 0;
 514
 515     /*
 516      * check for replaceable DNA ('map from') sequences
 517      */
 518     for (SequenceToSequenceMapping ssm : mappings)
 519     {
 520       SequenceI dna = ssm.fromSeq;
 521       if (dna instanceof SequenceDummy
 522               && dna.getName().equals(ds.getName()))
 523       {
 524         Mapping mapping = ssm.mapping;
 525         int mapStart = mapping.getMap().getFromLowest();
 526         int mapEnd = mapping.getMap().getFromHighest();
 527         boolean mappable = couldRealiseSequence(dna, ds, mapStart, mapEnd);
 528         if (mappable)
 529         {
 530           count++;
 531           if (doUpdate)
 532           {
 533             // TODO: new method ? ds.realise(dna);
 534             // might want to copy database refs as well
 535             ds.setSequenceFeatures(dna.getSequenceFeatures());
 536             // dnaSeqs[i] = ds;
 537             ssm.fromSeq = ds;
 538             System.out.println("Realised mapped sequence " + ds.getName());
 539           }
 540         }
 541       }
 542
 543       /*
 544        * check for replaceable protein ('map to') sequences
 545        */
 546       Mapping mapping = ssm.mapping;
 547       SequenceI prot = mapping.getTo();
 548       int mapStart = mapping.getMap().getToLowest();
 549       int mapEnd = mapping.getMap().getToHighest();
 550       boolean mappable = couldRealiseSequence(prot, ds, mapStart, mapEnd);
 551       if (mappable)
 552       {
 553         count++;
 554         if (doUpdate)
 555         {
 556           // TODO: new method ? ds.realise(dna);
 557           // might want to copy database refs as well
 558           ds.setSequenceFeatures(dna.getSequenceFeatures());
 559           ssm.mapping.setTo(ds);
 560         }
 561       }
 562     }
 563     return count;
 564   }
 565
 566   /**
 567    * Helper method to test whether a 'real' sequence could replace a 'dummy'
 568    * sequence in the map. The criteria are that they have the same name, and
 569    * that the mapped region overlaps the candidate sequence.
 570    *
 571    * @param existing
 572    * @param replacement
 573    * @param mapStart
 574    * @param mapEnd
 575    * @return
 576    */
 577   protected static boolean couldRealiseSequence(SequenceI existing,
 578           SequenceI replacement, int mapStart, int mapEnd)
 579   {
 580     if (existing instanceof SequenceDummy
 581             && !(replacement instanceof SequenceDummy)
 582             && existing.getName().equals(replacement.getName()))
 583     {
 584       int start = replacement.getStart();
 585       int end = replacement.getEnd();
 586       boolean mappingOverlapsSequence = (mapStart >= start && mapStart <= end)
 587               || (mapEnd >= start && mapEnd <= end);
 588       if (mappingOverlapsSequence)
 589       {
 590         return true;
 591       }
 592     }
 593     return false;
 594   }
 595
 596   /**
 597    * Change any mapping to the given sequence to be to its dataset sequence
 598    * instead. For use when mappings are created before their referenced
 599    * sequences are instantiated, for example when parsing GFF data.
 600    *
 601    * @param seq
 602    */
 603   public void updateToDataset(SequenceI seq)
 604   {
 605     if (seq == null || seq.getDatasetSequence() == null)
 606     {
 607       return;
 608     }
 609     SequenceI ds = seq.getDatasetSequence();
 610
 611     for (SequenceToSequenceMapping ssm : mappings)
 612     /*
 613      * 'from' sequences
 614      */
 615     {
 616       if (ssm.fromSeq == seq)
 617       {
 618         ssm.fromSeq = ds;
 619       }
 620
 621       /*
 622        * 'to' sequences
 623        */
 624       if (ssm.mapping.to == seq)
 625       {
 626         ssm.mapping.to = ds;
 627       }
 628     }
 629   }
 630
 631   /**
 632    * Answers true if this object contains no mappings
 633    *
 634    * @return
 635    */
 636   public boolean isEmpty()
 637   {
 638     return mappings.isEmpty();
 639   }
 640 }