JAL-3949 - refactor logging from jalview.bin.Cache to jalview.bin.Console
[jalview.git] / src / jalview / analysis / AlignmentUtils.java
index 77cc4d6..28f23b2 100644 (file)
  */
 package jalview.analysis;
 
-import static jalview.io.gff.GffConstants.CLINICAL_SIGNIFICANCE;
+import java.util.Locale;
 
+import java.util.ArrayList;
+import java.util.Arrays;
+import java.util.Collection;
+import java.util.Collections;
+import java.util.HashMap;
+import java.util.HashSet;
+import java.util.Iterator;
+import java.util.LinkedHashMap;
+import java.util.List;
+import java.util.Map;
+import java.util.Map.Entry;
+import java.util.NoSuchElementException;
+import java.util.Set;
+import java.util.SortedMap;
+import java.util.TreeMap;
+
+import jalview.bin.Console;
+import jalview.commands.RemoveGapColCommand;
 import jalview.datamodel.AlignedCodon;
 import jalview.datamodel.AlignedCodonFrame;
 import jalview.datamodel.AlignedCodonFrame.SequenceToSequenceMapping;
@@ -37,7 +55,6 @@ import jalview.datamodel.SequenceFeature;
 import jalview.datamodel.SequenceGroup;
 import jalview.datamodel.SequenceI;
 import jalview.datamodel.features.SequenceFeatures;
-import jalview.io.gff.Gff3Helper;
 import jalview.io.gff.SequenceOntologyI;
 import jalview.schemes.ResidueProperties;
 import jalview.util.Comparison;
@@ -45,25 +62,6 @@ import jalview.util.DBRefUtils;
 import jalview.util.IntRangeComparator;
 import jalview.util.MapList;
 import jalview.util.MappingUtils;
-import jalview.util.StringUtils;
-
-import java.io.UnsupportedEncodingException;
-import java.net.URLEncoder;
-import java.util.ArrayList;
-import java.util.Arrays;
-import java.util.Collection;
-import java.util.Collections;
-import java.util.HashMap;
-import java.util.HashSet;
-import java.util.Iterator;
-import java.util.LinkedHashMap;
-import java.util.List;
-import java.util.Map;
-import java.util.Map.Entry;
-import java.util.NoSuchElementException;
-import java.util.Set;
-import java.util.SortedMap;
-import java.util.TreeMap;
 
 /**
  * grab bag of useful alignment manipulation operations Expect these to be
@@ -186,9 +184,9 @@ public class AlignmentUtils
       // TODO use Character.toLowerCase to avoid creating String objects?
       char[] upstream = new String(ds
               .getSequence(s.getStart() - 1 - ustream_ds, s.getStart() - 1))
-                      .toLowerCase().toCharArray();
+                      .toLowerCase(Locale.ROOT).toCharArray();
       char[] downstream = new String(
-              ds.getSequence(s_end - 1, s_end + dstream_ds)).toLowerCase()
+              ds.getSequence(s_end - 1, s_end + dstream_ds)).toLowerCase(Locale.ROOT)
                       .toCharArray();
       char[] coreseq = s.getSequence();
       char[] nseq = new char[offset + upstream.length + downstream.length
@@ -467,7 +465,7 @@ public class AlignmentUtils
     if (cdnaLength != mappedLength && cdnaLength > 2)
     {
       String lastCodon = String.valueOf(cdnaSeqChars,
-              cdnaLength - CODON_LENGTH, CODON_LENGTH).toUpperCase();
+              cdnaLength - CODON_LENGTH, CODON_LENGTH).toUpperCase(Locale.ROOT);
       for (String stop : ResidueProperties.STOP_CODONS)
       {
         if (lastCodon.equals(stop))
@@ -484,7 +482,7 @@ public class AlignmentUtils
      */
     int startOffset = 0;
     if (cdnaLength != mappedLength && cdnaLength > 2
-            && String.valueOf(cdnaSeqChars, 0, CODON_LENGTH).toUpperCase()
+            && String.valueOf(cdnaSeqChars, 0, CODON_LENGTH).toUpperCase(Locale.ROOT)
                     .equals(ResidueProperties.START))
     {
       startOffset += CODON_LENGTH;
@@ -1602,11 +1600,12 @@ public class AlignmentUtils
       return false;
     }
     String name = seq2.getName();
