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                   || ssm.mapping.to == sourceAligned)
 326           {
 327             return sourceAligned;
 328           }
 329         }
 330       }
 331     }
 332
 333     /*
 334      * Then try mapped dna sequences.
 335      */
 336     for (SequenceToSequenceMapping ssm : mappings)
 337     {
 338       if (ssm.mapping.to == seq)
 339       {
 340         for (SequenceI sourceAligned : al.getSequences())
 341         {
 342           if (ssm.fromSeq == sourceAligned.getDatasetSequence())
 343           {
 344             return sourceAligned;
 345           }
 346         }
 347       }
 348     }
 349
 350     return null;
 351   }
 352
 353   /**
 354    * Returns the region in the target sequence's dataset that is mapped to the
 355    * given position (base 1) in the query sequence's dataset. The region is a
 356    * set of start/end position pairs.
 357    *
 358    * @param target
 359    * @param query
 360    * @param queryPos
 361    * @return
 362    */
 363   public int[] getMappedRegion(SequenceI target, SequenceI query,
 364           int queryPos)
 365   {
 366     SequenceI targetDs = target.getDatasetSequence() == null ? target
 367             : target.getDatasetSequence();
 368     SequenceI queryDs = query.getDatasetSequence() == null ? query : query
 369             .getDatasetSequence();
 370     if (targetDs == null || queryDs == null /*|| dnaToProt == null*/)
 371     {
 372       return null;
 373     }
 374     for (SequenceToSequenceMapping ssm : mappings)
 375     {
 376       /*
 377        * try mapping from target to query
 378        */
 379       if (ssm.fromSeq == targetDs && ssm.mapping.to == queryDs)
 380       {
 381         int[] codon = ssm.mapping.map.locateInFrom(queryPos, queryPos);
 382         if (codon != null)
 383         {
 384           return codon;
 385         }
 386       }
 387       /*
 388        * else try mapping from query to target
 389        */
 390       else if (ssm.fromSeq == queryDs && ssm.mapping.to == targetDs)
 391       {
 392         int[] codon = ssm.mapping.map.locateInTo(queryPos, queryPos);
 393         if (codon != null)
 394         {
 395           return codon;
 396         }
 397       }
 398     }
 399     return null;
 400   }
 401
 402   /**
 403    * Returns the mapped DNA codons for the given position in a protein sequence,
 404    * or null if no mapping is found. Returns a list of (e.g.) ['g', 'c', 't']
 405    * codons. There may be more than one codon mapped to the protein if (for
 406    * example), there are mappings to cDNA variants.
 407    *
 408    * @param protein
 409    *          the peptide dataset sequence
 410    * @param aaPos
 411    *          residue position (base 1) in the peptide sequence
 412    * @return
 413    */
 414   public List<char[]> getMappedCodons(SequenceI protein, int aaPos)
 415   {
 416     MapList ml = null;
 417     SequenceI dnaSeq = null;
 418     List<char[]> result = new ArrayList<char[]>();
 419
 420     for (SequenceToSequenceMapping ssm : mappings)
 421     {
 422       if (ssm.mapping.to == protein)
 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         char[] dna = dnaSeq.getSequence();
 438         int start = dnaSeq.getStart();
 439         result.add(new char[] { dna[codonPos[0] - start],
 440             dna[codonPos[1] - start], dna[codonPos[2] - start] });
 441       }
 442     }
 443     return result.isEmpty() ? null : result;
 444   }
 445
 446   /**
 447    * Returns any mappings found which are to (or from) the given sequence, and
 448    * to distinct sequences.
 449    *
 450    * @param seq
 451    * @return
 452    */
 453   public List<Mapping> getMappingsForSequence(SequenceI seq)
 454   {
 455     List<Mapping> result = new ArrayList<Mapping>();
 456     List<SequenceI> related = new ArrayList<SequenceI>();
 457     SequenceI seqDs = seq.getDatasetSequence();
 458     seqDs = seqDs != null ? seqDs : seq;
 459
 460     for (SequenceToSequenceMapping ssm : mappings)
 461     {
 462       final Mapping mapping = ssm.mapping;
 463       if (ssm.fromSeq == seqDs || mapping.to == seqDs)
 464       {
 465         if (!related.contains(mapping.to))
 466         {
 467           result.add(mapping);
 468           related.add(mapping.to);
 469         }
 470       }
 471     }
 472     return result;
 473   }
 474
 475   /**
 476    * Test whether the given sequence is substitutable for one or more dummy
 477    * sequences in this mapping
 478    *
 479    * @param map
 480    * @param seq
 481    * @return
 482    */
 483   public boolean isRealisableWith(SequenceI seq)
 484   {
 485     return realiseWith(seq, false) > 0;
 486   }
 487
 488   /**
 489    * Replace any matchable mapped dummy sequences with the given real one.
 490    * Returns the count of sequence mappings instantiated.
 491    *
 492    * @param seq
 493    * @return
 494    */
 495   public int realiseWith(SequenceI seq)
 496   {
 497     return realiseWith(seq, true);
 498   }
 499
 500   /**
 501    * Returns the number of mapped dummy sequences that could be replaced with
 502    * the given real sequence.
 503    *
 504    * @param seq
 505    *          a dataset sequence
 506    * @param doUpdate
 507    *          if true, performs replacements, else only counts
 508    * @return
 509    */
 510   protected int realiseWith(SequenceI seq, boolean doUpdate)
 511   {
 512     SequenceI ds = seq.getDatasetSequence() != null ? seq
 513             .getDatasetSequence() : seq;
 514     int count = 0;
 515
 516     /*
 517      * check for replaceable DNA ('map from') sequences
 518      */
 519     for (SequenceToSequenceMapping ssm : mappings)
 520     {
 521       SequenceI dna = ssm.fromSeq;
 522       if (dna instanceof SequenceDummy
 523               && dna.getName().equals(ds.getName()))
 524       {
 525         Mapping mapping = ssm.mapping;
 526         int mapStart = mapping.getMap().getFromLowest();
 527         int mapEnd = mapping.getMap().getFromHighest();
 528         boolean mappable = couldRealiseSequence(dna, ds, mapStart, mapEnd);
 529         if (mappable)
 530         {
 531           count++;
 532           if (doUpdate)
 533           {
 534             // TODO: new method ? ds.realise(dna);
 535             // might want to copy database refs as well
 536             ds.setSequenceFeatures(dna.getSequenceFeatures());
 537             // dnaSeqs[i] = ds;
 538             ssm.fromSeq = ds;
 539             System.out.println("Realised mapped sequence " + ds.getName());
 540           }
 541         }
 542       }
 543
 544       /*
 545        * check for replaceable protein ('map to') sequences
 546        */
 547       Mapping mapping = ssm.mapping;
 548       SequenceI prot = mapping.getTo();
 549       int mapStart = mapping.getMap().getToLowest();
 550       int mapEnd = mapping.getMap().getToHighest();
 551       boolean mappable = couldRealiseSequence(prot, ds, mapStart, mapEnd);
 552       if (mappable)
 553       {
 554         count++;
 555         if (doUpdate)
 556         {
 557           // TODO: new method ? ds.realise(dna);
 558           // might want to copy database refs as well
 559           ds.setSequenceFeatures(dna.getSequenceFeatures());
 560           ssm.mapping.setTo(ds);
 561         }
 562       }
 563     }
 564     return count;
 565   }
 566
 567   /**
 568    * Helper method to test whether a 'real' sequence could replace a 'dummy'
 569    * sequence in the map. The criteria are that they have the same name, and
 570    * that the mapped region overlaps the candidate sequence.
 571    *
 572    * @param existing
 573    * @param replacement
 574    * @param mapStart
 575    * @param mapEnd
 576    * @return
 577    */
 578   protected static boolean couldRealiseSequence(SequenceI existing,
 579           SequenceI replacement, int mapStart, int mapEnd)
 580   {
 581     if (existing instanceof SequenceDummy
 582             && !(replacement instanceof SequenceDummy)
 583             && existing.getName().equals(replacement.getName()))
 584     {
 585       int start = replacement.getStart();
 586       int end = replacement.getEnd();
 587       boolean mappingOverlapsSequence = (mapStart >= start && mapStart <= end)
 588               || (mapEnd >= start && mapEnd <= end);
 589       if (mappingOverlapsSequence)
 590       {
 591         return true;
 592       }
 593     }
 594     return false;
 595   }
 596
 597   /**
 598    * Change any mapping to the given sequence to be to its dataset sequence
 599    * instead. For use when mappings are created before their referenced
 600    * sequences are instantiated, for example when parsing GFF data.
 601    *
 602    * @param seq
 603    */
 604   public void updateToDataset(SequenceI seq)
 605   {
 606     if (seq == null || seq.getDatasetSequence() == null)
 607     {
 608       return;
 609     }
 610     SequenceI ds = seq.getDatasetSequence();
 611
 612     for (SequenceToSequenceMapping ssm : mappings)
 613     /*
 614      * 'from' sequences
 615      */
 616     {
 617       if (ssm.fromSeq == seq)
 618       {
 619         ssm.fromSeq = ds;
 620       }
 621
 622       /*
 623        * 'to' sequences
 624        */
 625       if (ssm.mapping.to == seq)
 626       {
 627         ssm.mapping.to = ds;
 628       }
 629     }
 630   }
 631
 632   /**
 633    * Answers true if this object contains no mappings
 634    *
 635    * @return
 636    */
 637   public boolean isEmpty()
 638   {
 639     return mappings.isEmpty();
 640   }
 641
 642   /**
 643    * Method for debug / inspection purposes only, may change in future
 644    */
 645   @Override
 646   public String toString()
 647   {
 648     return mappings == null ? "null" : mappings.toString();
 649   }
 650 }