JAL-3187 CDS has same dataset sequence as transcript
[jalview.git] / src / jalview / analysis / AlignmentUtils.java
index dfb185a..77cc4d6 100644 (file)
@@ -29,6 +29,7 @@ import jalview.datamodel.Alignment;
 import jalview.datamodel.AlignmentAnnotation;
 import jalview.datamodel.AlignmentI;
 import jalview.datamodel.DBRefEntry;
+import jalview.datamodel.GeneLociI;
 import jalview.datamodel.IncompleteCodonException;
 import jalview.datamodel.Mapping;
 import jalview.datamodel.Sequence;
@@ -36,6 +37,7 @@ 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;
@@ -72,12 +74,15 @@ import java.util.TreeMap;
  */
 public class AlignmentUtils
 {
-
   private static final int CODON_LENGTH = 3;
 
   private static final String SEQUENCE_VARIANT = "sequence_variant:";
 
-  private static final String ID = "ID";
+  /*
+   * the 'id' attribute is provided for variant features fetched from
+   * Ensembl using its REST service with JSON format 
+   */
+  public static final String VARIANT_ID = "id";
 
   /**
    * A data model to hold the 'normal' base value at a position, and an optional
@@ -105,6 +110,15 @@ public class AlignmentUtils
     {
       return variant == null ? null : variant.getFeatureGroup();
     }
+
+    /**
+     * toString for aid in the debugger only
+     */
+    @Override
+    public String toString()
+    {
+      return base + ":" + (variant == null ? "" : variant.getDescription());
+    }
   }
 
   /**
@@ -384,7 +398,7 @@ public class AlignmentUtils
    * Answers true if the mappings include one between the given (dataset)
    * sequences.
    */
