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