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