-    final DBRefEntry[] xrefs = seq1.getDBRefs();
+    final List<DBRefEntry> xrefs = seq1.getDBRefs();
     if (xrefs != null)
     {
-      for (DBRefEntry xref : xrefs)
+      for (int ix = 0, nx = xrefs.size(); ix < nx; ix++)
       {
+        DBRefEntry xref = xrefs.get(ix);
         String xrefName = xref.getSource() + "|" + xref.getAccessionId();
         // case-insensitive test, consistent with DBRefEntry.equalRef()
         if (xrefName.equalsIgnoreCase(name))
@@ -1736,15 +1735,8 @@ public class AlignmentUtils
 
           cdsSeqs.add(cdsSeq);
 
-          if (!dataset.getSequences().contains(cdsSeqDss))
-          {
-            // check if this sequence is a newly created one
-            // so needs adding to the dataset
-            dataset.addSequence(cdsSeqDss);
-          }
-
           /*
-           * add a mapping from CDS to the (unchanged) mapped to range
+           * build the mapping from CDS to protein
            */
           List<int[]> cdsRange = Collections
                   .singletonList(new int[]
@@ -1753,16 +1745,26 @@ public class AlignmentUtils
           MapList cdsToProteinMap = new MapList(cdsRange,
                   mapList.getToRanges(), mapList.getFromRatio(),
                   mapList.getToRatio());
-          AlignedCodonFrame cdsToProteinMapping = new AlignedCodonFrame();
-          cdsToProteinMapping.addMap(cdsSeqDss, proteinProduct,
-                  cdsToProteinMap);
 
-          /*
-           * guard against duplicating the mapping if repeating this action
-           */
-          if (!mappings.contains(cdsToProteinMapping))
+          if (!dataset.getSequences().contains(cdsSeqDss))
           {
-            mappings.add(cdsToProteinMapping);
+            /*
+             * if this sequence is a newly created one, add it to the dataset
+             * and made a CDS to protein mapping (if sequence already exists,
+             * CDS-to-protein mapping _is_ the transcript-to-protein mapping)
+             */
+            dataset.addSequence(cdsSeqDss);
+            AlignedCodonFrame cdsToProteinMapping = new AlignedCodonFrame();
+            cdsToProteinMapping.addMap(cdsSeqDss, proteinProduct,
+                  cdsToProteinMap);
+
+            /*
+             * guard against duplicating the mapping if repeating this action
+             */
+            if (!mappings.contains(cdsToProteinMapping))
+            {
+              mappings.add(cdsToProteinMapping);
+            }
           }
 
           propagateDBRefsToCDS(cdsSeqDss, dnaSeq.getDatasetSequence(),
@@ -1804,8 +1806,10 @@ public class AlignmentUtils
           // need to
           // synthesize an xref.
 
-          for (DBRefEntry primRef : dnaDss.getPrimaryDBRefs())
+          List<DBRefEntry> primrefs = dnaDss.getPrimaryDBRefs();
+          for (int ip = 0, np = primrefs.size(); ip < np; ip++)
           {
+                 DBRefEntry primRef = primrefs.get(ip);
             /*
              * create a cross-reference from CDS to the source sequence's
              * primary reference and vice versa
@@ -1819,7 +1823,6 @@ public class AlignmentUtils
 
             dnaSeq.addDBRef(new DBRefEntry(source, version, cdsSeq
                     .getName(), new Mapping(cdsSeqDss, dnaToCdsMap)));
-
             // problem here is that the cross-reference is synthesized -
             // cdsSeq.getName() may be like 'CDS|dnaaccession' or
             // 'CDS|emblcdsacc'
@@ -1832,7 +1835,6 @@ public class AlignmentUtils
                     .getInverse()));
             proteinProduct.addDBRef(proteinToCdsRef);
           }
-
           /*
            * transfer any features on dna that overlap the CDS
            */
@@ -1996,45 +1998,31 @@ public class AlignmentUtils
 
     SequenceI newSeq = null;
 
-    final MapList maplist = mapping.getMap();
-    if (maplist.isContiguous() && maplist.isFromForwardStrand())
-    {
-      /*
-       * just a subsequence, keep same dataset sequence
-       */
-      int start = maplist.getFromLowest();
-      int end = maplist.getFromHighest();
-      newSeq = seq.getSubSequence(start - 1, end);
-      newSeq.setName(seqId);
-    }
-    else
-    {
-      /*
-       * construct by splicing mapped from ranges
-       */
-      char[] seqChars = seq.getSequence();
-      List<int[]> fromRanges = maplist.getFromRanges();
-      int cdsWidth = MappingUtils.getLength(fromRanges);
-      char[] newSeqChars = new char[cdsWidth];
+    /*
+     * construct CDS sequence by splicing mapped from ranges
+     */
+    char[] seqChars = seq.getSequence();
+    List<int[]> fromRanges = mapping.getMap().getFromRanges();
+    int cdsWidth = MappingUtils.getLength(fromRanges);
+    char[] newSeqChars = new char[cdsWidth];
 
-      int newPos = 0;
-      for (int[] range : fromRanges)
+    int newPos = 0;
+    for (int[] range : fromRanges)
+    {
+      if (range[0] <= range[1])
       {
-        if (range[0] <= range[1])
-        {
-          // forward strand mapping - just copy the range
-          int length = range[1] - range[0] + 1;
-          System.arraycopy(seqChars, range[0] - 1, newSeqChars, newPos,
-                  length);
-          newPos += length;
-        }
-        else
+        // forward strand mapping - just copy the range
+        int length = range[1] - range[0] + 1;
+        System.arraycopy(seqChars, range[0] - 1, newSeqChars, newPos,
+                length);
+        newPos += length;
+      }
+      else
+      {
+        // reverse strand mapping - copy and complement one by one
+        for (int i = range[0]; i >= range[1]; i--)
         {
-          // reverse strand mapping - copy and complement one by one
-          for (int i = range[0]; i >= range[1]; i--)
-          {
-            newSeqChars[newPos++] = Dna.getComplement(seqChars[i - 1]);
-          }
+          newSeqChars[newPos++] = Dna.getComplement(seqChars[i - 1]);
         }
       }
 
@@ -2068,9 +2056,8 @@ public class AlignmentUtils
           }
           else
           {
-            System.err.println(
-                    "JAL-2154 regression: warning - found (and ignnored a duplicate CDS sequence):"
-                            + mtch.toString());
+            Console.error(
+                    "JAL-2154 regression: warning - found (and ignored) a duplicate CDS sequence:" + mtch.toString());
           }
         }
       }
@@ -2098,13 +2085,15 @@ public class AlignmentUtils
     List<DBRefEntry> direct = new ArrayList<>();
     HashSet<String> directSources = new HashSet<>();
 
-    if (contig.getDBRefs() != null)
+    List<DBRefEntry> refs = contig.getDBRefs();
+    if (refs != null)
     {
-      for (DBRefEntry dbr : contig.getDBRefs())
+      for (int ib = 0, nb = refs.size(); ib < nb; ib++)
       {
-        if (dbr.hasMap() && dbr.getMap().getMap().isTripletMap())
+        DBRefEntry dbr = refs.get(ib);
+        MapList map;
+        if (dbr.hasMap() && (map = dbr.getMap().getMap()).isTripletMap())
         {
-          MapList map = dbr.getMap().getMap();
           // check if map is the CDS mapping
           if (mapping.getMap().equals(map))
           {
@@ -2114,21 +2103,22 @@ public class AlignmentUtils
         }
       }
     }
-    DBRefEntry[] onSource = DBRefUtils.selectRefs(
+    List<DBRefEntry> onSource = DBRefUtils.selectRefs(
             proteinProduct.getDBRefs(),
             directSources.toArray(new String[0]));
     List<DBRefEntry> propagated = new ArrayList<>();
 
     // and generate appropriate mappings
-    for (DBRefEntry cdsref : direct)
+    for (int ic = 0, nc = direct.size(); ic < nc; ic++)
     {
+      DBRefEntry cdsref = direct.get(ic);
+      Mapping m = cdsref.getMap();
       // clone maplist and mapping
       MapList cdsposmap = new MapList(
               Arrays.asList(new int[][]
               { new int[] { cdsSeq.getStart(), cdsSeq.getEnd() } }),
-              cdsref.getMap().getMap().getToRanges(), 3, 1);
-      Mapping cdsmap = new Mapping(cdsref.getMap().getTo(),
-              cdsref.getMap().getMap());
+              m.getMap().getToRanges(), 3, 1);
+      Mapping cdsmap = new Mapping(m.getTo(), m.getMap());
 
       // create dbref
       DBRefEntry newref = new DBRefEntry(cdsref.getSource(),
@@ -2351,10 +2341,14 @@ public class AlignmentUtils
       int phase = 0;
       try
       {
-        phase = Integer.parseInt(sf.getPhase());
+       String s = sf.getPhase();
+       if (s != null) 
+       {
+               phase = Integer.parseInt(s);
+       }
       } catch (NumberFormatException e)
       {
-        // ignore
+        // leave as zero
       }
       /*
        * phase > 0 on first codon means 5' incomplete - skip to the start
@@ -2388,394 +2382,6 @@ public class AlignmentUtils
   }
 
   /**
-   * Maps exon features from dna to protein, and computes variants in peptide
-   * product generated by variants in dna, and adds them as sequence_variant
-   * features on the protein sequence. Returns the number of variant features
-   * added.
-   * 
-   * @param dnaSeq
-   * @param peptide
-   * @param dnaToProtein
-   */
-  public static int computeProteinFeatures(SequenceI dnaSeq,
-          SequenceI peptide, MapList dnaToProtein)
-  {
-    while (dnaSeq.getDatasetSequence() != null)
-    {
-      dnaSeq = dnaSeq.getDatasetSequence();
-    }
-    while (peptide.getDatasetSequence() != null)
-    {
-      peptide = peptide.getDatasetSequence();
-    }
-
-    transferFeatures(dnaSeq, peptide, dnaToProtein, SequenceOntologyI.EXON);
-
-    /*
-     * compute protein variants from dna variants and codon mappings;
-     * NB - alternatively we could retrieve this using the REST service e.g.
-     * http://rest.ensembl.org/overlap/translation
-     * /ENSP00000288602?feature=transcript_variation;content-type=text/xml
-     * which would be a bit slower but possibly more reliable
-     */
-
-    /*
-     * build a map with codon variations for each potentially varying peptide
-     */
-    LinkedHashMap<Integer, List<DnaVariant>[]> variants = buildDnaVariantsMap(
-            dnaSeq, dnaToProtein);
-
-    /*
-     * scan codon variations, compute peptide variants and add to peptide sequence
-     */
-    int count = 0;
-    for (Entry<Integer, List<DnaVariant>[]> variant : variants.entrySet())
-    {
-      int peptidePos = variant.getKey();
-      List<DnaVariant>[] codonVariants = variant.getValue();
-      count += computePeptideVariants(peptide, peptidePos, codonVariants);
-    }
-
-    return count;
-  }
-
-  /**
-   * Computes non-synonymous peptide variants from codon variants and adds them
-   * as sequence_variant features on the protein sequence (one feature per
-   * allele variant). Selected attributes (variant id, clinical significance)
-   * are copied over to the new features.
-   * 
-   * @param peptide
-   *          the protein sequence
-   * @param peptidePos
-   *          the position to compute peptide variants for
-   * @param codonVariants
-   *          a list of dna variants per codon position
-   * @return the number of features added
-   */
-  static int computePeptideVariants(SequenceI peptide, int peptidePos,
-          List<DnaVariant>[] codonVariants)
-  {
-    String residue = String
-            .valueOf(peptide.getCharAt(peptidePos - peptide.getStart()));
-    int count = 0;
-    String base1 = codonVariants[0].get(0).base;
-    String base2 = codonVariants[1].get(0).base;
-    String base3 = codonVariants[2].get(0).base;
-
-    /*
-     * variants in first codon base
-     */
-    for (DnaVariant var : codonVariants[0])
-    {
-      if (var.variant != null)
-      {
-        String alleles = (String) var.variant.getValue(Gff3Helper.ALLELES);
-        if (alleles != null)
-        {
-          for (String base : alleles.split(","))
-          {
-            if (!base1.equalsIgnoreCase(base))
-            {
-              String codon = base.toUpperCase() + base2.toLowerCase()
-                      + base3.toLowerCase();
-              String canonical = base1.toUpperCase() + base2.toLowerCase()
-                      + base3.toLowerCase();
-              if (addPeptideVariant(peptide, peptidePos, residue, var,
-                      codon, canonical))
-              {
-                count++;
-              }
-            }
-          }
-        }
-      }
-    }
-
-    /*
-     * variants in second codon base
-     */
-    for (DnaVariant var : codonVariants[1])
-    {
-      if (var.variant != null)
-      {
-        String alleles = (String) var.variant.getValue(Gff3Helper.ALLELES);
-        if (alleles != null)
-        {
-          for (String base : alleles.split(","))
-          {
-            if (!base2.equalsIgnoreCase(base))
-            {
-              String codon = base1.toLowerCase() + base.toUpperCase()
-                      + base3.toLowerCase();
-              String canonical = base1.toLowerCase() + base2.toUpperCase()
-                      + base3.toLowerCase();
-              if (addPeptideVariant(peptide, peptidePos, residue, var,
-                      codon, canonical))
-              {
-                count++;
-              }
-            }
-          }
-        }
-      }
-    }
-
-    /*
-     * variants in third codon base
-     */
-    for (DnaVariant var : codonVariants[2])
-    {
-      if (var.variant != null)
-      {
-        String alleles = (String) var.variant.getValue(Gff3Helper.ALLELES);
-        if (alleles != null)
-        {
-          for (String base : alleles.split(","))
-          {
-            if (!base3.equalsIgnoreCase(base))
-            {
-              String codon = base1.toLowerCase() + base2.toLowerCase()
-                      + base.toUpperCase();
-              String canonical = base1.toLowerCase() + base2.toLowerCase()
-                      + base3.toUpperCase();
-              if (addPeptideVariant(peptide, peptidePos, residue, var,
-                      codon, canonical))
-              {
-                count++;
-              }
-            }
-          }
-        }
-      }
-    }
-
-    return count;
-  }
-
-  /**
-   * Helper method that adds a peptide variant feature. ID and
-   * clinical_significance attributes of the dna variant (if present) are copied
-   * to the new feature.
-   * 
-   * @param peptide
-   * @param peptidePos
-   * @param residue
-   * @param var
-   * @param codon
-   *          the variant codon e.g. aCg
-   * @param canonical
-   *          the 'normal' codon e.g. aTg
-   * @return true if a feature was added, else false
-   */
-  static boolean addPeptideVariant(SequenceI peptide, int peptidePos,
-          String residue, DnaVariant var, String codon, String canonical)
-  {
-    /*
-     * get peptide translation of codon e.g. GAT -> D
-     * note that variants which are not single alleles,
-     * e.g. multibase variants or HGMD_MUTATION etc
-     * are currently ignored here
-     */
-    String trans = codon.contains("-") ? null
-            : (codon.length() > CODON_LENGTH ? null
-                    : ResidueProperties.codonTranslate(codon));
-    if (trans == null)
-    {
-      return false;
-    }
-    String desc = canonical + "/" + codon;
-    String featureType = "";
-    if (trans.equals(residue))
-    {
-      featureType = SequenceOntologyI.SYNONYMOUS_VARIANT;
-    }
-    else if (ResidueProperties.STOP.equals(trans))
-    {
-      featureType = SequenceOntologyI.STOP_GAINED;
-    }
-    else
-    {
-      String residue3Char = StringUtils
-              .toSentenceCase(ResidueProperties.aa2Triplet.get(residue));
-      String trans3Char = StringUtils
-              .toSentenceCase(ResidueProperties.aa2Triplet.get(trans));
-      desc = "p." + residue3Char + peptidePos + trans3Char;
-      featureType = SequenceOntologyI.NONSYNONYMOUS_VARIANT;
-    }
-    SequenceFeature sf = new SequenceFeature(featureType, desc, peptidePos,
-            peptidePos, var.getSource());
-
-    StringBuilder attributes = new StringBuilder(32);
-    String id = (String) var.variant.getValue(VARIANT_ID);
-    if (id != null)
-    {
-      if (id.startsWith(SEQUENCE_VARIANT))
-      {
-        id = id.substring(SEQUENCE_VARIANT.length());
-      }
-      sf.setValue(VARIANT_ID, id);
-      attributes.append(VARIANT_ID).append("=").append(id);
-      // TODO handle other species variants JAL-2064
-      StringBuilder link = new StringBuilder(32);
-      try
-      {
-        link.append(desc).append(" ").append(id).append(
-                "|http://www.ensembl.org/Homo_sapiens/Variation/Summary?v=")
-                .append(URLEncoder.encode(id, "UTF-8"));
-        sf.addLink(link.toString());
-      } catch (UnsupportedEncodingException e)
-      {
-        // as if
-      }
-    }
-    String clinSig = (String) var.variant.getValue(CLINICAL_SIGNIFICANCE);
-    if (clinSig != null)
-    {
-      sf.setValue(CLINICAL_SIGNIFICANCE, clinSig);
-      attributes.append(";").append(CLINICAL_SIGNIFICANCE).append("=")
-              .append(clinSig);
-    }
-    peptide.addSequenceFeature(sf);
-    if (attributes.length() > 0)
-    {
-      sf.setAttributes(attributes.toString());
-    }
-    return true;
-  }
-
-  /**
-   * Builds a map whose key is position in the protein sequence, and value is a
-   * list of the base and all variants for each corresponding codon position.
-   * <p>
-   * This depends on dna variants being held as a comma-separated list as
-   * property "alleles" on variant features.
-   * 
-   * @param dnaSeq
-   * @param dnaToProtein
-   * @return
-   */
-  @SuppressWarnings("unchecked")
-  static LinkedHashMap<Integer, List<DnaVariant>[]> buildDnaVariantsMap(
-          SequenceI dnaSeq, MapList dnaToProtein)
-  {
-    /*
-     * map from peptide position to all variants of the codon which codes for it
-     * LinkedHashMap ensures we keep the peptide features in sequence order
-     */
-    LinkedHashMap<Integer, List<DnaVariant>[]> variants = new LinkedHashMap<>();
-
-    List<SequenceFeature> dnaFeatures = dnaSeq.getFeatures()
-            .getFeaturesByOntology(SequenceOntologyI.SEQUENCE_VARIANT);
-    if (dnaFeatures.isEmpty())
-    {
-      return variants;
-    }
-
-    int dnaStart = dnaSeq.getStart();
-    int[] lastCodon = null;
-    int lastPeptidePostion = 0;
-
-    /*
-     * build a map of codon variations for peptides
-     */
-    for (SequenceFeature sf : dnaFeatures)
-    {
-      int dnaCol = sf.getBegin();
-      if (dnaCol != sf.getEnd())
-      {
-        // not handling multi-locus variant features
-        continue;
-      }
-
-      /*
-       * ignore variant if not a SNP
-       */
-      String alls = (String) sf.getValue(Gff3Helper.ALLELES);
-      if (alls == null)
-      {
-        continue; // non-SNP VCF variant perhaps - can't process this
-      }
-
-      String[] alleles = alls.toUpperCase().split(",");
-      boolean isSnp = true;
-      for (String allele : alleles)
-      {
-        if (allele.trim().length() > 1)
-        {
-          isSnp = false;
-        }
-      }
-      if (!isSnp)
-      {
-        continue;
-      }
-
-      int[] mapsTo = dnaToProtein.locateInTo(dnaCol, dnaCol);
-      if (mapsTo == null)
-      {
-        // feature doesn't lie within coding region
-        continue;
-      }
-      int peptidePosition = mapsTo[0];
-      List<DnaVariant>[] codonVariants = variants.get(peptidePosition);
-      if (codonVariants == null)
-      {
-        codonVariants = new ArrayList[CODON_LENGTH];
-        codonVariants[0] = new ArrayList<>();
-        codonVariants[1] = new ArrayList<>();
-        codonVariants[2] = new ArrayList<>();
-        variants.put(peptidePosition, codonVariants);
-      }
-
-      /*
-       * get this peptide's codon positions e.g. [3, 4, 5] or [4, 7, 10]
-       */
-      int[] codon = peptidePosition == lastPeptidePostion ? lastCodon
-              : MappingUtils.flattenRanges(dnaToProtein.locateInFrom(
-                      peptidePosition, peptidePosition));
-      lastPeptidePostion = peptidePosition;
-      lastCodon = codon;
-
-      /*
-       * save nucleotide (and any variant) for each codon position
-       */
-      for (int codonPos = 0; codonPos < CODON_LENGTH; codonPos++)
-      {
-        String nucleotide = String.valueOf(
-                dnaSeq.getCharAt(codon[codonPos] - dnaStart)).toUpperCase();
-        List<DnaVariant> codonVariant = codonVariants[codonPos];
-        if (codon[codonPos] == dnaCol)
-        {
-          if (!codonVariant.isEmpty()
-                  && codonVariant.get(0).variant == null)
-          {
-            /*
-             * already recorded base value, add this variant
-             */
-            codonVariant.get(0).variant = sf;
-          }
-          else
-          {
-            /*
-             * add variant with base value
-             */
-            codonVariant.add(new DnaVariant(nucleotide, sf));
-          }
-        }
-        else if (codonVariant.isEmpty())
-        {
-          /*
-           * record (possibly non-varying) base value
-           */
-          codonVariant.add(new DnaVariant(nucleotide));
-        }
-      }
-    }
-    return variants;
-  }
-
-  /**
    * Makes an alignment with a copy of the given sequences, adding in any
    * non-redundant sequences which are mapped to by the cross-referenced
    * sequences.
@@ -2795,19 +2401,25 @@ public class AlignmentUtils
     SequenceIdMatcher matcher = new SequenceIdMatcher(seqs);
     if (xrefs != null)
     {
-      for (SequenceI xref : xrefs)
+       // BH 2019.01.25 recoded to remove iterators
+       
+      for (int ix = 0, nx = xrefs.length; ix < nx; ix++)
       {
-        DBRefEntry[] dbrefs = xref.getDBRefs();
+       SequenceI xref = xrefs[ix];
+        List<DBRefEntry> dbrefs = xref.getDBRefs();
         if (dbrefs != null)
         {
-          for (DBRefEntry dbref : dbrefs)
+          for (int ir = 0, nir = dbrefs.size(); ir < nir; ir++)
           {
-            if (dbref.getMap() == null || dbref.getMap().getTo() == null
-                    || dbref.getMap().getTo().isProtein() != isProtein)
+            DBRefEntry dbref = dbrefs.get(ir);
+            Mapping map = dbref.getMap();
+            SequenceI mto;
+            if (map == null || (mto = map.getTo()) == null
+                    || mto.isProtein() != isProtein)
             {
               continue;
             }
-            SequenceI mappedTo = dbref.getMap().getTo();
+            SequenceI mappedTo = mto;
             SequenceI match = matcher.findIdMatch(mappedTo);
             if (match == null)
             {
@@ -2837,6 +2449,17 @@ public class AlignmentUtils
    */
   public static int alignAs(AlignmentI unaligned, AlignmentI aligned)
   {
+    /*
+     * easy case - aligning a copy of aligned sequences
+     */
+    if (alignAsSameSequences(unaligned, aligned))
+    {
+      return unaligned.getHeight();
+    }
+
+    /*
+     * fancy case - aligning via mappings between sequences
+     */
     List<SequenceI> unmapped = new ArrayList<>();
     Map<Integer, Map<SequenceI, Character>> columnMap = buildMappedColumnsMap(
             unaligned, aligned, unmapped);
@@ -2902,9 +2525,7 @@ public class AlignmentUtils
    *                    - 'guide' alignment containing sequences derived from same
    *                    dataset as unaligned
    * @return
-   * @deprecated probably obsolete and incomplete
    */
-  @Deprecated
   static boolean alignAsSameSequences(AlignmentI unaligned,
           AlignmentI aligned)
   {
@@ -2954,6 +2575,13 @@ public class AlignmentUtils
     {
       List<SequenceI> alignedSequences = alignedDatasets
               .get(seq.getDatasetSequence());
+      if (alignedSequences.isEmpty())
+      {
+        /*
+         * defensive check - shouldn't happen! (JAL-3536)
+         */
+        continue;
+      }
       SequenceI alignedSeq = alignedSequences.get(0);
 
       /*
@@ -2975,6 +2603,12 @@ public class AlignmentUtils
       }
     }
 
+    /*
+     * finally remove gapped columns (e.g. introns)
+     */
+    new RemoveGapColCommand("", unaligned.getSequencesArray(), 0,
+            unaligned.getWidth() - 1, unaligned);
+
     return true;
   }