JAL-2738 use GeneLocus extends DBRefEntry to hold chromosomal mappings
[jalview.git] / src / jalview / analysis / AlignmentUtils.java
index 6acac01..0bc8180 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;
@@ -73,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
@@ -127,7 +131,7 @@ public class AlignmentUtils
    */
   public static AlignmentI expandContext(AlignmentI core, int flankSize)
   {
-    List<SequenceI> sq = new ArrayList<SequenceI>();
+    List<SequenceI> sq = new ArrayList<>();
     int maxoffset = 0;
     for (SequenceI s : core.getSequences())
     {
@@ -257,7 +261,7 @@ public class AlignmentUtils
   public static Map<String, List<SequenceI>> getSequencesByName(
           AlignmentI al)
   {
-    Map<String, List<SequenceI>> theMap = new LinkedHashMap<String, List<SequenceI>>();
+    Map<String, List<SequenceI>> theMap = new LinkedHashMap<>();
     for (SequenceI seq : al.getSequences())
     {
       String name = seq.getName();
@@ -266,7 +270,7 @@ public class AlignmentUtils
         List<SequenceI> seqs = theMap.get(name);
         if (seqs == null)
         {
-          seqs = new ArrayList<SequenceI>();
+          seqs = new ArrayList<>();
           theMap.put(name, seqs);
         }
         seqs.add(seq);
@@ -293,8 +297,8 @@ public class AlignmentUtils
       return false;
     }
 
-    Set<SequenceI> mappedDna = new HashSet<SequenceI>();
-    Set<SequenceI> mappedProtein = new HashSet<SequenceI>();
+    Set<SequenceI> mappedDna = new HashSet<>();
+    Set<SequenceI> mappedProtein = new HashSet<>();
 
     /*
      * First pass - map sequences where cross-references exist. This include
@@ -394,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)
@@ -464,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))
         {
@@ -535,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;
       }
@@ -567,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;
       }
@@ -880,7 +886,7 @@ public class AlignmentUtils
       System.err.println("Wrong alignment type in alignProteinAsDna");
       return 0;
     }
-    List<SequenceI> unmappedProtein = new ArrayList<SequenceI>();
+    List<SequenceI> unmappedProtein = new ArrayList<>();
     Map<AlignedCodon, Map<SequenceI, AlignedCodon>> alignedCodons = buildCodonColumnsMap(
             protein, dna, unmappedProtein);
     return alignProteinAs(protein, alignedCodons, unmappedProtein);
@@ -1091,7 +1097,7 @@ public class AlignmentUtils
      * {dnaSequence, {proteinSequence, codonProduct}} at that position. The
      * comparator keeps the codon positions ordered.
      */
-    Map<AlignedCodon, Map<SequenceI, AlignedCodon>> alignedCodons = new TreeMap<AlignedCodon, Map<SequenceI, AlignedCodon>>(
+    Map<AlignedCodon, Map<SequenceI, AlignedCodon>> alignedCodons = new TreeMap<>(
             new CodonComparator());
 
     for (SequenceI dnaSeq : dna.getSequences())
@@ -1137,9 +1143,9 @@ public class AlignmentUtils
     // TODO delete this ugly hack once JAL-2022 is resolved
     // i.e. we can model startPhase > 0 (incomplete start codon)
 
-    List<SequenceI> sequencesChecked = new ArrayList<SequenceI>();
+    List<SequenceI> sequencesChecked = new ArrayList<>();
     AlignedCodon lastCodon = null;
-    Map<SequenceI, AlignedCodon> toAdd = new HashMap<SequenceI, AlignedCodon>();
+    Map<SequenceI, AlignedCodon> toAdd = new HashMap<>();
 
     for (Entry<AlignedCodon, Map<SequenceI, AlignedCodon>> entry : alignedCodons
             .entrySet())
@@ -1318,7 +1324,7 @@ public class AlignmentUtils
     Map<SequenceI, AlignedCodon> seqProduct = alignedCodons.get(codon);
     if (seqProduct == null)
     {
-      seqProduct = new HashMap<SequenceI, AlignedCodon>();
+      seqProduct = new HashMap<>();
       alignedCodons.put(codon, seqProduct);
     }
     seqProduct.put(protein, codon);
@@ -1455,7 +1461,7 @@ public class AlignmentUtils
       {
         continue;
       }
-      final List<AlignmentAnnotation> result = new ArrayList<AlignmentAnnotation>();
+      final List<AlignmentAnnotation> result = new ArrayList<>();
       for (AlignmentAnnotation dsann : datasetAnnotations)
       {
         /*
@@ -1637,17 +1643,17 @@ public class AlignmentUtils
       throw new IllegalArgumentException(
               "IMPLEMENTATION ERROR: dataset.getDataset() must be null!");
     }
-    List<SequenceI> foundSeqs = new ArrayList<SequenceI>();
-    List<SequenceI> cdsSeqs = new ArrayList<SequenceI>();
+    List<SequenceI> foundSeqs = new ArrayList<>();
+    List<SequenceI> cdsSeqs = new ArrayList<>();
     List<AlignedCodonFrame> mappings = dataset.getCodonFrames();
     HashSet<SequenceI> productSeqs = null;
     if (products != null)
     {
-      productSeqs = new HashSet<SequenceI>();
+      productSeqs = new HashSet<>();
       for (SequenceI seq : products)
       {
-        productSeqs.add(seq.getDatasetSequence() == null ? seq
-                : seq.getDatasetSequence());
+        productSeqs.add(seq.getDatasetSequence() == null ? seq : seq
+                .getDatasetSequence());
       }
     }
 
@@ -1740,9 +1746,8 @@ public class AlignmentUtils
           /*
            * add a mapping from CDS to the (unchanged) mapped to range
            */
-          List<int[]> cdsRange = Collections
-                  .singletonList(new int[]
-                  { 1, cdsSeq.getLength() });
+          List<int[]> cdsRange = Collections.singletonList(new int[] { 1,
+              cdsSeq.getLength() });
           MapList cdsToProteinMap = new MapList(cdsRange,
                   mapList.getToRanges(), mapList.getFromRatio(),
                   mapList.getToRatio());
@@ -1764,7 +1769,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);
@@ -1774,6 +1779,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
@@ -1792,26 +1804,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);
           }
 
@@ -1824,14 +1840,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.
@@ -1843,7 +1891,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,
@@ -1868,7 +1916,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;
+      }
     }
 
     /*
@@ -1893,10 +1949,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())
@@ -1999,21 +2057,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
-    List<DBRefEntry> direct = new ArrayList<DBRefEntry>();
-    HashSet<String> directSources = new HashSet<String>();
+    // 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())
@@ -2033,7 +2093,7 @@ public class AlignmentUtils
     DBRefEntry[] onSource = DBRefUtils.selectRefs(
             proteinProduct.getDBRefs(),
             directSources.toArray(new String[0]));
-    List<DBRefEntry> propagated = new ArrayList<DBRefEntry>();
+    List<DBRefEntry> propagated = new ArrayList<>();
 
     // and generate appropriate mappings
     for (DBRefEntry cdsref : direct)
@@ -2091,7 +2151,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;
@@ -2174,7 +2234,10 @@ public class AlignmentUtils
 
   /**
    * Returns a mapping from dna to protein by inspecting sequence features of
-   * type "CDS" on the dna.
+   * type "CDS" on the dna. A mapping is constructed if the total CDS feature
+   * length is 3 times the peptide length (optionally after dropping a trailing
+   * stop codon). This method does not check whether the CDS nucleotide sequence
+   * translates to the peptide sequence.
    * 
    * @param dnaSeq
    * @param proteinSeq
@@ -2186,6 +2249,16 @@ public class AlignmentUtils
     List<int[]> ranges = findCdsPositions(dnaSeq);
     int mappedDnaLength = MappingUtils.getLength(ranges);
 
+    /*
+     * if not a whole number of codons, truncate mapping
+     */
+    int codonRemainder = mappedDnaLength % CODON_LENGTH;
+    if (codonRemainder > 0)
+    {
+      mappedDnaLength -= codonRemainder;
+      MappingUtils.removeEndPositions(codonRemainder, ranges);
+    }
+
     int proteinLength = proteinSeq.getLength();
     int proteinStart = proteinSeq.getStart();
     int proteinEnd = proteinSeq.getEnd();
