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, aaResidue++)
 413     {
 414       String codon = String.valueOf(cdnaSeqChars, i, 3);
 415       final String translated = ResidueProperties.codonTranslate(
 416               codon);
 417       /*
 418        * ? allow X in protein to match untranslatable in dna ?
 419        */
 420       final char aaRes = aaSeqChars[aaResidue];
 421       if (translated == null && aaRes == 'X')
 422       {
 423         continue;
 424       }
 425       if (translated == null
 426               || !(aaRes == translated.charAt(0)))
 427       {
 428         return false;
 429       }
 430     }
 431     // fail if we didn't match all of the aa sequence
 432     return (aaResidue == aaSeqChars.length);
 433   }
 434
 435   /**
 436    * Align sequence 'seq' to match the alignment of a mapped sequence. Note this
 437    * currently assumes that we are aligning cDNA to match protein.
 438    *
 439    * @param seq
 440    *          the sequence to be realigned
 441    * @param al
 442    *          the alignment whose sequence alignment is to be 'copied'
 443    * @param gap
 444    *          character string represent a gap in the realigned sequence
 445    * @param preserveUnmappedGaps
 446    * @param preserveMappedGaps
 447    * @return true if the sequence was realigned, false if it could not be
 448    */
 449   public static boolean alignSequenceAs(SequenceI seq, AlignmentI al,
 450           String gap, boolean preserveMappedGaps,
 451           boolean preserveUnmappedGaps)
 452   {
 453     /*
 454      * Get any mappings from the source alignment to the target (dataset) sequence.
 455      */
 456     // TODO there may be one AlignedCodonFrame per dataset sequence, or one with
 457     // all mappings. Would it help to constrain this?
 458     List<AlignedCodonFrame> mappings = al.getCodonFrame(seq);
 459     if (mappings == null || mappings.isEmpty())
 460     {
 461       return false;
 462     }
 463
 464     /*
 465      * Locate the aligned source sequence whose dataset sequence is mapped. We
 466      * just take the first match here (as we can't align cDNA like more than one
 467      * protein sequence).
 468      */
 469     SequenceI alignFrom = null;
 470     AlignedCodonFrame mapping = null;
 471     for (AlignedCodonFrame mp : mappings)
 472     {
 473       alignFrom = mp.findAlignedSequence(seq.getDatasetSequence(), al);
 474       if (alignFrom != null)
 475       {
 476         mapping = mp;
 477         break;
 478       }
 479     }
 480
 481     if (alignFrom == null)
 482     {
 483       return false;
 484     }
 485     alignSequenceAs(seq, alignFrom, mapping, gap, al.getGapCharacter(),
 486             preserveMappedGaps, preserveUnmappedGaps);
 487     return true;
 488   }
 489
 490   /**
 491    * Align sequence 'alignTo' the same way as 'alignFrom', using the mapping to
 492    * match residues and codons. Flags control whether existing gaps in unmapped
 493    * (intron) and mapped (exon) regions are preserved or not. Gaps linking intro
 494    * and exon are only retained if both flags are set.
 495    *
 496    * @param alignTo
 497    * @param alignFrom
 498    * @param mapping
 499    * @param myGap
 500    * @param sourceGap
 501    * @param preserveUnmappedGaps
 502    * @param preserveMappedGaps
 503    */
 504   public static void alignSequenceAs(SequenceI alignTo,
 505           SequenceI alignFrom,
 506           AlignedCodonFrame mapping, String myGap, char sourceGap,
 507           boolean preserveMappedGaps, boolean preserveUnmappedGaps)
 508   {
 509     // TODO generalise to work for Protein-Protein, dna-dna, dna-protein
 510     final char[] thisSeq = alignTo.getSequence();
 511     final char[] thatAligned = alignFrom.getSequence();
 512     StringBuilder thisAligned = new StringBuilder(2 * thisSeq.length);
 513
 514     // aligned and dataset sequence positions, all base zero
 515     int thisSeqPos = 0;
 516     int sourceDsPos = 0;
 517
 518     int basesWritten = 0;
 519     char myGapChar = myGap.charAt(0);
 520     int ratio = myGap.length();
 521
 522     /*
 523      * Traverse the aligned protein sequence.
 524      */
 525     int sourceGapMappedLength = 0;
 526     boolean inExon = false;
 527     for (char sourceChar : thatAligned)
 528     {
 529       if (sourceChar == sourceGap)
 530       {
 531         sourceGapMappedLength += ratio;
 532         continue;
 533       }
 534
 535       /*
 536        * Found a residue. Locate its mapped codon (start) position.
 537        */
 538       sourceDsPos++;
 539       // Note mapping positions are base 1, our sequence positions base 0
 540       int[] mappedPos = mapping.getMappedRegion(alignTo, alignFrom,
 541               sourceDsPos);
 542       if (mappedPos == null)
 543       {
 544         /*
 545          * Abort realignment if unmapped protein. Or could ignore it??
 546          */
 547         System.err.println("Can't align: no codon mapping to residue "
 548                 + sourceDsPos + "(" + sourceChar + ")");
 549         return;
 550       }
 551
 552       int mappedCodonStart = mappedPos[0]; // position (1...) of codon start
 553       int mappedCodonEnd = mappedPos[mappedPos.length - 1]; // codon end pos
 554       StringBuilder trailingCopiedGap = new StringBuilder();
 555
 556       /*
 557        * Copy dna sequence up to and including this codon. Optionally, include
 558        * gaps before the codon starts (in introns) and/or after the codon starts
 559        * (in exons).
 560        *
 561        * Note this only works for 'linear' splicing, not reverse or interleaved.
 562        * But then 'align dna as protein' doesn't make much sense otherwise.
 563        */
 564       int intronLength = 0;
 565       while (basesWritten < mappedCodonEnd && thisSeqPos < thisSeq.length)
 566       {
 567         final char c = thisSeq[thisSeqPos++];
 568         if (c != myGapChar)
 569         {
 570           basesWritten++;
 571
 572           if (basesWritten < mappedCodonStart)
 573           {
 574             /*
 575              * Found an unmapped (intron) base. First add in any preceding gaps
 576              * (if wanted).
 577              */
 578             if (preserveUnmappedGaps && trailingCopiedGap.length() > 0)
 579             {
 580               thisAligned.append(trailingCopiedGap.toString());
 581               intronLength += trailingCopiedGap.length();
 582               trailingCopiedGap = new StringBuilder();
 583             }
 584             intronLength++;
 585             inExon = false;
 586           }
 587           else
 588           {
 589             final boolean startOfCodon = basesWritten == mappedCodonStart;
 590             int gapsToAdd = calculateGapsToInsert(preserveMappedGaps,
 591                     preserveUnmappedGaps, sourceGapMappedLength, inExon,
 592                     trailingCopiedGap.length(), intronLength, startOfCodon);
 593             for (int i = 0; i < gapsToAdd; i++)
 594             {
 595               thisAligned.append(myGapChar);
 596             }
 597             sourceGapMappedLength = 0;
 598             inExon = true;
 599           }
 600           thisAligned.append(c);
 601           trailingCopiedGap = new StringBuilder();
 602         }
 603         else
 604         {
 605           if (inExon && preserveMappedGaps)
 606           {
 607             trailingCopiedGap.append(myGapChar);
 608           }
 609           else if (!inExon && preserveUnmappedGaps)
 610           {
 611             trailingCopiedGap.append(myGapChar);
 612           }
 613         }
 614       }
 615     }
 616
 617     /*
 618      * At end of protein sequence. Copy any remaining dna sequence, optionally
 619      * including (intron) gaps. We do not copy trailing gaps in protein.
 620      */
 621     while (thisSeqPos < thisSeq.length)
 622     {
 623       final char c = thisSeq[thisSeqPos++];
 624       if (c != myGapChar || preserveUnmappedGaps)
 625       {
 626         thisAligned.append(c);
 627       }
 628     }
 629
 630     /*
 631      * All done aligning, set the aligned sequence.
 632      */
 633     alignTo.setSequence(new String(thisAligned));
 634   }
 635
 636   /**
 637    * Helper method to work out how many gaps to insert when realigning.
 638    *
 639    * @param preserveMappedGaps
 640    * @param preserveUnmappedGaps
 641    * @param sourceGapMappedLength
 642    * @param inExon
 643    * @param trailingCopiedGap
 644    * @param intronLength
 645    * @param startOfCodon
 646    * @return
 647    */
 648   protected static int calculateGapsToInsert(boolean preserveMappedGaps,
 649           boolean preserveUnmappedGaps, int sourceGapMappedLength,
 650           boolean inExon, int trailingGapLength,
 651           int intronLength, final boolean startOfCodon)
 652   {
 653     int gapsToAdd = 0;
 654     if (startOfCodon)
 655     {
 656       /*
 657        * Reached start of codon. Ignore trailing gaps in intron unless we are
 658        * preserving gaps in both exon and intron. Ignore them anyway if the
 659        * protein alignment introduces a gap at least as large as the intronic
 660        * region.
 661        */
 662       if (inExon && !preserveMappedGaps)
 663       {
 664         trailingGapLength = 0;
 665       }
 666       if (!inExon && !(preserveMappedGaps && preserveUnmappedGaps))
 667       {
 668         trailingGapLength = 0;
 669       }
 670       if (inExon)
 671       {
 672         gapsToAdd = Math.max(sourceGapMappedLength, trailingGapLength);
 673       }
 674       else
 675       {
 676         if (intronLength + trailingGapLength <= sourceGapMappedLength)
 677         {
 678           gapsToAdd = sourceGapMappedLength - intronLength;
 679         }
 680         else
 681         {
 682           gapsToAdd = Math.min(intronLength + trailingGapLength
 683                   - sourceGapMappedLength, trailingGapLength);
 684         }
 685       }
 686     }
 687     else
 688     {
 689       /*
 690        * second or third base of codon; check for any gaps in dna
 691        */
 692       if (!preserveMappedGaps)
 693       {
 694         trailingGapLength = 0;
 695       }
 696       gapsToAdd = Math.max(sourceGapMappedLength, trailingGapLength);
 697     }
 698     return gapsToAdd;
 699   }
 700
 701   /**
 702    * Returns a list of sequences mapped from the given sequences and aligned
 703    * (gapped) in the same way. For example, the cDNA for aligned protein, where
 704    * a single gap in protein generates three gaps in cDNA.
 705    *
 706    * @param sequences
 707    * @param gapCharacter
 708    * @param mappings
 709    * @return
 710    */
 711   public static List<SequenceI> getAlignedTranslation(
 712           List<SequenceI> sequences, char gapCharacter,
 713           Set<AlignedCodonFrame> mappings)
 714   {
 715     List<SequenceI> alignedSeqs = new ArrayList<SequenceI>();
 716
 717     for (SequenceI seq : sequences)
 718     {
 719       List<SequenceI> mapped = getAlignedTranslation(seq, gapCharacter,
 720               mappings);
 721       alignedSeqs.addAll(mapped);
 722     }
 723     return alignedSeqs;
 724   }
 725
 726   /**
 727    * Returns sequences aligned 'like' the source sequence, as mapped by the
 728    * given mappings. Normally we expect zero or one 'mapped' sequences, but this
 729    * will support 1-to-many as well.
 730    *
 731    * @param seq
 732    * @param gapCharacter
 733    * @param mappings
 734    * @return
 735    */
 736   protected static List<SequenceI> getAlignedTranslation(SequenceI seq,
 737           char gapCharacter, Set<AlignedCodonFrame> mappings)
 738   {
 739     List<SequenceI> result = new ArrayList<SequenceI>();
 740     for (AlignedCodonFrame mapping : mappings)
 741     {
 742       if (mapping.involvesSequence(seq))
 743       {
 744         SequenceI mapped = getAlignedTranslation(seq, gapCharacter, mapping);
 745         if (mapped != null)
 746         {
 747           result.add(mapped);
 748         }
 749       }
 750     }
 751     return result;
 752   }
 753
 754   /**
 755    * Returns the translation of 'seq' (as held in the mapping) with
 756    * corresponding alignment (gaps).
 757    *
 758    * @param seq
 759    * @param gapCharacter
 760    * @param mapping
 761    * @return
 762    */
 763   protected static SequenceI getAlignedTranslation(SequenceI seq,
 764           char gapCharacter, AlignedCodonFrame mapping)
 765   {
 766     String gap = String.valueOf(gapCharacter);
 767     boolean toDna = false;
 768     int fromRatio = 1;
 769     SequenceI mapTo = mapping.getDnaForAaSeq(seq);
 770     if (mapTo != null)
 771     {
 772       // mapping is from protein to nucleotide
 773       toDna = true;
 774       // should ideally get gap count ratio from mapping
 775       gap = String.valueOf(new char[]
 776       { gapCharacter, gapCharacter, gapCharacter });
 777     }
 778     else
 779     {
 780       // mapping is from nucleotide to protein
 781       mapTo = mapping.getAaForDnaSeq(seq);
 782       fromRatio = 3;
 783     }
 784     StringBuilder newseq = new StringBuilder(seq.getLength()
 785             * (toDna ? 3 : 1));
 786
 787     int residueNo = 0; // in seq, base 1
 788     int[] phrase = new int[fromRatio];
 789     int phraseOffset = 0;
 790     int gapWidth = 0;
 791     boolean first = true;
 792     final Sequence alignedSeq = new Sequence("", "");
 793
 794     for (char c : seq.getSequence())
 795     {
 796       if (c == gapCharacter)
 797       {
 798         gapWidth++;
 799         if (gapWidth >= fromRatio)
 800         {
 801           newseq.append(gap);
 802           gapWidth = 0;
 803         }
 804       }
 805       else
 806       {
 807         phrase[phraseOffset++] = residueNo + 1;
 808         if (phraseOffset == fromRatio)
 809         {
 810           /*
 811            * Have read a whole codon (or protein residue), now translate: map
 812            * source phrase to positions in target sequence add characters at
 813            * these positions to newseq Note mapping positions are base 1, our
 814            * sequence positions base 0.
 815            */
 816           SearchResults sr = new SearchResults();
 817           for (int pos : phrase)
 818           {
 819             mapping.markMappedRegion(seq, pos, sr);
 820           }
 821           newseq.append(sr.toString());
 822           if (first)
 823           {
 824             first = false;
 825             // Hack: Copy sequence dataset, name and description from
 826             // SearchResults.match[0].sequence
 827             // TODO? carry over sequence names from original 'complement'
 828             // alignment
 829             SequenceI mappedTo = sr.getResultSequence(0);
 830             alignedSeq.setName(mappedTo.getName());
 831             alignedSeq.setDescription(mappedTo.getDescription());
 832             alignedSeq.setDatasetSequence(mappedTo);
 833           }
 834           phraseOffset = 0;
 835         }
 836         residueNo++;
 837       }
 838     }
 839     alignedSeq.setSequence(newseq.toString());
 840     return alignedSeq;
 841   }
 842
 843   /**
 844    * Realigns the given protein to match the alignment of the dna, using codon
 845    * mappings to translate aligned codon positions to protein residues.
 846    *
 847    * @param protein
 848    *          the alignment whose sequences are realigned by this method
 849    * @param dna
 850    *          the dna alignment whose alignment we are 'copying'
 851    * @return the number of sequences that were realigned
 852    */
 853   public static int alignProteinAsDna(AlignmentI protein, AlignmentI dna)
 854   {
 855     Set<AlignedCodonFrame> mappings = protein.getCodonFrames();
 856
 857     /*
 858      * Map will hold, for each aligned codon position e.g. [3, 5, 6], a map of
 859      * {dnaSequence, {proteinSequence, codonProduct}} at that position. The
 860      * comparator keeps the codon positions ordered.
 861      */
 862     Map<AlignedCodon, Map<SequenceI, String>> alignedCodons = new TreeMap<AlignedCodon, Map<SequenceI, String>>(
 863             new CodonComparator());
 864     for (SequenceI dnaSeq : dna.getSequences())
 865     {
 866       for (AlignedCodonFrame mapping : mappings)
 867       {
 868         Mapping seqMap = mapping.getMappingForSequence(dnaSeq);
 869         SequenceI prot = mapping.findAlignedSequence(
 870                 dnaSeq.getDatasetSequence(), protein);
 871         if (prot != null)
 872         {
 873           addCodonPositions(dnaSeq, prot, protein.getGapCharacter(),
 874                   seqMap, alignedCodons);
 875         }
 876       }
 877     }
 878     return alignProteinAs(protein, alignedCodons);
 879   }
 880
 881   /**
 882    * Update the aligned protein sequences to match the codon alignments given in
 883    * the map.
 884    *
 885    * @param protein
 886    * @param alignedCodons
 887    *          an ordered map of codon positions (columns), with sequence/peptide
 888    *          values present in each column
 889    * @return
 890    */
 891   protected static int alignProteinAs(AlignmentI protein,
 892           Map<AlignedCodon, Map<SequenceI, String>> alignedCodons)
 893   {
 894     /*
 895      * Prefill aligned sequences with gaps before inserting aligned protein
 896      * residues.
 897      */
 898     int alignedWidth = alignedCodons.size();
 899     char[] gaps = new char[alignedWidth];
 900     Arrays.fill(gaps, protein.getGapCharacter());
 901     String allGaps = String.valueOf(gaps);
 902     for (SequenceI seq : protein.getSequences())
 903     {
 904       seq.setSequence(allGaps);
 905     }
 906
 907     int column = 0;
 908     for (AlignedCodon codon : alignedCodons.keySet())
 909     {
 910       final Map<SequenceI, String> columnResidues = alignedCodons.get(codon);
 911       for (Entry<SequenceI, String> entry : columnResidues
 912               .entrySet())
 913       {
 914         // place translated codon at its column position in sequence
 915         entry.getKey().getSequence()[column] = entry.getValue().charAt(0);
 916       }
 917       column++;
 918     }
 919     return 0;
 920   }
 921
 922   /**
 923    * Populate the map of aligned codons by traversing the given sequence
 924    * mapping, locating the aligned positions of mapped codons, and adding those
 925    * positions and their translation products to the map.
 926    *
 927    * @param dna
 928    *          the aligned sequence we are mapping from
 929    * @param protein
 930    *          the sequence to be aligned to the codons
 931    * @param gapChar
 932    *          the gap character in the dna sequence
 933    * @param seqMap
 934    *          a mapping to a sequence translation
 935    * @param alignedCodons
 936    *          the map we are building up
 937    */
 938   static void addCodonPositions(SequenceI dna, SequenceI protein,
 939           char gapChar,
 940           Mapping seqMap,
 941           Map<AlignedCodon, Map<SequenceI, String>> alignedCodons)
 942   {
 943     Iterator<AlignedCodon> codons = seqMap.getCodonIterator(dna, gapChar);
 944     while (codons.hasNext())
 945     {
 946       AlignedCodon codon = codons.next();
 947       Map<SequenceI, String> seqProduct = alignedCodons.get(codon);
 948       if (seqProduct == null)
 949       {
 950         seqProduct = new HashMap<SequenceI, String>();
 951         alignedCodons.put(codon, seqProduct);
 952       }
 953       seqProduct.put(protein, codon.product);
 954     }
 955   }
 956
 957   /**
 958    * Returns true if a cDNA/Protein mapping either exists, or could be made,
 959    * between at least one pair of sequences in the two alignments. Currently,
 960    * the logic is:
 961    * <ul>
 962    * <li>One alignment must be nucleotide, and the other protein</li>
 963    * <li>At least one pair of sequences must be already mapped, or mappable</li>
 964    * <li>Mappable means the nucleotide translation matches the protein sequence</li>
 965    * <li>The translation may ignore start and stop codons if present in the
 966    * nucleotide</li>
 967    * </ul>
 968    *
 969    * @param al1
 970    * @param al2
 971    * @return
 972    */
 973   public static boolean isMappable(AlignmentI al1, AlignmentI al2)
 974   {
 975     /*
 976      * Require one nucleotide and one protein
 977      */
 978     if (al1.isNucleotide() == al2.isNucleotide())
 979     {
 980       return false;
 981     }
 982     AlignmentI dna = al1.isNucleotide() ? al1 : al2;
 983     AlignmentI protein = dna == al1 ? al2 : al1;
 984     Set<AlignedCodonFrame> mappings = protein.getCodonFrames();
 985     for (SequenceI dnaSeq : dna.getSequences())
 986     {
 987       for (SequenceI proteinSeq : protein.getSequences())
 988       {
 989         if (isMappable(dnaSeq, proteinSeq, mappings))
 990         {
 991           return true;
 992         }
 993       }
 994     }
 995     return false;
 996   }
 997
 998   /**
 999    * Returns true if the dna sequence is mapped, or could be mapped, to the
1000    * protein sequence.
1001    *
1002    * @param dnaSeq
1003    * @param proteinSeq
1004    * @param mappings
1005    * @return
1006    */
1007   public static boolean isMappable(SequenceI dnaSeq, SequenceI proteinSeq,
1008           Set<AlignedCodonFrame> mappings)
1009   {
1010     SequenceI dnaDs = dnaSeq.getDatasetSequence() == null ? dnaSeq : dnaSeq.getDatasetSequence();
1011     SequenceI proteinDs = proteinSeq.getDatasetSequence() == null ? proteinSeq
1012             : proteinSeq.getDatasetSequence();
1013
1014     /*
1015      * Already mapped?
1016      */
1017     for (AlignedCodonFrame mapping : mappings) {
1018       if ( proteinDs == mapping.getAaForDnaSeq(dnaDs)) {
1019         return true;
1020       }
1021     }
1022
1023     /*
1024      * Just try to make a mapping (it is not yet stored), test whether
1025      * successful.
1026      */
1027     return mapProteinToCdna(proteinDs, dnaDs) != null;
1028   }
1029 }