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