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