@@ -2200,7 +2273,7 @@ public class AlignmentUtils
       proteinStart++;
       proteinLength--;
     }
-    List<int[]> proteinRange = new ArrayList<int[]>();
+    List<int[]> proteinRange = new ArrayList<>();
 
     /*
      * dna length should map to protein (or protein plus stop codon)
@@ -2209,8 +2282,12 @@ public class AlignmentUtils
     if (codesForResidues == (proteinLength + 1))
     {
       // assuming extra codon is for STOP and not in peptide
+      // todo: check trailing codon is indeed a STOP codon
       codesForResidues--;
+      mappedDnaLength -= CODON_LENGTH;
+      MappingUtils.removeEndPositions(CODON_LENGTH, ranges);
     }
+
     if (codesForResidues == proteinLength)
     {
       proteinRange.add(new int[] { proteinStart, proteinEnd });
@@ -2221,7 +2298,7 @@ public class AlignmentUtils
 
   /**
    * Returns a list of CDS ranges found (as sequence positions base 1), i.e. of
-   * start/end positions of sequence features of type "CDS" (or a sub-type of
+   * [start, end] positions of sequence features of type "CDS" (or a sub-type of
    * CDS in the Sequence Ontology). The ranges are sorted into ascending start
    * position order, so this method is only valid for linear CDS in the same
    * sense as the protein product.
@@ -2229,9 +2306,9 @@ 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<int[]>();
+    List<int[]> result = new ArrayList<>();
 
     List<SequenceFeature> sfs = dnaSeq.getFeatures().getFeaturesByOntology(
             SequenceOntologyI.CDS);
@@ -2240,7 +2317,6 @@ public class AlignmentUtils
       return result;
     }
     SequenceFeatures.sortFeatures(sfs, true);
-    int startPhase = 0;
 
     for (SequenceFeature sf : sfs)
     {
@@ -2258,7 +2334,7 @@ public class AlignmentUtils
        */
       int begin = sf.getBegin();
       int end = sf.getEnd();
