src/jalview/analysis/Dna.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.bin.Cache;
  24 import jalview.datamodel.AlignedCodonFrame;
  25 import jalview.datamodel.Alignment;
  26 import jalview.datamodel.AlignmentAnnotation;
  27 import jalview.datamodel.AlignmentI;
  28 import jalview.datamodel.Annotation;
  29 import jalview.datamodel.DBRefEntry;
  30 import jalview.datamodel.DBRefSource;
  31 import jalview.datamodel.FeatureProperties;
  32 import jalview.datamodel.GraphLine;
  33 import jalview.datamodel.Mapping;
  34 import jalview.datamodel.Sequence;
  35 import jalview.datamodel.SequenceFeature;
  36 import jalview.datamodel.SequenceI;
  37 import jalview.schemes.ResidueProperties;
  38 import jalview.util.Comparison;
  39 import jalview.util.DBRefUtils;
  40 import jalview.util.MapList;
  41 import jalview.util.ShiftList;
  42
  43 import java.util.ArrayList;
  44 import java.util.Arrays;
  45 import java.util.Hashtable;
  46 import java.util.List;
  47
  48 public class Dna
  49 {
  50   /**
  51    * Test whether codon positions cdp1 should align before, with, or after cdp2.
  52    * Returns zero if all positions match (or either argument is null). Returns
  53    * -1 if any position in the first codon precedes the corresponding position
  54    * in the second codon. Else returns +1 (some position in the second codon
  55    * precedes the corresponding position in the first).
  56    *
  57    * Note this is not necessarily symmetric, for example:
  58    * <ul>
  59    * <li>compareCodonPos([2,5,6], [3,4,5]) returns -1</li>
  60    * <li>compareCodonPos([3,4,5], [2,5,6]) also returns -1</li>
  61    * </ul>
  62    *
  63    * @param cdp1
  64    * @param cdp2
  65    * @return
  66    */
  67   public static int compareCodonPos(int[] cdp1, int[] cdp2)
  68   {
  69     if (cdp1 == null
  70             || cdp2 == null
  71             || (cdp1[0] == cdp2[0] && cdp1[1] == cdp2[1] && cdp1[2] == cdp2[2]))
  72     {
  73       return 0;
  74     }
  75     if (cdp1[0] < cdp2[0] || cdp1[1] < cdp2[1] || cdp1[2] < cdp2[2])
  76     {
  77       // one base in cdp1 precedes the corresponding base in the other codon
  78       return -1;
  79     }
  80     // one base in cdp1 appears after the corresponding base in the other codon.
  81     return 1;
  82   }
  83
  84   /**
  85    * DNA->mapped protein sequence alignment translation given set of sequences
  86    * 1. id distinct coding regions within selected region for each sequence 2.
  87    * generate peptides based on inframe (or given) translation or (optionally
  88    * and where specified) out of frame translations (annotated appropriately) 3.
  89    * align peptides based on codon alignment
  90    */
  91   /**
  92    * id potential products from dna 1. search for distinct products within
  93    * selected region for each selected sequence 2. group by associated DB type.
  94    * 3. return as form for input into above function
  95    */
  96   /**
  97    *
  98    */
  99   /**
 100    * create a new alignment of protein sequences by an inframe translation of
 101    * the provided NA sequences
 102    *
 103    * @param selection
 104    * @param seqstring
 105    * @param viscontigs
 106    * @param gapCharacter
 107    * @param annotations
 108    * @param aWidth
 109    * @param dataset
 110    *          destination dataset for translated sequences and mappings
 111    * @return
 112    */
 113   public static AlignmentI cdnaTranslate(SequenceI[] selection,
 114           String[] seqstring, int viscontigs[], char gapCharacter,
 115           AlignmentAnnotation[] annotations, int aWidth, Alignment dataset)
 116   {
 117     return cdnaTranslate(Arrays.asList(selection), seqstring, null,
 118             viscontigs, gapCharacter, annotations, aWidth, dataset);
 119   }
 120
 121   /**
 122    *
 123    * @param cdnaseqs
 124    * @param seqstring
 125    * @param product
 126    *          - array of DbRefEntry objects from which exon map in seqstring is
 127    *          derived
 128    * @param viscontigs
 129    * @param gapCharacter
 130    * @param annotations
 131    * @param aWidth
 132    * @param dataset
 133    * @return
 134    */
 135   public static AlignmentI cdnaTranslate(List<SequenceI> cdnaseqs,
 136           String[] seqstring, DBRefEntry[] product, int viscontigs[],
 137           char gapCharacter, AlignmentAnnotation[] annotations, int aWidth,
 138           Alignment dataset)
 139   {
 140     AlignedCodonFrame acf = new AlignedCodonFrame(aWidth);
 141
 142     /*
 143      * This array will be built up so that position i holds the codon positions
 144      * e.g. [7, 9, 10] that match column i (base 0) in the aligned translation.
 145      * Note this implies a contract that if two codons do not align exactly,
 146      * their translated products must occupy different column positions.
 147      */
 148     int[][] alignedCodons = new int[aWidth][];
 149
 150     int s;
 151     int sSize = cdnaseqs.size();
 152     List<SequenceI> pepseqs = new ArrayList<SequenceI>();
 153     for (s = 0; s < sSize; s++)
 154     {
 155       SequenceI newseq = translateCodingRegion(cdnaseqs.get(s),
 156               seqstring[s], viscontigs, acf, alignedCodons, pepseqs,
 157               gapCharacter,
 158               false);
 159
 160       if (newseq != null)
 161       {
 162         pepseqs.add(newseq);
 163         SequenceI ds = newseq;
 164         if (dataset != null)
 165         {
 166           while (ds.getDatasetSequence() != null)
 167           {
 168             ds = ds.getDatasetSequence();
 169           }
 170           dataset.addSequence(ds);
 171         }
 172       }
 173     }
 174
 175     SequenceI[] newseqs = pepseqs.toArray(new SequenceI[pepseqs.size()]);
 176     AlignmentI al = new Alignment(newseqs);
 177     // ensure we look aligned.
 178     al.padGaps();
 179     // link the protein translation to the DNA dataset
 180     al.setDataset(dataset);
 181     translateAlignedAnnotations(annotations, al, acf, alignedCodons);
 182     al.addCodonFrame(acf);
 183     return al;
 184   }
 185
 186   /**
 187    * fake the collection of DbRefs with associated exon mappings to identify if
 188    * a translation would generate distinct product in the currently selected
 189    * region.
 190    *
 191    * @param selection
 192    * @param viscontigs
 193    * @return
 194    */
 195   public static boolean canTranslate(SequenceI[] selection,
 196           int viscontigs[])
 197   {
 198     for (int gd = 0; gd < selection.length; gd++)
 199     {
 200       SequenceI dna = selection[gd];
 201       DBRefEntry[] dnarefs = DBRefUtils
 202               .selectRefs(dna.getDBRef(),
 203                       jalview.datamodel.DBRefSource.DNACODINGDBS);
 204       if (dnarefs != null)
 205       {
 206         // intersect with pep
 207         List<DBRefEntry> mappedrefs = new ArrayList<DBRefEntry>();
 208         DBRefEntry[] refs = dna.getDBRef();
 209         for (int d = 0; d < refs.length; d++)
 210         {
 211           if (refs[d].getMap() != null && refs[d].getMap().getMap() != null
 212                   && refs[d].getMap().getMap().getFromRatio() == 3
 213                   && refs[d].getMap().getMap().getToRatio() == 1)
 214           {
 215             mappedrefs.add(refs[d]); // add translated protein maps
 216           }
 217         }
 218         dnarefs = mappedrefs.toArray(new DBRefEntry[mappedrefs.size()]);
 219         for (int d = 0; d < dnarefs.length; d++)
 220         {
 221           Mapping mp = dnarefs[d].getMap();
 222           if (mp != null)
 223           {
 224             for (int vc = 0; vc < viscontigs.length; vc += 2)
 225             {
 226               int[] mpr = mp.locateMappedRange(viscontigs[vc],
 227                       viscontigs[vc + 1]);
 228               if (mpr != null)
 229               {
 230                 return true;
 231               }
 232             }
 233           }
 234         }
 235       }
 236     }
 237     return false;
 238   }
 239
 240   /**
 241    * Translate na alignment annotations onto translated amino acid alignment al
 242    * using codon mapping codons
 243    *
 244    * @param annotations
 245    * @param al
 246    * @param acf
 247    */
 248   protected static void translateAlignedAnnotations(
 249           AlignmentAnnotation[] annotations, AlignmentI al,
 250           AlignedCodonFrame acf, int[][] codons)
 251   {
 252     // Can only do this for columns with consecutive codons, or where
 253     // annotation is sequence associated.
 254
 255     if (annotations != null)
 256     {
 257       for (AlignmentAnnotation annotation : annotations)
 258       {
 259         /*
 260          * Skip hidden or autogenerated annotation. Also (for now), RNA
 261          * secondary structure annotation. If we want to show this against
 262          * protein we need a smarter way to 'translate' without generating
 263          * invalid (unbalanced) structure annotation.
 264          */
 265         if (annotation.autoCalculated || !annotation.visible
 266                 || annotation.isRNA())
 267         {
 268           continue;
 269         }
 270
 271         int aSize = acf.getaaWidth(); // aa alignment width.
 272         Annotation[] anots = (annotation.annotations == null) ? null
 273                 : new Annotation[aSize];
 274         if (anots != null)
 275         {
 276           for (int a = 0; a < aSize; a++)
 277           {
 278             // process through codon map.
 279             if (a < codons.length && codons[a] != null
 280                     && codons[a][0] == (codons[a][2] - 2))
 281             {
 282               anots[a] = getCodonAnnotation(codons[a],
 283                       annotation.annotations);
 284             }
 285           }
 286         }
 287
 288         AlignmentAnnotation aa = new AlignmentAnnotation(annotation.label,
 289                 annotation.description, anots);
 290         aa.graph = annotation.graph;
 291         aa.graphGroup = annotation.graphGroup;
 292         aa.graphHeight = annotation.graphHeight;
 293         if (annotation.getThreshold() != null)
 294         {
 295           aa.setThreshold(new GraphLine(annotation
 296                   .getThreshold()));
 297         }
 298         if (annotation.hasScore)
 299         {
 300           aa.setScore(annotation.getScore());
 301         }
 302
 303         final SequenceI seqRef = annotation.sequenceRef;
 304         if (seqRef != null)
 305         {
 306           SequenceI aaSeq = acf.getAaForDnaSeq(seqRef);
 307           if (aaSeq != null)
 308           {
 309             // aa.compactAnnotationArray(); // throw away alignment annotation
 310             // positioning
 311             aa.setSequenceRef(aaSeq);
 312             // rebuild mapping
 313             aa.createSequenceMapping(aaSeq, aaSeq.getStart(), true);
 314             aa.adjustForAlignment();
 315             aaSeq.addAlignmentAnnotation(aa);
 316           }
 317         }
 318         al.addAnnotation(aa);
 319       }
 320     }
 321   }
 322
 323   private static Annotation getCodonAnnotation(int[] is,
 324           Annotation[] annotations)
 325   {
 326     // Have a look at all the codon positions for annotation and put the first
 327     // one found into the translated annotation pos.
 328     int contrib = 0;
 329     Annotation annot = null;
 330     for (int p = 0; p < 3; p++)
 331     {
 332       if (annotations[is[p]] != null)
 333       {
 334         if (annot == null)
 335         {
 336           annot = new Annotation(annotations[is[p]]);
 337           contrib = 1;
 338         }
 339         else
 340         {
 341           // merge with last
 342           Annotation cpy = new Annotation(annotations[is[p]]);
 343           if (annot.colour == null)
 344           {
 345             annot.colour = cpy.colour;
 346           }
 347           if (annot.description == null || annot.description.length() == 0)
 348           {
 349             annot.description = cpy.description;
 350           }
 351           if (annot.displayCharacter == null)
 352           {
 353             annot.displayCharacter = cpy.displayCharacter;
 354           }
 355           if (annot.secondaryStructure == 0)
 356           {
 357             annot.secondaryStructure = cpy.secondaryStructure;
 358           }
 359           annot.value += cpy.value;
 360           contrib++;
 361         }
 362       }
 363     }
 364     if (contrib > 1)
 365     {
 366       annot.value /= contrib;
 367     }
 368     return annot;
 369   }
 370
 371   /**
 372    * Translate a na sequence
 373    *
 374    * @param selection
 375    *          sequence displayed under viscontigs visible columns
 376    * @param seqstring
 377    *          ORF read in some global alignment reference frame
 378    * @param viscontigs
 379    *          mapping from global reference frame to visible seqstring ORF read
 380    * @param acf
 381    *          Definition of global ORF alignment reference frame
 382    * @param alignedCodons
 383    * @param proteinSeqs
 384    * @param gapCharacter
 385    * @param starForStop
 386    *          when true stop codons will translate as '*', otherwise as 'X'
 387    * @return sequence ready to be added to alignment.
 388    */
 389   protected static SequenceI translateCodingRegion(SequenceI selection,
 390           String seqstring, int[] viscontigs, AlignedCodonFrame acf,
 391           int[][] alignedCodons, List<SequenceI> proteinSeqs,
 392           char gapCharacter,
 393           final boolean starForStop)
 394   {
 395     List<int[]> skip = new ArrayList<int[]>();
 396     int skipint[] = null;
 397     ShiftList vismapping = new ShiftList(); // map from viscontigs to seqstring
 398     // intervals
 399     int vc, scontigs[] = new int[viscontigs.length];
 400     int npos = 0;
 401     for (vc = 0; vc < viscontigs.length; vc += 2)
 402     {
 403       if (vc == 0)
 404       {
 405         vismapping.addShift(npos, viscontigs[vc]);
 406       }
 407       else
 408       {
 409         // hidden region
 410         vismapping.addShift(npos, viscontigs[vc] - viscontigs[vc - 1] + 1);
 411       }
 412       scontigs[vc] = viscontigs[vc];
 413       scontigs[vc + 1] = viscontigs[vc + 1];
 414     }
 415
 416     // allocate a roughly sized buffer for the protein sequence
 417     StringBuilder protein = new StringBuilder(seqstring.length() / 2);
 418     String seq = seqstring.replace('U', 'T').replace('u', 'T');
 419     char codon[] = new char[3];
 420     int cdp[] = new int[3], rf = 0, lastnpos = 0, nend;
 421     int aspos = 0;
 422     int resSize = 0;
 423     for (npos = 0, nend = seq.length(); npos < nend; npos++)
 424     {
 425       if (!Comparison.isGap(seq.charAt(npos)))
 426       {
 427         cdp[rf] = npos; // store position
 428         codon[rf++] = seq.charAt(npos); // store base
 429       }
 430       if (rf == 3)
 431       {
 432         /*
 433          * Filled up a reading frame...
 434          */
 435         String aa = ResidueProperties.codonTranslate(new String(codon));
 436         rf = 0;
 437         final String gapString = String.valueOf(gapCharacter);
 438         if (aa == null)
 439         {
 440           aa = gapString;
 441           if (skipint == null)
 442           {
 443             skipint = new int[]
 444             { cdp[0], cdp[2] };
 445           }
 446           skipint[1] = cdp[2];
 447         }
 448         else
 449         {
 450           if (skipint != null)
 451           {
 452             // edit scontigs
 453             skipint[0] = vismapping.shift(skipint[0]);
 454             skipint[1] = vismapping.shift(skipint[1]);
 455             for (vc = 0; vc < scontigs.length;)
 456             {
 457               if (scontigs[vc + 1] < skipint[0])
 458               {
 459                 // before skipint starts
 460                 vc += 2;
 461                 continue;
 462               }
 463               if (scontigs[vc] > skipint[1])
 464               {
 465                 // finished editing so
 466                 break;
 467               }
 468               // Edit the contig list to include the skipped region which did
 469               // not translate
 470               int[] t;
 471               // from : s1 e1 s2 e2 s3 e3
 472               // to s: s1 e1 s2 k0 k1 e2 s3 e3
 473               // list increases by one unless one boundary (s2==k0 or e2==k1)
 474               // matches, and decreases by one if skipint intersects whole
 475               // visible contig
 476               if (scontigs[vc] <= skipint[0])
 477               {
 478                 if (skipint[0] == scontigs[vc])
 479                 {
 480                   // skipint at start of contig
 481                   // shift the start of this contig
 482                   if (scontigs[vc + 1] > skipint[1])
 483                   {
 484                     scontigs[vc] = skipint[1];
 485                     vc += 2;
 486                   }
 487                   else
 488                   {
 489                     if (scontigs[vc + 1] == skipint[1])
 490                     {
 491                       // remove the contig
 492                       t = new int[scontigs.length - 2];
 493                       if (vc > 0)
 494                       {
 495                         System.arraycopy(scontigs, 0, t, 0, vc - 1);
 496                       }
 497                       if (vc + 2 < t.length)
 498                       {
 499                         System.arraycopy(scontigs, vc + 2, t, vc, t.length
 500                                 - vc + 2);
 501                       }
 502                       scontigs = t;
 503                     }
 504                     else
 505                     {
 506                       // truncate contig to before the skipint region
 507                       scontigs[vc + 1] = skipint[0] - 1;
 508                       vc += 2;
 509                     }
 510                   }
 511                 }
 512                 else
 513                 {
 514                   // scontig starts before start of skipint
 515                   if (scontigs[vc + 1] < skipint[1])
 516                   {
 517                     // skipint truncates end of scontig
 518                     scontigs[vc + 1] = skipint[0] - 1;
 519                     vc += 2;
 520                   }
 521                   else
 522                   {
 523                     // divide region to new contigs
 524                     t = new int[scontigs.length + 2];
 525                     System.arraycopy(scontigs, 0, t, 0, vc + 1);
 526                     t[vc + 1] = skipint[0];
 527                     t[vc + 2] = skipint[1];
 528                     System.arraycopy(scontigs, vc + 1, t, vc + 3,
 529                             scontigs.length - (vc + 1));
 530                     scontigs = t;
 531                     vc += 4;
 532                   }
 533                 }
 534               }
 535             }
 536             skip.add(skipint);
 537             skipint = null;
 538           }
 539           if (aa.equals("STOP"))
 540           {
 541             aa = starForStop ? "*" : "X";
 542           }
 543           resSize++;
 544         }
 545         // insert gaps prior to this codon - if necessary
 546         boolean findpos = true;
 547         while (findpos)
 548         {
 549           // expand the codons array if necessary
 550           alignedCodons = checkCodonFrameWidth(alignedCodons, aspos);
 551
 552           /*
 553            * Compare this codon's base positions with those currently aligned to
 554            * this column in the translation.
 555            */
 556           final int compareCodonPos = Dna.compareCodonPos(cdp,
 557                   alignedCodons[aspos]);
 558           // debug
 559           System.out.println(seq + "/" + aa + " codons: "
 560                   + Arrays.deepToString(alignedCodons));
 561           System.out
 562                   .println(("Compare " + Arrays.toString(cdp) + " at pos "
 563                           + aspos + " with "
 564                           + Arrays.toString(alignedCodons[aspos]) + " got " + compareCodonPos));
 565           // end debug
 566           switch (compareCodonPos)
 567           {
 568           case -1:
 569
 570             /*
 571              * This codon should precede the mapped positions - need to insert a
 572              * gap in all prior sequences.
 573              */
 574             insertAAGap(aspos, gapCharacter, alignedCodons, proteinSeqs);
 575             findpos = false;
 576             break;
 577
 578           case +1:
 579
 580             /*
 581              * This codon belongs after the aligned codons at aspos. Prefix it
 582              * with a gap and try the next position.
 583              */
 584             aa = gapString + aa;
 585             aspos++;
 586             break;
 587
 588           case 0:
 589
 590             /*
 591              * Exact match - codon 'belongs' at this translated position.
 592              */
 593             findpos = false;
 594           }
 595         }
 596         protein.append(aa);
 597         lastnpos = npos;
 598         if (alignedCodons[aspos] == null)
 599         {
 600           // mark this column as aligning to this aligned reading frame
 601           alignedCodons[aspos] = new int[]
 602           { cdp[0], cdp[1], cdp[2] };
 603         }
 604         else if (!Arrays.equals(alignedCodons[aspos], cdp))
 605         {
 606           throw new IllegalStateException("Tried to coalign "
 607                   + Arrays.asList(alignedCodons[aspos], cdp));
 608         }
 609         if (aspos >= acf.aaWidth)
 610         {
 611           // update maximum alignment width
 612           // (we can do this without calling checkCodonFrameWidth because it was
 613           // already done above)
 614           System.out.println("aspos " + aspos + " >= " + acf.aaWidth);
 615           acf.setAaWidth(aspos);
 616         }
 617         // ready for next translated reading frame alignment position (if any)
 618         aspos++;
 619       }
 620     }
 621     if (resSize > 0)
 622     {
 623       SequenceI newseq = new Sequence(selection.getName(),
 624               protein.toString());
 625       if (rf != 0)
 626       {
 627         final String errMsg = "trimming contigs for incomplete terminal codon.";
 628         if (Cache.log != null)
 629         {
 630           Cache.log.debug(errMsg);
 631         }
 632         else
 633         {
 634           System.err.println(errMsg);
 635         }
 636         // map and trim contigs to ORF region
 637         vc = scontigs.length - 1;
 638         lastnpos = vismapping.shift(lastnpos); // place npos in context of
 639         // whole dna alignment (rather
 640         // than visible contigs)
 641         // incomplete ORF could be broken over one or two visible contig
 642         // intervals.
 643         while (vc >= 0 && scontigs[vc] > lastnpos)
 644         {
 645           if (vc > 0 && scontigs[vc - 1] > lastnpos)
 646           {
 647             vc -= 2;
 648           }
 649           else
 650           {
 651             // correct last interval in list.
 652             scontigs[vc] = lastnpos;
 653           }
 654         }
 655
 656         if (vc > 0 && (vc + 1) < scontigs.length)
 657         {
 658           // truncate map list to just vc elements
 659           int t[] = new int[vc + 1];
 660           System.arraycopy(scontigs, 0, t, 0, vc + 1);
 661           scontigs = t;
 662         }
 663         if (vc <= 0)
 664         {
 665           scontigs = null;
 666         }
 667       }
 668       if (scontigs != null)
 669       {
 670         npos = 0;
 671         // map scontigs to actual sequence positions on selection
 672         for (vc = 0; vc < scontigs.length; vc += 2)
 673         {
 674           scontigs[vc] = selection.findPosition(scontigs[vc]); // not from 1!
 675           scontigs[vc + 1] = selection.findPosition(scontigs[vc + 1]); // exclusive
 676           if (scontigs[vc + 1] == selection.getEnd())
 677           {
 678             break;
 679           }
 680         }
 681         // trim trailing empty intervals.
 682         if ((vc + 2) < scontigs.length)
 683         {
 684           int t[] = new int[vc + 2];
 685           System.arraycopy(scontigs, 0, t, 0, vc + 2);
 686           scontigs = t;
 687         }
 688         /*
 689          * delete intervals in scontigs which are not translated. 1. map skip
 690          * into sequence position intervals 2. truncate existing ranges and add
 691          * new ranges to exclude untranslated regions. if (skip.size()>0) {
 692          * Vector narange = new Vector(); for (vc=0; vc<scontigs.length; vc++) {
 693          * narange.addElement(new int[] {scontigs[vc]}); } int sint=0,iv[]; vc =
 694          * 0; while (sint<skip.size()) { skipint = (int[]) skip.elementAt(sint);
 695          * do { iv = (int[]) narange.elementAt(vc); if (iv[0]>=skipint[0] &&
 696          * iv[0]<=skipint[1]) { if (iv[0]==skipint[0]) { // delete beginning of
 697          * range } else { // truncate range and create new one if necessary iv =
 698          * (int[]) narange.elementAt(vc+1); if (iv[0]<=skipint[1]) { // truncate
 699          * range iv[0] = skipint[1]; } else { } } } else if (iv[0]<skipint[0]) {
 700          * iv = (int[]) narange.elementAt(vc+1); } } while (iv[0]) } }
 701          */
 702         MapList map = new MapList(scontigs, new int[]
 703         { 1, resSize }, 3, 1);
 704
 705         transferCodedFeatures(selection, newseq, map, null, null);
 706
 707         /*
 708          * Construct a dataset sequence for our new peptide.
 709          */
 710         SequenceI rseq = newseq.deriveSequence();
 711
 712         /*
 713          * Store a mapping (between the dataset sequences for the two
 714          * sequences).
 715          */
 716         // SIDE-EFFECT: acf stores the aligned sequence reseq; to remove!
 717         acf.addMap(selection, rseq, map);
 718         return rseq;
 719       }
 720     }
 721     // register the mapping somehow
 722     //
 723     return null;
 724   }
 725
 726   /**
 727    * Insert a gap into the aligned proteins and the codon mapping array.
 728    *
 729    * @param pos
 730    * @param gapCharacter
 731    * @param alignedCodons
 732    * @param proteinSeqs
 733    */
 734   protected static void insertAAGap(int pos, char gapCharacter,
 735           int[][] alignedCodons, List<SequenceI> proteinSeqs)
 736   {
 737     System.out.println("insertAAGap " + pos + "/" + proteinSeqs.size());
 738     // aaWidth++;
 739     for (SequenceI seq : proteinSeqs)
 740     {
 741       seq.insertCharAt(pos, gapCharacter);
 742     }
 743
 744     // if (pos < aaWidth)
 745     // {
 746     // aaWidth++;
 747     System.arraycopy(alignedCodons, pos, alignedCodons, pos + 1,
 748             alignedCodons.length - pos - 1);
 749     alignedCodons[pos] = null; // clear so new codon position can be marked.
 750     // }
 751   }
 752
 753   /**
 754    * Check the codons array is big enough to accommodate the given position, if
 755    * not resize it.
 756    *
 757    * @param alignedCodons
 758    * @param aspos
 759    * @return the resized array (or the original if no resize needed)
 760    */
 761   protected static int[][] checkCodonFrameWidth(int[][] alignedCodons,
 762           int aspos)
 763   {
 764     // TODO why not codons.length < aspos ?
 765     // should codons expand if length is 2 or 3 and aslen==2 ?
 766     if (alignedCodons.length <= aspos + 1)
 767     {
 768       // probably never have to do this ?
 769       int[][] c = new int[alignedCodons.length + 10][];
 770       for (int i = 0; i < alignedCodons.length; i++)
 771       {
 772         c[i] = alignedCodons[i];
 773       }
 774       return c;
 775     }
 776     return alignedCodons;
 777   }
 778
 779   /**
 780    * Given a peptide newly translated from a dna sequence, copy over and set any
 781    * features on the peptide from the DNA. If featureTypes is null, all features
 782    * on the dna sequence are searched (rather than just the displayed ones), and
 783    * similarly for featureGroups.
 784    *
 785    * @param dna
 786    * @param pep
 787    * @param map
 788    * @param featureTypes
 789    *          hash who's keys are the displayed feature type strings
 790    * @param featureGroups
 791    *          hash where keys are feature groups and values are Boolean objects
 792    *          indicating if they are displayed.
 793    */
 794   private static void transferCodedFeatures(SequenceI dna, SequenceI pep,
 795           MapList map, Hashtable featureTypes, Hashtable featureGroups)
 796   {
 797     SequenceFeature[] sf = (dna.getDatasetSequence() != null ? dna
 798             .getDatasetSequence() : dna).getSequenceFeatures();
 799     Boolean fgstate;
 800     DBRefEntry[] dnarefs = DBRefUtils.selectRefs(dna.getDBRef(),
 801             DBRefSource.DNACODINGDBS);
 802     if (dnarefs != null)
 803     {
 804       // intersect with pep
 805       for (int d = 0; d < dnarefs.length; d++)
 806       {
 807         Mapping mp = dnarefs[d].getMap();
 808         if (mp != null)
 809         {
 810         }
 811       }
 812     }
 813     if (sf != null)
 814     {
 815       for (int f = 0; f < sf.length; f++)
 816       {
 817         fgstate = (featureGroups == null) ? null : ((Boolean) featureGroups
 818                 .get(sf[f].featureGroup));
 819         if ((featureTypes == null || featureTypes.containsKey(sf[f]
 820                 .getType())) && (fgstate == null || fgstate.booleanValue()))
 821         {
 822           if (FeatureProperties.isCodingFeature(null, sf[f].getType()))
 823           {
 824             // if (map.intersectsFrom(sf[f].begin, sf[f].end))
 825             {
 826
 827             }
 828           }
 829         }
 830       }
 831     }
 832   }
 833 }