src/jalview/analysis/AlignmentUtils.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.analysis;
  22
  23 import jalview.datamodel.AlignedCodon;
  24 import jalview.datamodel.AlignedCodonFrame;
  25 import jalview.datamodel.AlignmentAnnotation;
  26 import jalview.datamodel.AlignmentI;
  27 import jalview.datamodel.Mapping;
  28 import jalview.datamodel.SearchResults;
  29 import jalview.datamodel.Sequence;
  30 import jalview.datamodel.SequenceI;
  31 import jalview.schemes.ResidueProperties;
  32 import jalview.util.MapList;
  33
  34 import java.util.ArrayList;
  35 import java.util.Arrays;
  36 import java.util.HashMap;
  37 import java.util.Iterator;
  38 import java.util.LinkedHashMap;
  39 import java.util.List;
  40 import java.util.Map;
  41 import java.util.Map.Entry;
  42 import java.util.Set;
  43 import java.util.TreeMap;
  44
  45 /**
  46  * grab bag of useful alignment manipulation operations Expect these to be
  47  * refactored elsewhere at some point.
  48  *
  49  * @author jimp
  50  *
  51  */
  52 public class AlignmentUtils
  53 {
  54
  55   /**
  56    * Represents the 3 possible results of trying to map one alignment to
  57    * another.
  58    */
  59   public enum MappingResult
  60   {
  61     Mapped, NotMapped, AlreadyMapped
  62   }
  63
  64   /**
  65    * given an existing alignment, create a new alignment including all, or up to
  66    * flankSize additional symbols from each sequence's dataset sequence
  67    *
  68    * @param core
  69    * @param flankSize
  70    * @return AlignmentI
  71    */
  72   public static AlignmentI expandContext(AlignmentI core, int flankSize)
  73   {
  74     List<SequenceI> sq = new ArrayList<SequenceI>();
  75     int maxoffset = 0;
  76     for (SequenceI s : core.getSequences())
  77     {
  78       SequenceI newSeq = s.deriveSequence();
  79       if (newSeq.getStart() > maxoffset
  80               && newSeq.getDatasetSequence().getStart() < s.getStart())
  81       {
  82         maxoffset = newSeq.getStart();
  83       }
  84       sq.add(newSeq);
  85     }
  86     if (flankSize > -1)
  87     {
  88       maxoffset = flankSize;
  89     }
  90     // now add offset to create a new expanded alignment
  91     for (SequenceI s : sq)
  92     {
  93       SequenceI ds = s;
  94       while (ds.getDatasetSequence() != null)
  95       {
  96         ds = ds.getDatasetSequence();
  97       }
  98       int s_end = s.findPosition(s.getStart() + s.getLength());
  99       // find available flanking residues for sequence
 100       int ustream_ds = s.getStart() - ds.getStart(), dstream_ds = ds
 101               .getEnd() - s_end;
 102
 103       // build new flanked sequence
 104
 105       // compute gap padding to start of flanking sequence
 106       int offset = maxoffset - ustream_ds;
 107
 108       // padding is gapChar x ( maxoffset - min(ustream_ds, flank)
 109       if (flankSize >= 0)
 110       {
 111         if (flankSize < ustream_ds)
 112         {
 113           // take up to flankSize residues
 114           offset = maxoffset - flankSize;
 115           ustream_ds = flankSize;
 116         }
 117         if (flankSize < dstream_ds)
 118         {
 119           dstream_ds = flankSize;
 120         }
 121       }
 122       char[] upstream = new String(ds.getSequence(s.getStart() - 1
 123               - ustream_ds, s.getStart() - 1)).toLowerCase().toCharArray();
 124       char[] downstream = new String(ds.getSequence(s_end - 1, s_end + 1
 125               + dstream_ds)).toLowerCase().toCharArray();
 126       char[] coreseq = s.getSequence();
 127       char[] nseq = new char[offset + upstream.length + downstream.length
 128               + coreseq.length];
 129       char c = core.getGapCharacter();
 130       // TODO could lowercase the flanking regions
 131       int p = 0;
 132       for (; p < offset; p++)
 133       {
 134         nseq[p] = c;
 135       }
 136       // s.setSequence(new String(upstream).toLowerCase()+new String(coreseq) +
 137       // new String(downstream).toLowerCase());
 138       System.arraycopy(upstream, 0, nseq, p, upstream.length);
 139       System.arraycopy(coreseq, 0, nseq, p + upstream.length,
 140               coreseq.length);
 141       System.arraycopy(downstream, 0, nseq, p + coreseq.length
 142               + upstream.length, downstream.length);
 143       s.setSequence(new String(nseq));
 144       s.setStart(s.getStart() - ustream_ds);
 145       s.setEnd(s_end + downstream.length);
 146     }
 147     AlignmentI newAl = new jalview.datamodel.Alignment(
 148             sq.toArray(new SequenceI[0]));
 149     for (SequenceI s : sq)
 150     {
 151       if (s.getAnnotation() != null)
 152       {
 153         for (AlignmentAnnotation aa : s.getAnnotation())
 154         {
 155           newAl.addAnnotation(aa);
 156         }
 157       }
 158     }
 159     newAl.setDataset(core.getDataset());
 160     return newAl;
 161   }
 162
 163   /**
 164    * Returns the index (zero-based position) of a sequence in an alignment, or
 165    * -1 if not found.
 166    *
 167    * @param al
 168    * @param seq
 169    * @return
 170    */
 171   public static int getSequenceIndex(AlignmentI al, SequenceI seq)
 172   {
 173     int result = -1;
 174     int pos = 0;
 175     for (SequenceI alSeq : al.getSequences())
 176     {
 177       if (alSeq == seq)
 178       {
 179         result = pos;
 180         break;
 181       }
 182       pos++;
 183     }
 184     return result;
 185   }
 186
 187   /**
 188    * Returns a map of lists of sequences in the alignment, keyed by sequence
 189    * name. For use in mapping between different alignment views of the same
 190    * sequences.
 191    *
 192    * @see jalview.datamodel.AlignmentI#getSequencesByName()
 193    */
 194   public static Map<String, List<SequenceI>> getSequencesByName(
 195           AlignmentI al)
 196   {
 197     Map<String, List<SequenceI>> theMap = new LinkedHashMap<String, List<SequenceI>>();
 198     for (SequenceI seq : al.getSequences())
 199     {
 200       String name = seq.getName();
 201       if (name != null)
 202       {
 203         List<SequenceI> seqs = theMap.get(name);
 204         if (seqs == null)
 205         {
 206           seqs = new ArrayList<SequenceI>();
 207           theMap.put(name, seqs);
 208         }
 209         seqs.add(seq);
 210       }
 211     }
 212     return theMap;
 213   }
 214
 215   /**
 216    * Build mapping of protein to cDNA alignment. Mappings are made between
 217    * sequences where the cDNA translates to the protein sequence. Any new
 218    * mappings are added to the protein alignment. Has a 3-valued result: either
 219    * Mapped (at least one sequence mapping was created), AlreadyMapped (all
 220    * possible sequence mappings already exist), or NotMapped (no possible
 221    * sequence mappings exist).
 222    *
 223    * @param proteinAlignment
 224    * @param cdnaAlignment
 225    * @return
 226    */
 227   public static MappingResult mapProteinToCdna(
 228           final AlignmentI proteinAlignment,
 229           final AlignmentI cdnaAlignment)
 230   {
 231     if (proteinAlignment == null || cdnaAlignment == null)
 232     {
 233       return MappingResult.NotMapped;
 234     }
 235
 236     boolean mappingPossible = false;
 237     boolean mappingPerformed = false;
 238
 239     List<SequenceI> mapped = new ArrayList<SequenceI>();
 240
 241     List<SequenceI> thisSeqs = proteinAlignment.getSequences();
 242
 243     for (SequenceI aaSeq : thisSeqs)
 244     {
 245       AlignedCodonFrame acf = new AlignedCodonFrame();
 246
 247       for (SequenceI cdnaSeq : cdnaAlignment.getSequences())
 248       {
 249         /*
 250          * Heuristic rule: don't map more than one AA sequence to the same cDNA;
 251          * map progressively assuming that alignments have mappable sequences in
 252          * the same respective order
 253          */
 254         if (mapped.contains(cdnaSeq))
 255         {
 256           continue;
 257         }
 258         if (!mappingExists(proteinAlignment.getCodonFrames(),
 259                 aaSeq.getDatasetSequence(), cdnaSeq.getDatasetSequence()))
 260         {
 261           MapList map = mapProteinToCdna(aaSeq, cdnaSeq);
 262           if (map != null)
 263           {
 264             acf.addMap(cdnaSeq, aaSeq, map);
 265             mappingPerformed = true;
 266             mapped.add(cdnaSeq);
 267
 268             /*
 269              * Heuristic rule #2: don't map AA sequence to more than one cDNA
 270              */
 271             break;
 272           }
 273         }
 274       }
 275       proteinAlignment.addCodonFrame(acf);
 276     }
 277
 278     /*
 279      * If at least one mapping was possible but none was done, then the
 280      * alignments are already as mapped as they can be.
 281      */
 282     if (mappingPossible && !mappingPerformed)
 283     {
 284       return MappingResult.AlreadyMapped;
 285     }
 286     else
 287     {
 288       return mappingPerformed ? MappingResult.Mapped
 289               : MappingResult.NotMapped;
 290     }
 291   }
 292
 293   /**
 294    * Answers true if the mappings include one between the given (dataset)
 295    * sequences.
 296    */
 297   public static boolean mappingExists(Set<AlignedCodonFrame> set,
 298           SequenceI aaSeq, SequenceI cdnaSeq)
 299   {
 300     if (set != null)
 301     {
 302       for (AlignedCodonFrame acf : set)
 303       {
 304         if (cdnaSeq == acf.getDnaForAaSeq(aaSeq))
 305         {
 306           return true;
 307         }
 308       }
 309     }
 310     return false;
 311   }
 312
 313   /**
 314    * Build a mapping (if possible) of a protein to a cDNA sequence. The cDNA
 315    * must be three times the length of the protein, possibly after ignoring
 316    * start and/or stop codons, and must translate to the protein. Returns null
 317    * if no mapping is determined.
 318    *
 319    * @param proteinSeqs
 320    * @param cdnaSeq
 321    * @return
 322    */
 323   public static MapList mapProteinToCdna(SequenceI proteinSeq,
 324           SequenceI cdnaSeq)
 325   {
 326     /*
 327      * Here we handle either dataset sequence set (desktop) or absent (applet).
 328      * Use only the char[] form of the sequence to avoid creating possibly large
 329      * String objects.
 330      */
 331     final SequenceI proteinDataset = proteinSeq.getDatasetSequence();
 332     char[] aaSeqChars = proteinDataset != null ? proteinDataset
 333             .getSequence() : proteinSeq.getSequence();
 334     final SequenceI cdnaDataset = cdnaSeq.getDatasetSequence();
 335     char[] cdnaSeqChars = cdnaDataset != null ? cdnaDataset.getSequence()
 336             : cdnaSeq.getSequence();
 337     if (aaSeqChars == null || cdnaSeqChars == null)
 338     {
 339       return null;
 340     }
 341
 342     /*
 343      * cdnaStart/End, proteinStartEnd are base 1 (for dataset sequence mapping)
 344      */
 345     final int mappedLength = 3 * aaSeqChars.length;
 346     int cdnaLength = cdnaSeqChars.length;
 347     int cdnaStart = 1;
 348     int cdnaEnd = cdnaLength;
 349     final int proteinStart = 1;
 350     final int proteinEnd = aaSeqChars.length;
 351
 352     /*
 353      * If lengths don't match, try ignoring stop codon.
 354      */
 355     if (cdnaLength != mappedLength && cdnaLength > 2)
 356     {
 357       String lastCodon = String.valueOf(cdnaSeqChars, cdnaLength - 3, 3)
 358               .toUpperCase();
 359       for (String stop : ResidueProperties.STOP)
 360       {
 361         if (lastCodon.equals(stop))
 362         {
 363           cdnaEnd -= 3;
 364           cdnaLength -= 3;
 365           break;
 366         }
 367       }
 368     }
 369
 370     /*
 371      * If lengths still don't match, try ignoring start codon.
 372      */
 373     if (cdnaLength != mappedLength
 374             && cdnaLength > 2
 375             && String.valueOf(cdnaSeqChars, 0, 3).toUpperCase()
 376                     .equals(
 377                     ResidueProperties.START))
 378     {
 379       cdnaStart += 3;
 380       cdnaLength -= 3;
 381     }
 382
 383     if (cdnaLength != mappedLength)
 384     {
 385       return null;
 386     }
 387     if (!translatesAs(cdnaSeqChars, cdnaStart - 1, aaSeqChars))
 388     {
 389       return null;
 390     }
 391     MapList map = new MapList(new int[]
 392     { cdnaStart, cdnaEnd }, new int[]
 393     { proteinStart, proteinEnd }, 3, 1);
 394     return map;
 395   }
 396
 397   /**
 398    * Test whether the given cdna sequence, starting at the given offset,
 399    * translates to the given amino acid sequence, using the standard translation
 400    * table. Designed to fail fast i.e. as soon as a mismatch position is found.
 401    *
 402    * @param cdnaSeqChars
 403    * @param cdnaStart
 404    * @param aaSeqChars
 405    * @return
 406    */
 407   protected static boolean translatesAs(char[] cdnaSeqChars, int cdnaStart,
 408           char[] aaSeqChars)
 409   {
 410     int aaResidue = 0;
 411     for (int i = cdnaStart; i < cdnaSeqChars.length - 2
 412             && aaResidue < aaSeqChars.length; i += 3)
 413     {
 414       String codon = String.valueOf(cdnaSeqChars, i, 3);
 415       final String translated = ResidueProperties.codonTranslate(
 416               codon);
 417       if (translated == null
 418               || !(aaSeqChars[aaResidue] == translated.charAt(0)))
 419       {
 420         return false;
 421       }
 422       aaResidue++;
 423     }
 424     // fail if we didn't match all of the aa sequence
 425     return (aaResidue == aaSeqChars.length);
 426   }
 427
 428   /**
 429    * Align sequence 'seq' to match the alignment of a mapped sequence. Note this
 430    * currently assumes that we are aligning cDNA to match protein.
 431    *
 432    * @param seq
 433    *          the sequence to be realigned
 434    * @param al
 435    *          the alignment whose sequence alignment is to be 'copied'
 436    * @param gap
 437    *          character string represent a gap in the realigned sequence
 438    * @param preserveUnmappedGaps
 439    * @param preserveMappedGaps
 440    * @return true if the sequence was realigned, false if it could not be
 441    */
 442   public static boolean alignSequenceAs(SequenceI seq, AlignmentI al,
 443           String gap, boolean preserveMappedGaps,
 444           boolean preserveUnmappedGaps)
 445   {
 446     /*
 447      * Get any mappings from the source alignment to the target (dataset) sequence.
 448      */
 449     // TODO there may be one AlignedCodonFrame per dataset sequence, or one with
 450     // all mappings. Would it help to constrain this?
 451     List<AlignedCodonFrame> mappings = al.getCodonFrame(seq);
 452     if (mappings == null || mappings.isEmpty())
 453     {
 454       return false;
 455     }
 456
 457     /*
 458      * Locate the aligned source sequence whose dataset sequence is mapped. We
 459      * just take the first match here (as we can't align cDNA like more than one
 460      * protein sequence).
 461      */
 462     SequenceI alignFrom = null;
 463     AlignedCodonFrame mapping = null;
 464     for (AlignedCodonFrame mp : mappings)
 465     {
 466       alignFrom = mp.findAlignedSequence(seq.getDatasetSequence(), al);
 467       if (alignFrom != null)
 468       {
 469         mapping = mp;
 470         break;
 471       }
 472     }
 473
 474     if (alignFrom == null)
 475     {
 476       return false;
 477     }
 478     alignSequenceAs(seq, alignFrom, mapping, gap, al.getGapCharacter(),
 479             preserveMappedGaps, preserveUnmappedGaps);
 480     return true;
 481   }
 482
 483   /**
 484    * Align sequence 'alignTo' the same way as 'alignFrom', using the mapping to
 485    * match residues and codons. Flags control whether existing gaps in unmapped
 486    * (intron) and mapped (exon) regions are preserved or not. Gaps linking intro
 487    * and exon are only retained if both flags are set.
 488    *
 489    * @param alignTo
 490    * @param alignFrom
 491    * @param mapping
 492    * @param myGap
 493    * @param sourceGap
 494    * @param preserveUnmappedGaps
 495    * @param preserveMappedGaps
 496    */
 497   public static void alignSequenceAs(SequenceI alignTo,
 498           SequenceI alignFrom,
 499           AlignedCodonFrame mapping, String myGap, char sourceGap,
 500           boolean preserveMappedGaps, boolean preserveUnmappedGaps)
 501   {
 502     // TODO generalise to work for Protein-Protein, dna-dna, dna-protein
 503     final char[] thisSeq = alignTo.getSequence();
 504     final char[] thatAligned = alignFrom.getSequence();
 505     StringBuilder thisAligned = new StringBuilder(2 * thisSeq.length);
 506
 507     // aligned and dataset sequence positions, all base zero
 508     int thisSeqPos = 0;
 509     int sourceDsPos = 0;
 510
 511     int basesWritten = 0;
 512     char myGapChar = myGap.charAt(0);
 513     int ratio = myGap.length();
 514
 515     /*
 516      * Traverse the aligned protein sequence.
 517      */
 518     int sourceGapMappedLength = 0;
 519     boolean inExon = false;
 520     for (char sourceChar : thatAligned)
 521     {
 522       if (sourceChar == sourceGap)
 523       {
 524         sourceGapMappedLength += ratio;
 525         continue;
 526       }
 527
 528       /*
 529        * Found a residue. Locate its mapped codon (start) position.
 530        */
 531       sourceDsPos++;
 532       // Note mapping positions are base 1, our sequence positions base 0
 533       int[] mappedPos = mapping.getMappedRegion(alignTo, alignFrom,
 534               sourceDsPos);
 535       if (mappedPos == null)
 536       {
 537         /*
 538          * Abort realignment if unmapped protein. Or could ignore it??
 539          */
 540         System.err.println("Can't align: no codon mapping to residue "
 541                 + sourceDsPos + "(" + sourceChar + ")");
 542         return;
 543       }
 544
 545       int mappedCodonStart = mappedPos[0]; // position (1...) of codon start
 546       int mappedCodonEnd = mappedPos[mappedPos.length - 1]; // codon end pos
 547       StringBuilder trailingCopiedGap = new StringBuilder();
 548
 549       /*
 550        * Copy dna sequence up to and including this codon. Optionally, include
 551        * gaps before the codon starts (in introns) and/or after the codon starts
 552        * (in exons).
 553        *
 554        * Note this only works for 'linear' splicing, not reverse or interleaved.
 555        * But then 'align dna as protein' doesn't make much sense otherwise.
 556        */
 557       int intronLength = 0;
 558       while (basesWritten < mappedCodonEnd && thisSeqPos < thisSeq.length)
 559       {
 560         final char c = thisSeq[thisSeqPos++];
 561         if (c != myGapChar)
 562         {
 563           basesWritten++;
 564
 565           if (basesWritten < mappedCodonStart)
 566           {
 567             /*
 568              * Found an unmapped (intron) base. First add in any preceding gaps
 569              * (if wanted).
 570              */
 571             if (preserveUnmappedGaps && trailingCopiedGap.length() > 0)
 572             {
 573               thisAligned.append(trailingCopiedGap.toString());
 574               intronLength += trailingCopiedGap.length();
 575               trailingCopiedGap = new StringBuilder();
 576             }
 577             intronLength++;
 578             inExon = false;
 579           }
 580           else
 581           {
 582             final boolean startOfCodon = basesWritten == mappedCodonStart;
 583             int gapsToAdd = calculateGapsToInsert(preserveMappedGaps,
 584                     preserveUnmappedGaps, sourceGapMappedLength, inExon,
 585                     trailingCopiedGap.length(), intronLength, startOfCodon);
 586             for (int i = 0; i < gapsToAdd; i++)
 587             {
 588               thisAligned.append(myGapChar);
 589             }
 590             sourceGapMappedLength = 0;
 591             inExon = true;
 592           }
 593           thisAligned.append(c);
 594           trailingCopiedGap = new StringBuilder();
 595         }
 596         else
 597         {
 598           if (inExon && preserveMappedGaps)
 599           {
 600             trailingCopiedGap.append(myGapChar);
 601           }
 602           else if (!inExon && preserveUnmappedGaps)
 603           {
 604             trailingCopiedGap.append(myGapChar);
 605           }
 606         }
 607       }
 608     }
 609
 610     /*
 611      * At end of protein sequence. Copy any remaining dna sequence, optionally
 612      * including (intron) gaps. We do not copy trailing gaps in protein.
 613      */
 614     while (thisSeqPos < thisSeq.length)
 615     {
 616       final char c = thisSeq[thisSeqPos++];
 617       if (c != myGapChar || preserveUnmappedGaps)
 618       {
 619         thisAligned.append(c);
 620       }
 621     }
 622
 623     /*
 624      * All done aligning, set the aligned sequence.
 625      */
 626     alignTo.setSequence(new String(thisAligned));
 627   }
 628
 629   /**
 630    * Helper method to work out how many gaps to insert when realigning.
 631    *
 632    * @param preserveMappedGaps
 633    * @param preserveUnmappedGaps
 634    * @param sourceGapMappedLength
 635    * @param inExon
 636    * @param trailingCopiedGap
 637    * @param intronLength
 638    * @param startOfCodon
 639    * @return
 640    */
 641   protected static int calculateGapsToInsert(boolean preserveMappedGaps,
 642           boolean preserveUnmappedGaps, int sourceGapMappedLength,
 643           boolean inExon, int trailingGapLength,
 644           int intronLength, final boolean startOfCodon)
 645   {
 646     int gapsToAdd = 0;
 647     if (startOfCodon)
 648     {
 649       /*
 650        * Reached start of codon. Ignore trailing gaps in intron unless we are
 651        * preserving gaps in both exon and intron. Ignore them anyway if the
 652        * protein alignment introduces a gap at least as large as the intronic
 653        * region.
 654        */
 655       if (inExon && !preserveMappedGaps)
 656       {
 657         trailingGapLength = 0;
 658       }
 659       if (!inExon && !(preserveMappedGaps && preserveUnmappedGaps))
 660       {
 661         trailingGapLength = 0;
 662       }
 663       if (inExon)
 664       {
 665         gapsToAdd = Math.max(sourceGapMappedLength, trailingGapLength);
 666       }
 667       else
 668       {
 669         if (intronLength + trailingGapLength <= sourceGapMappedLength)
 670         {
 671           gapsToAdd = sourceGapMappedLength - intronLength;
 672         }
 673         else
 674         {
 675           gapsToAdd = Math.min(intronLength + trailingGapLength
 676                   - sourceGapMappedLength, trailingGapLength);
 677         }
 678       }
 679     }
 680     else
 681     {
 682       /*
 683        * second or third base of codon; check for any gaps in dna
 684        */
 685       if (!preserveMappedGaps)
 686       {
 687         trailingGapLength = 0;
 688       }
 689       gapsToAdd = Math.max(sourceGapMappedLength, trailingGapLength);
 690     }
 691     return gapsToAdd;
 692   }
 693
 694   /**
 695    * Returns a list of sequences mapped from the given sequences and aligned
 696    * (gapped) in the same way. For example, the cDNA for aligned protein, where
 697    * a single gap in protein generates three gaps in cDNA.
 698    *
 699    * @param sequences
 700    * @param gapCharacter
 701    * @param mappings
 702    * @return
 703    */
 704   public static List<SequenceI> getAlignedTranslation(
 705           List<SequenceI> sequences, char gapCharacter,
 706           Set<AlignedCodonFrame> mappings)
 707   {
 708     List<SequenceI> alignedSeqs = new ArrayList<SequenceI>();
 709
 710     for (SequenceI seq : sequences)
 711     {
 712       List<SequenceI> mapped = getAlignedTranslation(seq, gapCharacter,
 713               mappings);
 714       alignedSeqs.addAll(mapped);
 715     }
 716     return alignedSeqs;
 717   }
 718
 719   /**
 720    * Returns sequences aligned 'like' the source sequence, as mapped by the
 721    * given mappings. Normally we expect zero or one 'mapped' sequences, but this
 722    * will support 1-to-many as well.
 723    *
 724    * @param seq
 725    * @param gapCharacter
 726    * @param mappings
 727    * @return
 728    */
 729   protected static List<SequenceI> getAlignedTranslation(SequenceI seq,
 730           char gapCharacter, Set<AlignedCodonFrame> mappings)
 731   {
 732     List<SequenceI> result = new ArrayList<SequenceI>();
 733     for (AlignedCodonFrame mapping : mappings)
 734     {
 735       if (mapping.involvesSequence(seq))
 736       {
 737         SequenceI mapped = getAlignedTranslation(seq, gapCharacter, mapping);
 738         if (mapped != null)
 739         {
 740           result.add(mapped);
 741         }
 742       }
 743     }
 744     return result;
 745   }
 746
 747   /**
 748    * Returns the translation of 'seq' (as held in the mapping) with
 749    * corresponding alignment (gaps).
 750    *
 751    * @param seq
 752    * @param gapCharacter
 753    * @param mapping
 754    * @return
 755    */
 756   protected static SequenceI getAlignedTranslation(SequenceI seq,
 757           char gapCharacter, AlignedCodonFrame mapping)
 758   {
 759     String gap = String.valueOf(gapCharacter);
 760     boolean toDna = false;
 761     int fromRatio = 1;
 762     SequenceI mapTo = mapping.getDnaForAaSeq(seq);
 763     if (mapTo != null)
 764     {
 765       // mapping is from protein to nucleotide
 766       toDna = true;
 767       // should ideally get gap count ratio from mapping
 768       gap = String.valueOf(new char[]
 769       { gapCharacter, gapCharacter, gapCharacter });
 770     }
 771     else
 772     {
 773       // mapping is from nucleotide to protein
 774       mapTo = mapping.getAaForDnaSeq(seq);
 775       fromRatio = 3;
 776     }
 777     StringBuilder newseq = new StringBuilder(seq.getLength()
 778             * (toDna ? 3 : 1));
 779
 780     int residueNo = 0; // in seq, base 1
 781     int[] phrase = new int[fromRatio];
 782     int phraseOffset = 0;
 783     int gapWidth = 0;
 784     boolean first = true;
 785     final Sequence alignedSeq = new Sequence("", "");
 786
 787     for (char c : seq.getSequence())
 788     {
 789       if (c == gapCharacter)
 790       {
 791         gapWidth++;
 792         if (gapWidth >= fromRatio)
 793         {
 794           newseq.append(gap);
 795           gapWidth = 0;
 796         }
 797       }
 798       else
 799       {
 800         phrase[phraseOffset++] = residueNo + 1;
 801         if (phraseOffset == fromRatio)
 802         {
 803           /*
 804            * Have read a whole codon (or protein residue), now translate: map
 805            * source phrase to positions in target sequence add characters at
 806            * these positions to newseq Note mapping positions are base 1, our
 807            * sequence positions base 0.
 808            */
 809           SearchResults sr = new SearchResults();
 810           for (int pos : phrase)
 811           {
 812             mapping.markMappedRegion(seq, pos, sr);
 813           }
 814           newseq.append(sr.toString());
 815           if (first)
 816           {
 817             first = false;
 818             // Hack: Copy sequence dataset, name and description from
 819             // SearchResults.match[0].sequence
 820             // TODO? carry over sequence names from original 'complement'
 821             // alignment
 822             SequenceI mappedTo = sr.getResultSequence(0);
 823             alignedSeq.setName(mappedTo.getName());
 824             alignedSeq.setDescription(mappedTo.getDescription());
 825             alignedSeq.setDatasetSequence(mappedTo);
 826           }
 827           phraseOffset = 0;
 828         }
 829         residueNo++;
 830       }
 831     }
 832     alignedSeq.setSequence(newseq.toString());
 833     return alignedSeq;
 834   }
 835
 836   /**
 837    * Realigns the given protein to match the alignment of the dna, using codon
 838    * mappings to translate aligned codon positions to protein residues.
 839    *
 840    * @param protein
 841    *          the alignment whose sequences are realigned by this method
 842    * @param dna
 843    *          the dna alignment whose alignment we are 'copying'
 844    * @return the number of sequences that were realigned
 845    */
 846   public static int alignProteinAsDna(AlignmentI protein, AlignmentI dna)
 847   {
 848     Set<AlignedCodonFrame> mappings = protein.getCodonFrames();
 849
 850     /*
 851      * Map will hold, for each aligned codon position e.g. [3, 5, 6], a map of
 852      * {dnaSequence, {proteinSequence, codonProduct}} at that position. The
 853      * comparator keeps the codon positions ordered.
 854      */
 855     Map<AlignedCodon, Map<SequenceI, String>> alignedCodons = new TreeMap<AlignedCodon, Map<SequenceI, String>>(
 856             new CodonComparator());
 857     for (SequenceI dnaSeq : dna.getSequences())
 858     {
 859       for (AlignedCodonFrame mapping : mappings)
 860       {
 861         Mapping seqMap = mapping.getMappingForSequence(dnaSeq);
 862         SequenceI prot = mapping.findAlignedSequence(
 863                 dnaSeq.getDatasetSequence(), protein);
 864         if (prot != null)
 865         {
 866           addCodonPositions(dnaSeq, prot, protein.getGapCharacter(),
 867                   seqMap, alignedCodons);
 868         }
 869       }
 870     }
 871     return alignProteinAs(protein, alignedCodons);
 872   }
 873
 874   /**
 875    * Update the aligned protein sequences to match the codon alignments given in
 876    * the map.
 877    *
 878    * @param protein
 879    * @param alignedCodons
 880    *          an ordered map of codon positions (columns), with sequence/peptide
 881    *          values present in each column
 882    * @return
 883    */
 884   protected static int alignProteinAs(AlignmentI protein,
 885           Map<AlignedCodon, Map<SequenceI, String>> alignedCodons)
 886   {
 887     /*
 888      * Prefill aligned sequences with gaps before inserting aligned protein
 889      * residues.
 890      */
 891     int alignedWidth = alignedCodons.size();
 892     char[] gaps = new char[alignedWidth];
 893     Arrays.fill(gaps, protein.getGapCharacter());
 894     String allGaps = String.valueOf(gaps);
 895     for (SequenceI seq : protein.getSequences())
 896     {
 897       seq.setSequence(allGaps);
 898     }
 899
 900     int column = 0;
 901     for (AlignedCodon codon : alignedCodons.keySet())
 902     {
 903       final Map<SequenceI, String> columnResidues = alignedCodons.get(codon);
 904       for (Entry<SequenceI, String> entry : columnResidues
 905               .entrySet())
 906       {
 907         // place translated codon at its column position in sequence
 908         entry.getKey().getSequence()[column] = entry.getValue().charAt(0);
 909       }
 910       column++;
 911     }
 912     return 0;
 913   }
 914
 915   /**
 916    * Populate the map of aligned codons by traversing the given sequence
 917    * mapping, locating the aligned positions of mapped codons, and adding those
 918    * positions and their translation products to the map.
 919    *
 920    * @param dna
 921    *          the aligned sequence we are mapping from
 922    * @param protein
 923    *          the sequence to be aligned to the codons
 924    * @param gapChar
 925    *          the gap character in the dna sequence
 926    * @param seqMap
 927    *          a mapping to a sequence translation
 928    * @param alignedCodons
 929    *          the map we are building up
 930    */
 931   static void addCodonPositions(SequenceI dna, SequenceI protein,
 932           char gapChar,
 933           Mapping seqMap,
 934           Map<AlignedCodon, Map<SequenceI, String>> alignedCodons)
 935   {
 936     Iterator<AlignedCodon> codons = seqMap.getCodonIterator(dna, gapChar);
 937     while (codons.hasNext())
 938     {
 939       AlignedCodon codon = codons.next();
 940       Map<SequenceI, String> seqProduct = alignedCodons.get(codon);
 941       if (seqProduct == null)
 942       {
 943         seqProduct = new HashMap<SequenceI, String>();
 944         alignedCodons.put(codon, seqProduct);
 945       }
 946       seqProduct.put(protein, codon.product);
 947     }
 948   }
 949
 950   /**
 951    * Returns true if a cDNA/Protein mapping either exists, or could be made,
 952    * between at least one pair of sequences in the two alignments. Currently,
 953    * the logic is:
 954    * <ul>
 955    * <li>One alignment must be nucleotide, and the other protein</li>
 956    * <li>At least one pair of sequences must be already mapped, or mappable</li>
 957    * <li>Mappable means the nucleotide translation matches the protein sequence</li>
 958    * <li>The translation may ignore start and stop codons if present in the
 959    * nucleotide</li>
 960    * </ul>
 961    *
 962    * @param al1
 963    * @param al2
 964    * @return
 965    */
 966   public static boolean isMappable(AlignmentI al1, AlignmentI al2)
 967   {
 968     /*
 969      * Require one nucleotide and one protein
 970      */
 971     if (al1.isNucleotide() == al2.isNucleotide())
 972     {
 973       return false;
 974     }
 975     AlignmentI dna = al1.isNucleotide() ? al1 : al2;
 976     AlignmentI protein = dna == al1 ? al2 : al1;
 977     Set<AlignedCodonFrame> mappings = protein.getCodonFrames();
 978     for (SequenceI dnaSeq : dna.getSequences())
 979     {
 980       for (SequenceI proteinSeq : protein.getSequences())
 981       {
 982         if (isMappable(dnaSeq, proteinSeq, mappings))
 983         {
 984           return true;
 985         }
 986       }
 987     }
 988     return false;
 989   }
 990
 991   /**
 992    * Returns true if the dna sequence is mapped, or could be mapped, to the
 993    * protein sequence.
 994    *
 995    * @param dnaSeq
 996    * @param proteinSeq
 997    * @param mappings
 998    * @return
 999    */
1000   public static boolean isMappable(SequenceI dnaSeq, SequenceI proteinSeq,
1001           Set<AlignedCodonFrame> mappings)
1002   {
1003     SequenceI dnaDs = dnaSeq.getDatasetSequence() == null ? dnaSeq : dnaSeq.getDatasetSequence();
1004     SequenceI proteinDs = proteinSeq.getDatasetSequence() == null ? proteinSeq
1005             : proteinSeq.getDatasetSequence();
1006
1007     /*
1008      * Already mapped?
1009      */
1010     for (AlignedCodonFrame mapping : mappings) {
1011       if ( proteinDs == mapping.getAaForDnaSeq(dnaDs)) {
1012         return true;
1013       }
1014     }
1015
1016     /*
1017      * Just try to make a mapping (it is not yet stored), test whether
1018      * successful.
1019      */
1020     return mapProteinToCdna(proteinDs, dnaDs) != null;
1021   }
1022 }