src/jalview/analysis/AlignmentUtils.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer (Version 2.8.2)
   3  * Copyright (C) 2014 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 which have the same name and compatible lengths. Any new mappings
 218    * are added to the protein alignment. Has a 3-valued result: either Mapped
 219    * (at least one sequence mapping was created), AlreadyMapped (all possible
 220    * sequence mappings already exist), or NotMapped (no possible sequence
 221    * 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> thisSeqs = proteinAlignment.getSequences();
 240
 241     /*
 242      * Build a look-up of cDNA sequences by name, for matching purposes.
 243      */
 244     Map<String, List<SequenceI>> cdnaSeqs = cdnaAlignment
 245             .getSequencesByName();
 246
 247     for (SequenceI aaSeq : thisSeqs)
 248     {
 249       AlignedCodonFrame acf = new AlignedCodonFrame();
 250       List<SequenceI> candidates = cdnaSeqs.get(aaSeq.getName());
 251       if (candidates == null)
 252       {
 253         /*
 254          * No cDNA sequence with matching name, so no mapping possible for this
 255          * protein sequence
 256          */
 257         continue;
 258       }
 259       mappingPossible = true;
 260       for (SequenceI cdnaSeq : candidates)
 261       {
 262         if (!mappingExists(proteinAlignment.getCodonFrames(),
 263                 aaSeq.getDatasetSequence(), cdnaSeq.getDatasetSequence()))
 264         {
 265           MapList map = mapProteinToCdna(aaSeq, cdnaSeq);
 266           if (map != null)
 267           {
 268             acf.addMap(cdnaSeq, aaSeq, map);
 269             mappingPerformed = true;
 270           }
 271         }
 272       }
 273       proteinAlignment.addCodonFrame(acf);
 274     }
 275
 276     /*
 277      * If at least one mapping was possible but none was done, then the
 278      * alignments are already as mapped as they can be.
 279      */
 280     if (mappingPossible && !mappingPerformed)
 281     {
 282       return MappingResult.AlreadyMapped;
 283     }
 284     else
 285     {
 286       return mappingPerformed ? MappingResult.Mapped
 287               : MappingResult.NotMapped;
 288     }
 289   }
 290
 291   /**
 292    * Answers true if the mappings include one between the given (dataset)
 293    * sequences.
 294    */
 295   public static boolean mappingExists(Set<AlignedCodonFrame> set,
 296           SequenceI aaSeq, SequenceI cdnaSeq)
 297   {
 298     if (set != null)
 299     {
 300       for (AlignedCodonFrame acf : set)
 301       {
 302         if (cdnaSeq == acf.getDnaForAaSeq(aaSeq))
 303         {
 304           return true;
 305         }
 306       }
 307     }
 308     return false;
 309   }
 310
 311   /**
 312    * Build a mapping (if possible) of a protein to a cDNA sequence. The cDNA
 313    * must be three times the length of the protein, possibly after ignoring
 314    * start and/or stop codons. Returns null if no mapping is determined.
 315    *
 316    * @param proteinSeqs
 317    * @param cdnaSeq
 318    * @return
 319    */
 320   public static MapList mapProteinToCdna(SequenceI proteinSeq,
 321           SequenceI cdnaSeq)
 322   {
 323     /*
 324      * Here we handle either dataset sequence set (desktop) or absent (applet)
 325      */
 326     final SequenceI proteinDataset = proteinSeq.getDatasetSequence();
 327     String aaSeqString = proteinDataset != null ? proteinDataset
 328             .getSequenceAsString() : proteinSeq.getSequenceAsString();
 329     final SequenceI cdnaDataset = cdnaSeq.getDatasetSequence();
 330     String cdnaSeqString = cdnaDataset != null ? cdnaDataset
 331             .getSequenceAsString() : cdnaSeq.getSequenceAsString();
 332     if (aaSeqString == null || cdnaSeqString == null)
 333     {
 334       return null;
 335     }
 336
 337     final int mappedLength = 3 * aaSeqString.length();
 338     int cdnaLength = cdnaSeqString.length();
 339     int cdnaStart = 1;
 340     int cdnaEnd = cdnaLength;
 341     final int proteinStart = 1;
 342     final int proteinEnd = aaSeqString.length();
 343
 344     /*
 345      * If lengths don't match, try ignoring stop codon.
 346      */
 347     if (cdnaLength != mappedLength)
 348     {
 349       for (Object stop : ResidueProperties.STOP)
 350       {
 351         if (cdnaSeqString.toUpperCase().endsWith((String) stop))
 352         {
 353           cdnaEnd -= 3;
 354           cdnaLength -= 3;
 355           break;
 356         }
 357       }
 358     }
 359
 360     /*
 361      * If lengths still don't match, try ignoring start codon.
 362      */
 363     if (cdnaLength != mappedLength
 364             && cdnaSeqString.toUpperCase().startsWith(
 365                     ResidueProperties.START))
 366     {
 367       cdnaStart += 3;
 368       cdnaLength -= 3;
 369     }
 370
 371     if (cdnaLength == mappedLength)
 372     {
 373       MapList map = new MapList(new int[]
 374       { cdnaStart, cdnaEnd }, new int[]
 375       { proteinStart, proteinEnd }, 3, 1);
 376       return map;
 377     }
 378     else
 379     {
 380       return null;
 381     }
 382   }
 383
 384   /**
 385    * Align sequence 'seq' to match the alignment of a mapped sequence. Note this
 386    * currently assumes that we are aligning cDNA to match protein.
 387    *
 388    * @param seq
 389    *          the sequence to be realigned
 390    * @param al
 391    *          the alignment whose sequence alignment is to be 'copied'
 392    * @param gap
 393    *          character string represent a gap in the realigned sequence
 394    * @param preserveUnmappedGaps
 395    * @param preserveMappedGaps
 396    * @return true if the sequence was realigned, false if it could not be
 397    */
 398   public static boolean alignSequenceAs(SequenceI seq, AlignmentI al,
 399           String gap, boolean preserveMappedGaps,
 400           boolean preserveUnmappedGaps)
 401   {
 402     /*
 403      * Get any mappings from the source alignment to the target (dataset) sequence.
 404      */
 405     // TODO there may be one AlignedCodonFrame per dataset sequence, or one with
 406     // all mappings. Would it help to constrain this?
 407     List<AlignedCodonFrame> mappings = al.getCodonFrame(seq);
 408     if (mappings == null || mappings.isEmpty())
 409     {
 410       return false;
 411     }
 412
 413     /*
 414      * Locate the aligned source sequence whose dataset sequence is mapped. We
 415      * just take the first match here (as we can't align cDNA like more than one
 416      * protein sequence).
 417      */
 418     SequenceI alignFrom = null;
 419     AlignedCodonFrame mapping = null;
 420     for (AlignedCodonFrame mp : mappings)
 421     {
 422       alignFrom = mp.findAlignedSequence(seq.getDatasetSequence(), al);
 423       if (alignFrom != null)
 424       {
 425         mapping = mp;
 426         break;
 427       }
 428     }
 429
 430     if (alignFrom == null)
 431     {
 432       return false;
 433     }
 434     alignSequenceAs(seq, alignFrom, mapping, gap, al.getGapCharacter(),
 435             preserveMappedGaps, preserveUnmappedGaps);
 436     return true;
 437   }
 438
 439   /**
 440    * Align sequence 'alignTo' the same way as 'alignFrom', using the mapping to
 441    * match residues and codons. Flags control whether existing gaps in unmapped
 442    * (intron) and mapped (exon) regions are preserved or not. Gaps linking intro
 443    * and exon are only retained if both flags are set.
 444    *
 445    * @param alignTo
 446    * @param alignFrom
 447    * @param mapping
 448    * @param myGap
 449    * @param sourceGap
 450    * @param preserveUnmappedGaps
 451    * @param preserveMappedGaps
 452    */
 453   public static void alignSequenceAs(SequenceI alignTo,
 454           SequenceI alignFrom,
 455           AlignedCodonFrame mapping, String myGap, char sourceGap,
 456           boolean preserveMappedGaps, boolean preserveUnmappedGaps)
 457   {
 458     // TODO generalise to work for Protein-Protein, dna-dna, dna-protein
 459     final char[] thisSeq = alignTo.getSequence();
 460     final char[] thatAligned = alignFrom.getSequence();
 461     StringBuilder thisAligned = new StringBuilder(2 * thisSeq.length);
 462
 463     // aligned and dataset sequence positions, all base zero
 464     int thisSeqPos = 0;
 465     int sourceDsPos = 0;
 466
 467     int basesWritten = 0;
 468     char myGapChar = myGap.charAt(0);
 469     int ratio = myGap.length();
 470
 471     /*
 472      * Traverse the aligned protein sequence.
 473      */
 474     int sourceGapMappedLength = 0;
 475     boolean inExon = false;
 476     for (char sourceChar : thatAligned)
 477     {
 478       if (sourceChar == sourceGap)
 479       {
 480         sourceGapMappedLength += ratio;
 481         continue;
 482       }
 483
 484       /*
 485        * Found a residue. Locate its mapped codon (start) position.
 486        */
 487       sourceDsPos++;
 488       // Note mapping positions are base 1, our sequence positions base 0
 489       int[] mappedPos = mapping.getMappedRegion(alignTo, alignFrom,
 490               sourceDsPos);
 491       if (mappedPos == null)
 492       {
 493         /*
 494          * Abort realignment if unmapped protein. Or could ignore it??
 495          */
 496         System.err.println("Can't align: no codon mapping to residue "
 497                 + sourceDsPos + "(" + sourceChar + ")");
 498         return;
 499       }
 500
 501       int mappedCodonStart = mappedPos[0]; // position (1...) of codon start
 502       int mappedCodonEnd = mappedPos[mappedPos.length - 1]; // codon end pos
 503       StringBuilder trailingCopiedGap = new StringBuilder();
 504
 505       /*
 506        * Copy dna sequence up to and including this codon. Optionally, include
 507        * gaps before the codon starts (in introns) and/or after the codon starts
 508        * (in exons).
 509        *
 510        * Note this only works for 'linear' splicing, not reverse or interleaved.
 511        * But then 'align dna as protein' doesn't make much sense otherwise.
 512        */
 513       int intronLength = 0;
 514       while (basesWritten < mappedCodonEnd && thisSeqPos < thisSeq.length)
 515       {
 516         final char c = thisSeq[thisSeqPos++];
 517         if (c != myGapChar)
 518         {
 519           basesWritten++;
 520
 521           if (basesWritten < mappedCodonStart)
 522           {
 523             /*
 524              * Found an unmapped (intron) base. First add in any preceding gaps
 525              * (if wanted).
 526              */
 527             if (preserveUnmappedGaps && trailingCopiedGap.length() > 0)
 528             {
 529               thisAligned.append(trailingCopiedGap.toString());
 530               intronLength += trailingCopiedGap.length();
 531               trailingCopiedGap = new StringBuilder();
 532             }
 533             intronLength++;
 534             inExon = false;
 535           }
 536           else
 537           {
 538             final boolean startOfCodon = basesWritten == mappedCodonStart;
 539             int gapsToAdd = calculateGapsToInsert(preserveMappedGaps,
 540                     preserveUnmappedGaps, sourceGapMappedLength, inExon,
 541                     trailingCopiedGap.length(), intronLength, startOfCodon);
 542             for (int i = 0; i < gapsToAdd; i++)
 543             {
 544               thisAligned.append(myGapChar);
 545             }
 546             sourceGapMappedLength = 0;
 547             inExon = true;
 548           }
 549           thisAligned.append(c);
 550           trailingCopiedGap = new StringBuilder();
 551         }
 552         else
 553         {
 554           if (inExon && preserveMappedGaps)
 555           {
 556             trailingCopiedGap.append(myGapChar);
 557           }
 558           else if (!inExon && preserveUnmappedGaps)
 559           {
 560             trailingCopiedGap.append(myGapChar);
 561           }
 562         }
 563       }
 564     }
 565
 566     /*
 567      * At end of protein sequence. Copy any remaining dna sequence, optionally
 568      * including (intron) gaps. We do not copy trailing gaps in protein.
 569      */
 570     while (thisSeqPos < thisSeq.length)
 571     {
 572       final char c = thisSeq[thisSeqPos++];
 573       if (c != myGapChar || preserveUnmappedGaps)
 574       {
 575         thisAligned.append(c);
 576       }
 577     }
 578
 579     /*
 580      * All done aligning, set the aligned sequence.
 581      */
 582     alignTo.setSequence(new String(thisAligned));
 583   }
 584
 585   /**
 586    * Helper method to work out how many gaps to insert when realigning.
 587    *
 588    * @param preserveMappedGaps
 589    * @param preserveUnmappedGaps
 590    * @param sourceGapMappedLength
 591    * @param inExon
 592    * @param trailingCopiedGap
 593    * @param intronLength
 594    * @param startOfCodon
 595    * @return
 596    */
 597   protected static int calculateGapsToInsert(boolean preserveMappedGaps,
 598           boolean preserveUnmappedGaps, int sourceGapMappedLength,
 599           boolean inExon, int trailingGapLength,
 600           int intronLength, final boolean startOfCodon)
 601   {
 602     int gapsToAdd = 0;
 603     if (startOfCodon)
 604     {
 605       /*
 606        * Reached start of codon. Ignore trailing gaps in intron unless we are
 607        * preserving gaps in both exon and intron. Ignore them anyway if the
 608        * protein alignment introduces a gap at least as large as the intronic
 609        * region.
 610        */
 611       if (inExon && !preserveMappedGaps)
 612       {
 613         trailingGapLength = 0;
 614       }
 615       if (!inExon && !(preserveMappedGaps && preserveUnmappedGaps))
 616       {
 617         trailingGapLength = 0;
 618       }
 619       if (inExon)
 620       {
 621         gapsToAdd = Math.max(sourceGapMappedLength, trailingGapLength);
 622       }
 623       else
 624       {
 625         if (intronLength + trailingGapLength <= sourceGapMappedLength)
 626         {
 627           gapsToAdd = sourceGapMappedLength - intronLength;
 628         }
 629         else
 630         {
 631           gapsToAdd = Math.min(intronLength + trailingGapLength
 632                   - sourceGapMappedLength, trailingGapLength);
 633         }
 634       }
 635     }
 636     else
 637     {
 638       /*
 639        * second or third base of codon; check for any gaps in dna
 640        */
 641       if (!preserveMappedGaps)
 642       {
 643         trailingGapLength = 0;
 644       }
 645       gapsToAdd = Math.max(sourceGapMappedLength, trailingGapLength);
 646     }
 647     return gapsToAdd;
 648   }
 649
 650   /**
 651    * Returns a list of sequences mapped from the given sequences and aligned
 652    * (gapped) in the same way. For example, the cDNA for aligned protein, where
 653    * a single gap in protein generates three gaps in cDNA.
 654    *
 655    * @param sequences
 656    * @param gapCharacter
 657    * @param mappings
 658    * @return
 659    */
 660   public static List<SequenceI> getAlignedTranslation(
 661           List<SequenceI> sequences, char gapCharacter,
 662           Set<AlignedCodonFrame> mappings)
 663   {
 664     List<SequenceI> alignedSeqs = new ArrayList<SequenceI>();
 665
 666     for (SequenceI seq : sequences)
 667     {
 668       List<SequenceI> mapped = getAlignedTranslation(seq, gapCharacter,
 669               mappings);
 670       alignedSeqs.addAll(mapped);
 671     }
 672     return alignedSeqs;
 673   }
 674
 675   /**
 676    * Returns sequences aligned 'like' the source sequence, as mapped by the
 677    * given mappings. Normally we expect zero or one 'mapped' sequences, but this
 678    * will support 1-to-many as well.
 679    *
 680    * @param seq
 681    * @param gapCharacter
 682    * @param mappings
 683    * @return
 684    */
 685   protected static List<SequenceI> getAlignedTranslation(SequenceI seq,
 686           char gapCharacter, Set<AlignedCodonFrame> mappings)
 687   {
 688     List<SequenceI> result = new ArrayList<SequenceI>();
 689     for (AlignedCodonFrame mapping : mappings)
 690     {
 691       if (mapping.involvesSequence(seq))
 692       {
 693         SequenceI mapped = getAlignedTranslation(seq, gapCharacter, mapping);
 694         if (mapped != null)
 695         {
 696           result.add(mapped);
 697         }
 698       }
 699     }
 700     return result;
 701   }
 702
 703   /**
 704    * Returns the translation of 'seq' (as held in the mapping) with
 705    * corresponding alignment (gaps).
 706    *
 707    * @param seq
 708    * @param gapCharacter
 709    * @param mapping
 710    * @return
 711    */
 712   protected static SequenceI getAlignedTranslation(SequenceI seq,
 713           char gapCharacter, AlignedCodonFrame mapping)
 714   {
 715     String gap = String.valueOf(gapCharacter);
 716     boolean toDna = false;
 717     int fromRatio = 1;
 718     SequenceI mapTo = mapping.getDnaForAaSeq(seq);
 719     if (mapTo != null)
 720     {
 721       // mapping is from protein to nucleotide
 722       toDna = true;
 723       // should ideally get gap count ratio from mapping
 724       gap = String.valueOf(new char[]
 725       { gapCharacter, gapCharacter, gapCharacter });
 726     }
 727     else
 728     {
 729       // mapping is from nucleotide to protein
 730       mapTo = mapping.getAaForDnaSeq(seq);
 731       fromRatio = 3;
 732     }
 733     StringBuilder newseq = new StringBuilder(seq.getLength()
 734             * (toDna ? 3 : 1));
 735
 736     int residueNo = 0; // in seq, base 1
 737     int[] phrase = new int[fromRatio];
 738     int phraseOffset = 0;
 739     int gapWidth = 0;
 740     boolean first = true;
 741     final Sequence alignedSeq = new Sequence("", "");
 742
 743     for (char c : seq.getSequence())
 744     {
 745       if (c == gapCharacter)
 746       {
 747         gapWidth++;
 748         if (gapWidth >= fromRatio)
 749         {
 750           newseq.append(gap);
 751           gapWidth = 0;
 752         }
 753       }
 754       else
 755       {
 756         phrase[phraseOffset++] = residueNo + 1;
 757         if (phraseOffset == fromRatio)
 758         {
 759           /*
 760            * Have read a whole codon (or protein residue), now translate: map
 761            * source phrase to positions in target sequence add characters at
 762            * these positions to newseq Note mapping positions are base 1, our
 763            * sequence positions base 0.
 764            */
 765           SearchResults sr = new SearchResults();
 766           for (int pos : phrase)
 767           {
 768             mapping.markMappedRegion(seq, pos, sr);
 769           }
 770           newseq.append(sr.toString());
 771           if (first)
 772           {
 773             first = false;
 774             // Hack: Copy sequence dataset, name and description from
 775             // SearchResults.match[0].sequence
 776             // TODO? carry over sequence names from original 'complement'
 777             // alignment
 778             SequenceI mappedTo = sr.getResultSequence(0);
 779             alignedSeq.setName(mappedTo.getName());
 780             alignedSeq.setDescription(mappedTo.getDescription());
 781             alignedSeq.setDatasetSequence(mappedTo);
 782           }
 783           phraseOffset = 0;
 784         }
 785         residueNo++;
 786       }
 787     }
 788     alignedSeq.setSequence(newseq.toString());
 789     return alignedSeq;
 790   }
 791
 792   /**
 793    * Realigns the given protein to match the alignment of the dna, using codon
 794    * mappings to translate aligned codon positions to protein residues.
 795    *
 796    * @param protein
 797    *          the alignment whose sequences are realigned by this method
 798    * @param dna
 799    *          the dna alignment whose alignment we are 'copying'
 800    * @return the number of sequences that were realigned
 801    */
 802   public static int alignProteinAsDna(AlignmentI protein, AlignmentI dna)
 803   {
 804     Set<AlignedCodonFrame> mappings = protein.getCodonFrames();
 805
 806     /*
 807      * Map will hold, for each aligned codon position e.g. [3, 5, 6], a map of
 808      * {dnaSequence, {proteinSequence, codonProduct}} at that position. The
 809      * comparator keeps the codon positions ordered.
 810      */
 811     Map<AlignedCodon, Map<SequenceI, String>> alignedCodons = new TreeMap<AlignedCodon, Map<SequenceI, String>>(
 812             new CodonComparator());
 813     for (SequenceI dnaSeq : dna.getSequences())
 814     {
 815       for (AlignedCodonFrame mapping : mappings)
 816       {
 817         Mapping seqMap = mapping.getMappingForSequence(dnaSeq);
 818         SequenceI prot = mapping.findAlignedSequence(
 819                 dnaSeq.getDatasetSequence(), protein);
 820         if (prot != null)
 821         {
 822           addCodonPositions(dnaSeq, prot, protein.getGapCharacter(),
 823                   seqMap, alignedCodons);
 824         }
 825       }
 826     }
 827     return alignProteinAs(protein, alignedCodons);
 828   }
 829
 830   /**
 831    * Update the aligned protein sequences to match the codon alignments given in
 832    * the map.
 833    *
 834    * @param protein
 835    * @param alignedCodons
 836    *          an ordered map of codon positions (columns), with sequence/peptide
 837    *          values present in each column
 838    * @return
 839    */
 840   protected static int alignProteinAs(AlignmentI protein,
 841           Map<AlignedCodon, Map<SequenceI, String>> alignedCodons)
 842   {
 843     /*
 844      * Prefill aligned sequences with gaps before inserting aligned protein
 845      * residues.
 846      */
 847     int alignedWidth = alignedCodons.size();
 848     char[] gaps = new char[alignedWidth];
 849     Arrays.fill(gaps, protein.getGapCharacter());
 850     String allGaps = String.valueOf(gaps);
 851     for (SequenceI seq : protein.getSequences())
 852     {
 853       seq.setSequence(allGaps);
 854     }
 855
 856     int column = 0;
 857     for (AlignedCodon codon : alignedCodons.keySet())
 858     {
 859       final Map<SequenceI, String> columnResidues = alignedCodons.get(codon);
 860       for (Entry<SequenceI, String> entry : columnResidues
 861               .entrySet())
 862       {
 863         // place translated codon at its column position in sequence
 864         entry.getKey().getSequence()[column] = entry.getValue().charAt(0);
 865       }
 866       column++;
 867     }
 868     return 0;
 869   }
 870
 871   /**
 872    * Populate the map of aligned codons by traversing the given sequence
 873    * mapping, locating the aligned positions of mapped codons, and adding those
 874    * positions and their translation products to the map.
 875    *
 876    * @param dna
 877    *          the aligned sequence we are mapping from
 878    * @param protein
 879    *          the sequence to be aligned to the codons
 880    * @param gapChar
 881    *          the gap character in the dna sequence
 882    * @param seqMap
 883    *          a mapping to a sequence translation
 884    * @param alignedCodons
 885    *          the map we are building up
 886    */
 887   static void addCodonPositions(SequenceI dna, SequenceI protein,
 888           char gapChar,
 889           Mapping seqMap,
 890           Map<AlignedCodon, Map<SequenceI, String>> alignedCodons)
 891   {
 892     Iterator<AlignedCodon> codons = seqMap.getCodonIterator(dna, gapChar);
 893     while (codons.hasNext())
 894     {
 895       AlignedCodon codon = codons.next();
 896       Map<SequenceI, String> seqProduct = alignedCodons.get(codon);
 897       if (seqProduct == null)
 898       {
 899         seqProduct = new HashMap<SequenceI, String>();
 900         alignedCodons.put(codon, seqProduct);
 901       }
 902       seqProduct.put(protein, codon.product);
 903     }
 904   }
 905 }