-  public static boolean mappingExists(List<AlignedCodonFrame> mappings,
+  protected static boolean mappingExists(List<AlignedCodonFrame> mappings,
           SequenceI aaSeq, SequenceI cdnaSeq)
   {
     if (mappings != null)
@@ -454,7 +468,7 @@ public class AlignmentUtils
     {
       String lastCodon = String.valueOf(cdnaSeqChars,
               cdnaLength - CODON_LENGTH, CODON_LENGTH).toUpperCase();
-      for (String stop : ResidueProperties.STOP)
+      for (String stop : ResidueProperties.STOP_CODONS)
       {
         if (lastCodon.equals(stop))
         {
@@ -525,7 +539,8 @@ public class AlignmentUtils
        * allow * in protein to match untranslatable in dna
        */
       final char aaRes = aaSeqChars[aaPos];
-      if ((translated == null || "STOP".equals(translated)) && aaRes == '*')
+      if ((translated == null || ResidueProperties.STOP.equals(translated))
+              && aaRes == '*')
       {
         continue;
       }
@@ -557,7 +572,8 @@ public class AlignmentUtils
     if (dnaPos == cdnaSeqChars.length - CODON_LENGTH)
     {
       String codon = String.valueOf(cdnaSeqChars, dnaPos, CODON_LENGTH);
-      if ("STOP".equals(ResidueProperties.codonTranslate(codon)))
+      if (ResidueProperties.STOP
+              .equals(ResidueProperties.codonTranslate(codon)))
       {
         return true;
       }
@@ -1636,8 +1652,8 @@ public class AlignmentUtils
       productSeqs = new HashSet<>();
       for (SequenceI seq : products)
       {
-        productSeqs.add(seq.getDatasetSequence() == null ? seq
-                : seq.getDatasetSequence());
+        productSeqs.add(seq.getDatasetSequence() == null ? seq : seq
+                .getDatasetSequence());
       }
     }
 
@@ -1732,7 +1748,8 @@ public class AlignmentUtils
            */
           List<int[]> cdsRange = Collections
                   .singletonList(new int[]
-                  { 1, cdsSeq.getLength() });
+                  { cdsSeq.getStart(),
+                      cdsSeq.getLength() + cdsSeq.getStart() - 1 });
           MapList cdsToProteinMap = new MapList(cdsRange,
                   mapList.getToRanges(), mapList.getFromRatio(),
                   mapList.getToRatio());
@@ -1754,7 +1771,7 @@ public class AlignmentUtils
            * add another mapping from original 'from' range to CDS
            */
           AlignedCodonFrame dnaToCdsMapping = new AlignedCodonFrame();
-          MapList dnaToCdsMap = new MapList(mapList.getFromRanges(),
+          final MapList dnaToCdsMap = new MapList(mapList.getFromRanges(),
                   cdsRange, 1, 1);
           dnaToCdsMapping.addMap(dnaSeq.getDatasetSequence(), cdsSeqDss,
                   dnaToCdsMap);
@@ -1764,6 +1781,13 @@ public class AlignmentUtils
           }
 
           /*
+           * transfer dna chromosomal loci (if known) to the CDS
+           * sequence (via the mapping)
+           */
+          final MapList cdsToDnaMap = dnaToCdsMap.getInverse();
+          transferGeneLoci(dnaSeq, cdsToDnaMap, cdsSeq);
+
+          /*
            * add DBRef with mapping from protein to CDS
            * (this enables Get Cross-References from protein alignment)
            * This is tricky because we can't have two DBRefs with the
@@ -1782,26 +1806,30 @@ public class AlignmentUtils
 
           for (DBRefEntry primRef : dnaDss.getPrimaryDBRefs())
           {
-            // creates a complementary cross-reference to the source sequence's
-            // primary reference.
-
-            DBRefEntry cdsCrossRef = new DBRefEntry(primRef.getSource(),
-                    primRef.getSource() + ":" + primRef.getVersion(),
-                    primRef.getAccessionId());
-            cdsCrossRef
-                    .setMap(new Mapping(dnaDss, new MapList(dnaToCdsMap)));
+            /*
+             * create a cross-reference from CDS to the source sequence's
+             * primary reference and vice versa
+             */
+            String source = primRef.getSource();
+            String version = primRef.getVersion();
+            DBRefEntry cdsCrossRef = new DBRefEntry(source, source + ":"
+                    + version, primRef.getAccessionId());
+            cdsCrossRef.setMap(new Mapping(dnaDss, new MapList(cdsToDnaMap)));
             cdsSeqDss.addDBRef(cdsCrossRef);
 
+            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'
             // assuming cds version same as dna ?!?
 
-            DBRefEntry proteinToCdsRef = new DBRefEntry(primRef.getSource(),
-                    primRef.getVersion(), cdsSeq.getName());
+            DBRefEntry proteinToCdsRef = new DBRefEntry(source, version,
+                    cdsSeq.getName());
             //
-            proteinToCdsRef.setMap(
-                    new Mapping(cdsSeqDss, cdsToProteinMap.getInverse()));
+            proteinToCdsRef.setMap(new Mapping(cdsSeqDss, cdsToProteinMap
+                    .getInverse()));
             proteinProduct.addDBRef(proteinToCdsRef);
           }
 
@@ -1814,14 +1842,46 @@ public class AlignmentUtils
       }
     }
 
-    AlignmentI cds = new Alignment(
-            cdsSeqs.toArray(new SequenceI[cdsSeqs.size()]));
+    AlignmentI cds = new Alignment(cdsSeqs.toArray(new SequenceI[cdsSeqs
+            .size()]));
     cds.setDataset(dataset);
 
     return cds;
   }
 
   /**
+   * Tries to transfer gene loci (dbref to chromosome positions) from fromSeq to
+   * toSeq, mediated by the given mapping between the sequences
+   * 
+   * @param fromSeq
+   * @param targetToFrom
+   *          Map
+   * @param targetSeq
+   */
+  protected static void transferGeneLoci(SequenceI fromSeq,
+          MapList targetToFrom, SequenceI targetSeq)
+  {
+    if (targetSeq.getGeneLoci() != null)
+    {
+      // already have - don't override
+      return;
+    }
+    GeneLociI fromLoci = fromSeq.getGeneLoci();
+    if (fromLoci == null)
+    {
+      return;
+    }
+
+    MapList newMap = targetToFrom.traverse(fromLoci.getMapping());
+
+    if (newMap != null)
+    {
+      targetSeq.setGeneLoci(fromLoci.getSpeciesId(),
+              fromLoci.getAssemblyId(), fromLoci.getChromosomeId(), newMap);
+    }
+  }
+
+  /**
    * A helper method that finds a CDS sequence in the alignment dataset that is
    * mapped to the given protein sequence, and either is, or has a mapping from,
    * the given dna sequence.
@@ -1833,7 +1893,7 @@ public class AlignmentUtils
    * @param seqMappings
    *          the set of mappings involving dnaSeq
    * @param aMapping
-   *          an initial candidate from seqMappings
+   *          a transcript-to-peptide mapping
    * @return
    */
   static SequenceI findCdsForProtein(List<AlignedCodonFrame> mappings,
@@ -1858,7 +1918,15 @@ public class AlignmentUtils
     if (mappedFromLength == dnaLength
             || mappedFromLength == dnaLength - CODON_LENGTH)
     {
-      return seqDss;
+      /*
+       * if sequence has CDS features, this is a transcript with no UTR
+       * - do not take this as the CDS sequence! (JAL-2789)
+       */
+      if (seqDss.getFeatures().getFeaturesByOntology(SequenceOntologyI.CDS)
+              .isEmpty())
+      {
+        return seqDss;
+      }
     }
 
     /*
@@ -1883,10 +1951,12 @@ public class AlignmentUtils
           {
             /*
             * found a 3:1 mapping to the protein product which covers
-            * the whole dna sequence i.e. is from CDS; finally check it
-            * is from the dna start sequence
+            * the whole dna sequence i.e. is from CDS; finally check the CDS
+            * is mapped from the given dna start sequence
             */
             SequenceI cdsSeq = map.getFromSeq();
+            // todo this test is weak if seqMappings contains multiple mappings;
+            // we get away with it if transcript:cds relationship is 1:1
             List<AlignedCodonFrame> dnaToCdsMaps = MappingUtils
                     .findMappingsForSequence(cdsSeq, seqMappings);
             if (!dnaToCdsMaps.isEmpty())
@@ -1916,39 +1986,61 @@ public class AlignmentUtils
   static SequenceI makeCdsSequence(SequenceI seq, Mapping mapping,
           AlignmentI dataset)
   {
-    char[] seqChars = seq.getSequence();
-    List<int[]> fromRanges = mapping.getMap().getFromRanges();
-    int cdsWidth = MappingUtils.getLength(fromRanges);
-    char[] newSeqChars = new char[cdsWidth];
+    /*
+     * construct CDS sequence name as "CDS|" with 'from id' held in the mapping
+     * if set (e.g. EMBL protein_id), else sequence name appended
+     */
+    String mapFromId = mapping.getMappedFromId();
+    final String seqId = "CDS|"
+            + (mapFromId != null ? mapFromId : seq.getName());
 
-    int newPos = 0;
-    for (int[] range : fromRanges)
+    SequenceI newSeq = null;
+
+    final MapList maplist = mapping.getMap();
+    if (maplist.isContiguous() && maplist.isFromForwardStrand())
     {
-      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
+      /*
+       * 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];
+
+      int newPos = 0;
+      for (int[] range : fromRanges)
       {
-        // reverse strand mapping - copy and complement one by one
-        for (int i = range[0]; i >= range[1]; i--)
+        if (range[0] <= range[1])
         {
-          newSeqChars[newPos++] = Dna.getComplement(seqChars[i - 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
+        {
+          // 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]);
+          }
         }
       }
+
+      newSeq = new Sequence(seqId, newSeqChars, 1, newPos);
     }
 
-    /*
-     * assign 'from id' held in the mapping if set (e.g. EMBL protein_id),
-     * else generate a sequence name
-     */
-    String mapFromId = mapping.getMappedFromId();
-    String seqId = "CDS|" + (mapFromId != null ? mapFromId : seq.getName());
-    SequenceI newSeq = new Sequence(seqId, newSeqChars, 1, newPos);
     if (dataset != null)
     {
       SequenceI[] matches = dataset.findSequenceMatch(newSeq.getName());
@@ -1989,21 +2081,23 @@ public class AlignmentUtils
   }
 
   /**
-   * add any DBRefEntrys to cdsSeq from contig that have a Mapping congruent to
+   * Adds any DBRefEntrys to cdsSeq from contig that have a Mapping congruent to
    * the given mapping.
    * 
    * @param cdsSeq
    * @param contig
+   * @param proteinProduct
    * @param mapping
-   * @return list of DBRefEntrys added.
+   * @return list of DBRefEntrys added
    */
-  public static List<DBRefEntry> propagateDBRefsToCDS(SequenceI cdsSeq,
+  protected static List<DBRefEntry> propagateDBRefsToCDS(SequenceI cdsSeq,
           SequenceI contig, SequenceI proteinProduct, Mapping mapping)
   {
 
-    // gather direct refs from contig congrent with mapping
+    // gather direct refs from contig congruent with mapping
     List<DBRefEntry> direct = new ArrayList<>();
     HashSet<String> directSources = new HashSet<>();
+
     if (contig.getDBRefs() != null)
     {
       for (DBRefEntry dbr : contig.getDBRefs())
@@ -2081,7 +2175,7 @@ public class AlignmentUtils
    *          subtypes in the Sequence Ontology)
    * @param omitting
    */
-  public static int transferFeatures(SequenceI fromSeq, SequenceI toSeq,
+  protected static int transferFeatures(SequenceI fromSeq, SequenceI toSeq,
           MapList mapping, String select, String... omitting)
   {
     SequenceI copyTo = toSeq;
@@ -2089,6 +2183,10 @@ public class AlignmentUtils
     {
       copyTo = copyTo.getDatasetSequence();
     }
+    if (fromSeq == copyTo || fromSeq.getDatasetSequence() == copyTo)
+    {
+      return 0; // shared dataset sequence
+    }
 
     /*
      * get features, optionally restricted by an ontology term
@@ -2180,12 +2278,13 @@ public class AlignmentUtils
     int mappedDnaLength = MappingUtils.getLength(ranges);
 
     /*
-     * if not a whole number of codons, something is wrong,
-     * abort mapping
+     * if not a whole number of codons, truncate mapping
      */
-    if (mappedDnaLength % CODON_LENGTH > 0)
+    int codonRemainder = mappedDnaLength % CODON_LENGTH;
+    if (codonRemainder > 0)
     {
-      return null;
+      mappedDnaLength -= codonRemainder;
+      MappingUtils.removeEndPositions(codonRemainder, ranges);
     }
 
     int proteinLength = proteinSeq.getLength();
@@ -2235,7 +2334,7 @@ public class AlignmentUtils
    * @param dnaSeq
    * @return
    */
-  public static List<int[]> findCdsPositions(SequenceI dnaSeq)
+  protected static List<int[]> findCdsPositions(SequenceI dnaSeq)
   {
     List<int[]> result = new ArrayList<>();
 
@@ -2357,7 +2456,8 @@ public class AlignmentUtils
   static int computePeptideVariants(SequenceI peptide, int peptidePos,
           List<DnaVariant>[] codonVariants)
   {
-    String residue = String.valueOf(peptide.getCharAt(peptidePos - 1));
+    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;
@@ -2370,15 +2470,22 @@ public class AlignmentUtils
     {
       if (var.variant != null)
       {
-        String alleles = (String) var.variant.getValue("alleles");
+        String alleles = (String) var.variant.getValue(Gff3Helper.ALLELES);
         if (alleles != null)
         {
           for (String base : alleles.split(","))
           {
-            String codon = base + base2 + base3;
-            if (addPeptideVariant(peptide, peptidePos, residue, var, codon))
+            if (!base1.equalsIgnoreCase(base))
             {
-              count++;
+              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++;
+              }
             }
           }
         }
@@ -2392,15 +2499,22 @@ public class AlignmentUtils
     {
       if (var.variant != null)
       {
-        String alleles = (String) var.variant.getValue("alleles");
+        String alleles = (String) var.variant.getValue(Gff3Helper.ALLELES);
         if (alleles != null)
         {
           for (String base : alleles.split(","))
           {
-            String codon = base1 + base + base3;
-            if (addPeptideVariant(peptide, peptidePos, residue, var, codon))
+            if (!base2.equalsIgnoreCase(base))
             {
-              count++;
+              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++;
+              }
             }
           }
         }
@@ -2414,15 +2528,22 @@ public class AlignmentUtils
     {
       if (var.variant != null)
       {
-        String alleles = (String) var.variant.getValue("alleles");
+        String alleles = (String) var.variant.getValue(Gff3Helper.ALLELES);
         if (alleles != null)
         {
           for (String base : alleles.split(","))
           {
-            String codon = base1 + base2 + base;
-            if (addPeptideVariant(peptide, peptidePos, residue, var, codon))
+            if (!base3.equalsIgnoreCase(base))
             {
-              count++;
+              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++;
+              }
             }
           }
         }
@@ -2433,20 +2554,22 @@ public class AlignmentUtils
   }
 
   /**
-   * Helper method that adds a peptide variant feature, provided the given codon
-   * translates to a value different to the current residue (is a non-synonymous
-   * variant). ID and clinical_significance attributes of the dna variant (if
-   * present) are copied to the new feature.
+   * 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 residue, DnaVariant var, String codon, String canonical)
   {
     /*
      * get peptide translation of codon e.g. GAT -> D
@@ -2454,62 +2577,79 @@ public class AlignmentUtils
      * e.g. multibase variants or HGMD_MUTATION etc
      * are currently ignored here
      */
-    String trans = codon.contains("-") ? "-"
+    String trans = codon.contains("-") ? null
             : (codon.length() > CODON_LENGTH ? null
                     : ResidueProperties.codonTranslate(codon));
-    if (trans != null && !trans.equals(residue))
+    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));
-      String desc = "p." + residue3Char + peptidePos + trans3Char;
-      SequenceFeature sf = new SequenceFeature(
-              SequenceOntologyI.SEQUENCE_VARIANT, desc, peptidePos,
-              peptidePos, var.getSource());
-      StringBuilder attributes = new StringBuilder(32);
-      String id = (String) var.variant.getValue(ID);
-      if (id != null)
-      {
-        if (id.startsWith(SEQUENCE_VARIANT))
-        {
-          id = id.substring(SEQUENCE_VARIANT.length());
-        }
-        sf.setValue(ID, id);
-        attributes.append(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)
+      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))
       {
-        sf.setValue(CLINICAL_SIGNIFICANCE, clinSig);
-        attributes.append(";").append(CLINICAL_SIGNIFICANCE).append("=")
-                .append(clinSig);
+        id = id.substring(SEQUENCE_VARIANT.length());
       }
-      peptide.addSequenceFeature(sf);
-      if (attributes.length() > 0)
+      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
       {
-        sf.setAttributes(attributes.toString());
+        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
       }
-      return true;
     }
-    return false;
+    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
+   * 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
@@ -2547,6 +2687,30 @@ public class AlignmentUtils
         // 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)
       {
@@ -2565,21 +2729,6 @@ public class AlignmentUtils
       }
 
       /*
-       * extract dna variants to a string array
-       */
-      String alls = (String) sf.getValue("alleles");
-      if (alls == null)
-      {
-        continue;
-      }
-      String[] alleles = alls.toUpperCase().split(",");
-      int i = 0;
-      for (String allele : alleles)
-      {
-        alleles[i++] = allele.trim(); // lose any space characters "A, G"
-      }
-
-      /*
        * get this peptide's codon positions e.g. [3, 4, 5] or [4, 7, 10]
        */
       int[] codon = peptidePosition == lastPeptidePostion ? lastCodon
@@ -2688,17 +2837,6 @@ 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);
@@ -2759,12 +2897,14 @@ public class AlignmentUtils
    * true; else returns false
    * 
    * @param unaligned
-   *          - sequences to be aligned based on aligned
+   *                    - sequences to be aligned based on aligned
    * @param aligned
-   *          - 'guide' alignment containing sequences derived from same dataset
-   *          as unaligned
+   *                    - 'guide' alignment containing sequences derived from same
+   *                    dataset as unaligned
    * @return
+   * @deprecated probably obsolete and incomplete
    */
+  @Deprecated
   static boolean alignAsSameSequences(AlignmentI unaligned,
           AlignmentI aligned)
   {
@@ -2786,15 +2926,22 @@ public class AlignmentUtils
     }
 
     /*
-     * first pass - check whether all sequences to be aligned share a dataset
-     * sequence with an aligned sequence
+     * first pass - check whether all sequences to be aligned share a 
+     * dataset sequence with an aligned sequence; also note the leftmost
+     * ungapped column from which to copy
      */
+    int leftmost = Integer.MAX_VALUE;
     for (SequenceI seq : unaligned.getSequences())
     {
-      if (!alignedDatasets.containsKey(seq.getDatasetSequence()))
+      final SequenceI ds = seq.getDatasetSequence();
+      if (!alignedDatasets.containsKey(ds))
       {
         return false;
       }
+      SequenceI alignedSeq = alignedDatasets.get(ds)
+              .get(0);
+      int startCol = alignedSeq.findIndex(seq.getStart()); // 1..
+      leftmost = Math.min(leftmost, startCol);
     }
 
     /*
@@ -2802,13 +2949,25 @@ public class AlignmentUtils
      * heuristic rule: pair off sequences in order for the case where 
      * more than one shares the same dataset sequence 
      */
+    final char gapCharacter = aligned.getGapCharacter();
     for (SequenceI seq : unaligned.getSequences())
     {
       List<SequenceI> alignedSequences = alignedDatasets
               .get(seq.getDatasetSequence());
-      // TODO: getSequenceAsString() will be deprecated in the future
-      // TODO: need to leave to SequenceI implementor to update gaps
-      seq.setSequence(alignedSequences.get(0).getSequenceAsString());
+      SequenceI alignedSeq = alignedSequences.get(0);
+
+      /*
+       * gap fill for leading (5') UTR if any
+       */
+      // TODO this copies intron columns - wrong!
+      int startCol = alignedSeq.findIndex(seq.getStart()); // 1..
+      int endCol = alignedSeq.findIndex(seq.getEnd());
+      char[] seqchars = new char[endCol - leftmost + 1];
+      Arrays.fill(seqchars, gapCharacter);
+      char[] toCopy = alignedSeq.getSequence(startCol - 1, endCol);
+      System.arraycopy(toCopy, 0, seqchars, startCol - leftmost,
+              toCopy.length);
+      seq.setSequence(String.valueOf(seqchars));
       if (alignedSequences.size() > 0)
       {
         // pop off aligned sequences (except the last one)