-      if (result.isEmpty())
+      if (result.isEmpty() && phase > 0)
       {
         begin += phase;
         if (begin > end)
@@ -2273,16 +2349,6 @@ public class AlignmentUtils
     }
 
     /*
-     * remove 'startPhase' positions (usually 0) from the first range 
-     * so we begin at the start of a complete codon
-     */
-    if (!result.isEmpty())
-    {
-      // TODO JAL-2022 correctly model start phase > 0
-      result.get(0)[0] += startPhase;
-    }
-
-    /*
      * Finally sort ranges by start position. This avoids a dependency on 
      * keeping features in order on the sequence (if they are in order anyway,
      * the sort will have almost no work to do). The implicit assumption is CDS
@@ -2380,10 +2446,17 @@ public class AlignmentUtils
         {
           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++;
+              }
             }
           }
         }
@@ -2402,10 +2475,17 @@ public class AlignmentUtils
         {
           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++;
+              }
             }
           }
         }
@@ -2424,10 +2504,17 @@ public class AlignmentUtils
         {
           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++;
+              }
             }
           }
         }
@@ -2438,20 +2525,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
@@ -2459,62 +2548,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
@@ -2528,7 +2634,7 @@ public class AlignmentUtils
      * 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<Integer, List<DnaVariant>[]>();
+    LinkedHashMap<Integer, List<DnaVariant>[]> variants = new LinkedHashMap<>();
 
     List<SequenceFeature> dnaFeatures = dnaSeq.getFeatures()
             .getFeaturesByOntology(SequenceOntologyI.SEQUENCE_VARIANT);
@@ -2554,18 +2660,26 @@ public class AlignmentUtils
       }
 
       /*
-       * extract dna variants to a string array
+       * 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(",");
-      int i = 0;
+      boolean isSnp = true;
       for (String allele : alleles)
       {
-        alleles[i++] = allele.trim(); // lose any space characters "A, G"
+        if (allele.trim().length() > 1)
+        {
+          isSnp = false;
+        }
+      }
+      if (!isSnp)
+      {
+        continue;
       }
 
       int[] mapsTo = dnaToProtein.locateInTo(dnaCol, dnaCol);
@@ -2579,9 +2693,9 @@ public class AlignmentUtils
       if (codonVariants == null)
       {
         codonVariants = new ArrayList[CODON_LENGTH];
-        codonVariants[0] = new ArrayList<DnaVariant>();
-        codonVariants[1] = new ArrayList<DnaVariant>();
-        codonVariants[2] = new ArrayList<DnaVariant>();
+        codonVariants[0] = new ArrayList<>();
+        codonVariants[1] = new ArrayList<>();
+        codonVariants[2] = new ArrayList<>();
         variants.put(peptidePosition, codonVariants);
       }
 
@@ -2705,7 +2819,7 @@ public class AlignmentUtils
     /*
      * fancy case - aligning via mappings between sequences
      */
-    List<SequenceI> unmapped = new ArrayList<SequenceI>();
+    List<SequenceI> unmapped = new ArrayList<>();
     Map<Integer, Map<SequenceI, Character>> columnMap = buildMappedColumnsMap(
             unaligned, aligned, unmapped);
     int width = columnMap.size();
@@ -2780,7 +2894,7 @@ public class AlignmentUtils
     }
 
     // map from dataset sequence to alignment sequence(s)
-    Map<SequenceI, List<SequenceI>> alignedDatasets = new HashMap<SequenceI, List<SequenceI>>();
+    Map<SequenceI, List<SequenceI>> alignedDatasets = new HashMap<>();
     for (SequenceI seq : aligned.getSequences())
     {
       SequenceI ds = seq.getDatasetSequence();
@@ -2843,7 +2957,7 @@ public class AlignmentUtils
      * {unalignedSequence, characterPerSequence} at that position.
      * TreeMap keeps the entries in ascending column order. 
      */
-    SortedMap<Integer, Map<SequenceI, Character>> map = new TreeMap<Integer, Map<SequenceI, Character>>();
+    SortedMap<Integer, Map<SequenceI, Character>> map = new TreeMap<>();
 
     /*
      * record any sequences that have no mapping so can't be realigned
@@ -2948,7 +3062,7 @@ public class AlignmentUtils
             Map<SequenceI, Character> seqsMap = map.get(fromCol);
             if (seqsMap == null)
             {
-              seqsMap = new HashMap<SequenceI, Character>();
+              seqsMap = new HashMap<>();
               map.put(fromCol, seqsMap);
             }
             seqsMap.put(seq, seq.getCharAt(mappedCharPos - toStart));