JAL-2110 ensuring mapped dbrefs between protein and cds
authorgmungoc <g.m.carstairs@dundee.ac.uk>
Wed, 6 Jul 2016 12:08:28 +0000 (13:08 +0100)
committergmungoc <g.m.carstairs@dundee.ac.uk>
Wed, 6 Jul 2016 12:08:28 +0000 (13:08 +0100)
src/jalview/analysis/AlignmentUtils.java
src/jalview/analysis/CrossRef.java
src/jalview/datamodel/AlignedCodonFrame.java
src/jalview/datamodel/Alignment.java
src/jalview/ext/ensembl/EnsemblSeqProxy.java
test/jalview/analysis/AlignmentUtilsTests.java

index aa7cb18..35aa91d 100644 (file)
@@ -22,6 +22,7 @@ package jalview.analysis;
 
 import static jalview.io.gff.GffConstants.CLINICAL_SIGNIFICANCE;
 
+import jalview.api.DBRefEntryI;
 import jalview.datamodel.AlignedCodon;
 import jalview.datamodel.AlignedCodonFrame;
 import jalview.datamodel.AlignedCodonFrame.SequenceToSequenceMapping;
@@ -39,6 +40,7 @@ import jalview.io.gff.SequenceOntologyFactory;
 import jalview.io.gff.SequenceOntologyI;
 import jalview.schemes.ResidueProperties;
 import jalview.util.Comparison;
+import jalview.util.DBRefUtils;
 import jalview.util.MapList;
 import jalview.util.MappingUtils;
 import jalview.util.StringUtils;
@@ -850,6 +852,11 @@ public class AlignmentUtils
    */
   public static int alignProteinAsDna(AlignmentI protein, AlignmentI dna)
   {
+    if (protein.isNucleotide() || !dna.isNucleotide())
+    {
+      System.err.println("Wrong alignment type in alignProteinAsDna");
+      return 0;
+    }
     List<SequenceI> unmappedProtein = new ArrayList<SequenceI>();
     Map<AlignedCodon, Map<SequenceI, AlignedCodon>> alignedCodons = buildCodonColumnsMap(
             protein, dna, unmappedProtein);
@@ -857,6 +864,68 @@ public class AlignmentUtils
   }
 
   /**
+   * Realigns the given dna to match the alignment of the protein, using codon
+   * mappings to translate aligned peptide positions to codons.
+   * 
+   * @param dna
+   *          the alignment whose sequences are realigned by this method
+   * @param protein
+   *          the protein alignment whose alignment we are 'copying'
+   * @return the number of sequences that were realigned
+   */
+  public static int alignCdsAsProtein(AlignmentI dna, AlignmentI protein)
+  {
+    if (protein.isNucleotide() || !dna.isNucleotide())
+    {
+      System.err.println("Wrong alignment type in alignProteinAsDna");
+      return 0;
+    }
+    // todo: implement this
+    List<AlignedCodonFrame> mappings = protein.getCodonFrames();
+    int alignedCount = 0;
+    for (SequenceI dnaSeq : dna.getSequences())
+    {
+      if (alignCdsSequenceAsProtein(dnaSeq, protein, mappings))
+      {
+        alignedCount++;
+      }
+    }
+    return alignedCount;
+  }
+
+  /**
+   * Helper method to align (if possible) the dna sequence (cds only) to match
+   * the alignment of a mapped protein sequence
+   * 
+   * @param dnaSeq
+   * @param protein
+   * @param mappings
+   * @return
+   */
+  static boolean alignCdsSequenceAsProtein(SequenceI dnaSeq,
+          AlignmentI protein, List<AlignedCodonFrame> mappings)
+  {
+    List<AlignedCodonFrame> dnaMappings = MappingUtils
+            .findMappingsForSequence(dnaSeq, mappings);
+    for (AlignedCodonFrame mapping : dnaMappings)
+    {
+      SequenceI peptide = mapping.findAlignedSequence(dnaSeq, protein);
+      if (peptide != null)
+      {
+        Mapping map = mapping.getMappingBetween(dnaSeq, peptide);
+        if (map != null)
+        {
+          /*
+           * traverse peptide adding gaps or codons to new cds sequence
+           */
+          MapList mapList = map.getMap();
+        }
+      }
+    }
+    return false;
+  }
+
+  /**
    * Builds a map whose key is an aligned codon position (3 alignment column
    * numbers base 0), and whose value is a map from protein sequence to each
    * protein's peptide residue for that codon. The map generates an ordering of
@@ -1404,9 +1473,9 @@ public class AlignmentUtils
    * added to the alignment dataset.
    * 
    * @param dna
-   *          aligned dna sequences
+   *          aligned nucleotide (dna or cds) sequences
    * @param dataset
-   *          - throws error if not given a dataset
+   *          the alignment dataset the sequences belong to
    * @param products
    *          (optional) to restrict results to CDS that map to specified
    *          protein products
@@ -1414,20 +1483,21 @@ public class AlignmentUtils
    *         sequences (or null if no mappings are found)
    */
   public static AlignmentI makeCdsAlignment(SequenceI[] dna,
-          AlignmentI dataset, AlignmentI products)
+          AlignmentI dataset, SequenceI[] products)
   {
-    if (dataset.getDataset() != null)
+    if (dataset == null || dataset.getDataset() != null)
     {
-      throw new Error(
+      throw new IllegalArgumentException(
               "IMPLEMENTATION ERROR: dataset.getDataset() must be null!");
     }
+    List<SequenceI> foundSeqs = new ArrayList<SequenceI>();
     List<SequenceI> cdsSeqs = new ArrayList<SequenceI>();
     List<AlignedCodonFrame> mappings = dataset.getCodonFrames();
     HashSet<SequenceI> productSeqs = null;
     if (products != null)
     {
       productSeqs = new HashSet<SequenceI>();
-      for (SequenceI seq : products.getSequences())
+      for (SequenceI seq : products)
       {
         productSeqs.add(seq.getDatasetSequence() == null ? seq : seq
                 .getDatasetSequence());
@@ -1435,24 +1505,32 @@ public class AlignmentUtils
     }
 
     /*
-     * construct CDS sequences from the (cds-to-protein) mappings made earlier;
-     * this makes it possible to model multiple products from dna (e.g. EMBL); 
-     * however it does mean we don't have the EMBL protein_id (a property on 
-     * the CDS features) in order to make the CDS sequence name :-( 
+     * Construct CDS sequences from mappings on the alignment dataset.
+     * The logic is:
+     * - find the protein product(s) mapped to from each dna sequence
+     * - if the mapping covers the whole dna sequence (give or take start/stop
+     *   codon), take the dna as the CDS sequence
+     * - else search dataset mappings for a suitable dna sequence, i.e. one
+     *   whose whole sequence is mapped to the protein 
+     * - if no sequence found, construct one from the dna sequence and mapping
+     *   (and add it to dataset so it is found if this is repeated)
      */
-    for (SequenceI seq : dna)
+    for (SequenceI dnaSeq : dna)
     {
-      SequenceI seqDss = seq.getDatasetSequence() == null ? seq : seq
-              .getDatasetSequence();
+      SequenceI dnaDss = dnaSeq.getDatasetSequence() == null ? dnaSeq
+              : dnaSeq.getDatasetSequence();
+
       List<AlignedCodonFrame> seqMappings = MappingUtils
-              .findMappingsForSequence(seq, mappings);
+              .findMappingsForSequence(dnaSeq, mappings);
       for (AlignedCodonFrame mapping : seqMappings)
       {
-        List<Mapping> mappingsFromSequence = mapping.getMappingsFromSequence(seq);
+        List<Mapping> mappingsFromSequence = mapping
+                .getMappingsFromSequence(dnaSeq);
 
         for (Mapping aMapping : mappingsFromSequence)
         {
-          if (aMapping.getMap().getFromRatio() == 1)
+          MapList mapList = aMapping.getMap();
+          if (mapList.getFromRatio() == 1)
           {
             /*
              * not a dna-to-protein mapping (likely dna-to-cds)
@@ -1461,55 +1539,33 @@ public class AlignmentUtils
           }
 
           /*
-           * check for an existing CDS sequence i.e. a 3:1 mapping to 
-           * the dna mapping's product
+           * skip if mapping is not to one of the target set of proteins
            */
-          SequenceI cdsSeq = null;
-
-          // TODO better mappings collection data model so we can do
-          // a direct lookup instead of double loops to find mappings
-
           SequenceI proteinProduct = aMapping.getTo();
-
-          /*
-           * skip if not mapped to one of a specified set of proteins
-           */
           if (productSeqs != null && !productSeqs.contains(proteinProduct))
           {
             continue;
           }
 
-          for (AlignedCodonFrame acf : MappingUtils
-                  .findMappingsForSequence(proteinProduct, mappings))
+          /*
+           * try to locate the CDS from the dataset mappings;
+           * guard against duplicate results (for the case that protein has
+           * dbrefs to both dna and cds sequences)
+           */
+          SequenceI cdsSeq = findCdsForProtein(mappings, dnaSeq,
+                  seqMappings, aMapping);
+          if (cdsSeq != null)
           {
-            for (SequenceToSequenceMapping map : acf.getMappings())
+            if (!foundSeqs.contains(cdsSeq))
             {
-              if (map.getMapping().getMap().getFromRatio() == 3
-                      && proteinProduct == map.getMapping().getTo()
-                      && seqDss != map.getFromSeq())
+              foundSeqs.add(cdsSeq);
+              SequenceI derivedSequence = cdsSeq.deriveSequence();
+              cdsSeqs.add(derivedSequence);
+              if (!dataset.getSequences().contains(cdsSeq))
               {
-                /*
-                 * found a 3:1 mapping to the protein product which is not
-                 * from the dna sequence...assume it is from the CDS sequence
-                 * TODO mappings data model that brings together related
-                 * dna-cds-protein mappings in one object
-                 */
-                cdsSeq = map.getFromSeq();
+                dataset.addSequence(cdsSeq);
               }
             }
-          }
-          if (cdsSeq != null)
-          {
-            /*
-             * mappings are always to dataset sequences so create an aligned
-             * sequence to own it; add the dataset sequence to the dataset
-             */
-            SequenceI derivedSequence = cdsSeq.deriveSequence();
-            cdsSeqs.add(derivedSequence);
-            if (!dataset.getSequences().contains(cdsSeq))
-            {
-              dataset.addSequence(cdsSeq);
-            }
             continue;
           }
 
@@ -1517,7 +1573,7 @@ public class AlignmentUtils
            * didn't find mapped CDS sequence - construct it and add
            * its dataset sequence to the dataset
            */
-          cdsSeq = makeCdsSequence(seq.getDatasetSequence(), aMapping);
+          cdsSeq = makeCdsSequence(dnaSeq.getDatasetSequence(), aMapping);
           SequenceI cdsSeqDss = cdsSeq.createDatasetSequence();
           cdsSeqs.add(cdsSeq);
           if (!dataset.getSequences().contains(cdsSeqDss))
@@ -1530,11 +1586,10 @@ public class AlignmentUtils
            */
           List<int[]> cdsRange = Collections.singletonList(new int[] { 1,
               cdsSeq.getLength() });
-          MapList map = new MapList(cdsRange, aMapping.getMap()
-                  .getToRanges(), aMapping.getMap().getFromRatio(),
-                  aMapping.getMap().getToRatio());
+          MapList cdsToProteinMap = new MapList(cdsRange, mapList.getToRanges(),
+                  mapList.getFromRatio(), mapList.getToRatio());
           AlignedCodonFrame cdsToProteinMapping = new AlignedCodonFrame();
-          cdsToProteinMapping.addMap(cdsSeq, proteinProduct, map);
+          cdsToProteinMapping.addMap(cdsSeq, proteinProduct, cdsToProteinMap);
 
           /*
            * guard against duplicating the mapping if repeating this action
@@ -1545,22 +1600,59 @@ public class AlignmentUtils
           }
 
           /*
+           * copy protein's dbrefs to CDS sequence
+           * this enables Get Cross-References from CDS alignment
+           */
+          DBRefEntry[] proteinRefs = DBRefUtils.selectDbRefs(false,
+                  proteinProduct.getDBRefs());
+          if (proteinRefs != null)
+          {
+            for (DBRefEntry ref : proteinRefs)
+            {
+              DBRefEntry cdsToProteinRef = new DBRefEntry(ref);
+              cdsToProteinRef.setMap(new Mapping(proteinProduct,
+                      cdsToProteinMap));
+              cdsSeqDss.addDBRef(cdsToProteinRef);
+            }
+          }
+
+          /*
            * add another mapping from original 'from' range to CDS
            */
-          AlignedCodonFrame dnaToProteinMapping = new AlignedCodonFrame();
-          map = new MapList(aMapping.getMap().getFromRanges(), cdsRange, 1,
+          AlignedCodonFrame dnaToCdsMapping = new AlignedCodonFrame();
+          MapList dnaToCdsMap = new MapList(mapList.getFromRanges(),
+                  cdsRange, 1,
                   1);
-          dnaToProteinMapping.addMap(seq.getDatasetSequence(), cdsSeq, map);
-          if (!mappings.contains(dnaToProteinMapping))
+          dnaToCdsMapping.addMap(dnaSeq.getDatasetSequence(), cdsSeq,
+                  dnaToCdsMap);
+          if (!mappings.contains(dnaToCdsMapping))
           {
-            mappings.add(dnaToProteinMapping);
+            mappings.add(dnaToCdsMapping);
           }
 
+          /*
+           * 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
+           * same source and accession, so need a different accession for
+           * the CDS from the dna sequence
+           */
+          DBRefEntryI dnaRef = dnaDss.getSourceDBRef();
+          if (dnaRef != null)
+          {
+            // assuming cds version same as dna ?!?
+            DBRefEntry proteinToCdsRef = new DBRefEntry(dnaRef.getSource(),
+                    dnaRef.getVersion(), cdsSeq.getName());
+            proteinToCdsRef.setMap(new Mapping(cdsSeqDss, cdsToProteinMap
+                    .getInverse()));
+            proteinProduct.addDBRef(proteinToCdsRef);
+          }
 
           /*
            * transfer any features on dna that overlap the CDS
            */
-          transferFeatures(seq, cdsSeq, map, null, SequenceOntologyI.CDS);
+          transferFeatures(dnaSeq, cdsSeq, cdsToProteinMap, null,
+                  SequenceOntologyI.CDS);
         }
       }
     }
@@ -1573,6 +1665,83 @@ public class AlignmentUtils
   }
 
   /**
+   * 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.
+   * 
+   * @param mappings
+   *          set of all mappings on the dataset
+   * @param dnaSeq
+   *          a dna (or cds) sequence we are searching from
+   * @param seqMappings
+   *          the set of mappings involving dnaSeq
+   * @param aMapping
+   *          an initial candidate from seqMappings
+   * @return
+   */
+  static SequenceI findCdsForProtein(List<AlignedCodonFrame> mappings,
+          SequenceI dnaSeq, List<AlignedCodonFrame> seqMappings,
+          Mapping aMapping)
+  {
+    /*
+     * TODO a better dna-cds-protein mapping data representation to allow easy
+     * navigation; until then this clunky looping around lists of mappings
+     */
+    SequenceI seqDss = dnaSeq.getDatasetSequence() == null ? dnaSeq
+            : dnaSeq.getDatasetSequence();
+    SequenceI proteinProduct = aMapping.getTo();
+
+    /*
+     * is this mapping from the whole dna sequence (i.e. CDS)?
+     * allowing for possible stop codon on dna but not peptide
+     */
+    int mappedFromLength = MappingUtils.getLength(aMapping.getMap()
+            .getFromRanges());
+    int dnaLength = seqDss.getLength();
+    if (mappedFromLength == dnaLength || mappedFromLength == dnaLength - 3)
+    {
+      return seqDss;
+    }
+
+    /*
+     * looks like we found the dna-to-protein mapping; search for the
+     * corresponding cds-to-protein mapping
+     */
+    List<AlignedCodonFrame> mappingsToPeptide = MappingUtils
+            .findMappingsForSequence(proteinProduct, mappings);
+    for (AlignedCodonFrame acf : mappingsToPeptide)
+    {
+      for (SequenceToSequenceMapping map : acf.getMappings())
+      {
+        Mapping mapping = map.getMapping();
+        if (mapping != aMapping && mapping.getMap().getFromRatio() == 3
+                && proteinProduct == mapping.getTo()
+                && seqDss != map.getFromSeq())
+        {
+          mappedFromLength = MappingUtils.getLength(mapping.getMap()
+                  .getFromRanges());
+          if (mappedFromLength == map.getFromSeq().getLength())
+          {
+            /*
+            * 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
+            */
+            SequenceI cdsSeq = map.getFromSeq();
+            List<AlignedCodonFrame> dnaToCdsMaps = MappingUtils
+                    .findMappingsForSequence(cdsSeq, seqMappings);
+            if (!dnaToCdsMaps.isEmpty())
+            {
+              return cdsSeq;
+            }
+          }
+        }
+      }
+    }
+    return null;
+  }
+
+  /**
    * Helper method that makes a CDS sequence as defined by the mappings from the
    * given sequence i.e. extracts the 'mapped from' ranges (which may be on
    * forward or reverse strand).
@@ -1609,8 +1778,15 @@ public class AlignmentUtils
       }
     }
 
-    SequenceI newSeq = new Sequence(seq.getName() + "|"
-            + mapping.getTo().getName(), 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);
+    // newSeq.setDescription(mapFromId);
+
     return newSeq;
   }
 
index 2b5a0e2..df5ca10 100644 (file)
@@ -295,12 +295,14 @@ public class CrossRef
               if (fromDna)
               {
                 // map is from dna seq to a protein product
-                cf.addMap(dss, rsq, xref.getMap().getMap());
+                cf.addMap(dss, rsq, xref.getMap().getMap(), xref.getMap()
+                        .getMappedFromId());
               }
               else
               {
                 // map should be from protein seq to its coding dna
-                cf.addMap(rsq, dss, xref.getMap().getMap().getInverse());
+                cf.addMap(rsq, dss, xref.getMap().getMap().getInverse(),
+                        xref.getMap().getMappedFromId());
               }
             }
           }
@@ -621,14 +623,14 @@ public class CrossRef
   void updateDbrefMappings(SequenceI mapFrom, DBRefEntry[] xrefs,
           SequenceI[] retrieved, AlignedCodonFrame acf, boolean fromDna)
   {
-    SequenceIdMatcher matcher = new SequenceIdMatcher(retrieved);
+    SequenceIdMatcher idMatcher = new SequenceIdMatcher(retrieved);
     for (DBRefEntry xref : xrefs)
     {
       if (!xref.hasMap())
       {
         String targetSeqName = xref.getSource() + "|"
                 + xref.getAccessionId();
-        SequenceI[] matches = matcher.findAllIdMatches(targetSeqName);
+        SequenceI[] matches = idMatcher.findAllIdMatches(targetSeqName);
         if (matches == null)
         {
           return;
@@ -705,6 +707,20 @@ public class CrossRef
       return false;
     }
     xref.setMap(new Mapping(mapTo, mapping));
+
+    /*
+     * and add a reverse DbRef with the inverse mapping
+     */
+    if (mapFrom.getDatasetSequence() != null
+            && mapFrom.getDatasetSequence().getSourceDBRef() != null)
+    {
+      DBRefEntry dbref = new DBRefEntry(mapFrom.getDatasetSequence()
+              .getSourceDBRef());
+      dbref.setMap(new Mapping(mapFrom.getDatasetSequence(), mapping
+              .getInverse()));
+      mapTo.addDBRef(dbref);
+    }
+
     if (fromDna)
     {
       AlignmentUtils.computeProteinFeatures(mapFrom, mapTo, mapping);
index bb705b6..326cc4e 100644 (file)
@@ -128,6 +128,21 @@ public class AlignedCodonFrame
    */
   public void addMap(SequenceI dnaseq, SequenceI aaseq, MapList map)
   {
+    addMap(dnaseq, aaseq, map, null);
+  }
+
+  /**
+   * Adds a mapping between the dataset sequences for the associated dna and
+   * protein sequence objects
+   * 
+   * @param dnaseq
+   * @param aaseq
+   * @param map
+   * @param mapFromId
+   */
+  public void addMap(SequenceI dnaseq, SequenceI aaseq, MapList map,
+          String mapFromId)
+  {
     // JBPNote DEBUG! THIS !
     // dnaseq.transferAnnotation(aaseq, mp);
     // aaseq.transferAnnotation(dnaseq, new Mapping(map.getInverse()));
@@ -155,6 +170,7 @@ public class AlignedCodonFrame
      * otherwise, add a new sequence mapping
      */
     Mapping mp = new Mapping(toSeq, map);
+    mp.setMappedFromId(mapFromId);
     mappings.add(new SequenceToSequenceMapping(fromSeq, mp));
   }
 
index 2b194cc..f5665db 100755 (executable)
@@ -30,7 +30,6 @@ import java.util.Collections;
 import java.util.Enumeration;
 import java.util.HashSet;
 import java.util.Hashtable;
-import java.util.Iterator;
 import java.util.List;
 import java.util.Map;
 import java.util.Set;
@@ -1724,6 +1723,10 @@ public class Alignment implements AlignmentI
     {
       return AlignmentUtils.alignProteinAsDna(this, al);
     }
+    else if (thatIsProtein && thisIsNucleotide)
+    {
+      return AlignmentUtils.alignCdsAsProtein(this, al);
+    }
     return AlignmentUtils.alignAs(this, al);
   }
 
index c86469f..31552af 100644 (file)
@@ -5,7 +5,6 @@ import jalview.analysis.Dna;
 import jalview.datamodel.Alignment;
 import jalview.datamodel.AlignmentI;
 import jalview.datamodel.DBRefEntry;
-import jalview.datamodel.DBRefSource;
 import jalview.datamodel.Mapping;
 import jalview.datamodel.SequenceFeature;
 import jalview.datamodel.SequenceI;
@@ -315,13 +314,6 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
     for (DBRefEntry xref : xrefs)
     {
       seq.addDBRef(xref);
-      /*
-       * Save any Uniprot xref to be the reference for SIFTS mapping
-       */
-      if (DBRefSource.UNIPROT.equals(xref.getSource()))
-      {
-        seq.setSourceDBRef(xref);
-      }
     }
 
     /*
@@ -330,6 +322,7 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
     DBRefEntry self = new DBRefEntry(getDbSource(),
             getEnsemblDataVersion(), seq.getName());
     seq.addDBRef(self);
+    seq.setSourceDBRef(self);
   }
 
   /**
@@ -387,8 +380,9 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
         if (ids.contains(name)
                 || ids.contains(name.replace("ENSP", "ENST")))
         {
-          DBRefUtils.parseToDbRef(sq, getDbSource(),
+          DBRefEntry dbref = DBRefUtils.parseToDbRef(sq, getDbSource(),
                   getEnsemblDataVersion(), name);
+          sq.setSourceDBRef(dbref);
         }
       }
       if (alignment == null)
index d704ec6..727cf22 100644 (file)
@@ -22,6 +22,7 @@ package jalview.analysis;
 
 import static org.testng.AssertJUnit.assertEquals;
 import static org.testng.AssertJUnit.assertFalse;
+import static org.testng.AssertJUnit.assertNotNull;
 import static org.testng.AssertJUnit.assertNull;
 import static org.testng.AssertJUnit.assertSame;
 import static org.testng.AssertJUnit.assertTrue;
@@ -973,10 +974,17 @@ public class AlignmentUtilsTests
   @Test(groups = { "Functional" })
   public void testMakeCdsAlignment()
   {
+    /*
+     * scenario:
+     *     dna1 --> [4, 6] [10,12]        --> pep1 
+     *     dna2 --> [1, 3] [7, 9] [13,15] --> pep1 
+     */
     SequenceI dna1 = new Sequence("dna1", "aaaGGGcccTTTaaa");
     SequenceI dna2 = new Sequence("dna2", "GGGcccTTTaaaCCC");
     SequenceI pep1 = new Sequence("pep1", "GF");
     SequenceI pep2 = new Sequence("pep2", "GFP");
+    pep1.addDBRef(new DBRefEntry("UNIPROT", "0", "pep1"));
+    pep2.addDBRef(new DBRefEntry("UNIPROT", "0", "pep2"));
     dna1.createDatasetSequence();
     dna2.createDatasetSequence();
     pep1.createDatasetSequence();
@@ -984,6 +992,19 @@ public class AlignmentUtilsTests
     AlignmentI dna = new Alignment(new SequenceI[] { dna1, dna2 });
     dna.setDataset(null);
 
+    /*
+     * need a sourceDbRef if we are to construct dbrefs to the CDS
+     * sequence
+     */
+    DBRefEntry dbref = new DBRefEntry("ENSEMBL", "0", "dna1");
+    dna1.getDatasetSequence().setSourceDBRef(dbref);
+    dbref = new DBRefEntry("ENSEMBL", "0", "dna2");
+    dna2.getDatasetSequence().setSourceDBRef(dbref);
+
+    /*
+     * CDS sequences are 'discovered' from dna-to-protein mappings on the alignment
+     * dataset (e.g. added from dbrefs by CrossRef.findXrefSequences)
+     */
     MapList map = new MapList(new int[] { 4, 6, 10, 12 },
             new int[] { 1, 2 }, 3, 1);
     AlignedCodonFrame acf = new AlignedCodonFrame();
@@ -1001,6 +1022,9 @@ public class AlignmentUtilsTests
     AlignmentI cds = AlignmentUtils.makeCdsAlignment(new SequenceI[] {
         dna1, dna2 }, dna.getDataset(), null);
 
+    /*
+     * verify cds sequences
+     */
     assertEquals(2, cds.getSequences().size());
     assertEquals("GGGTTT", cds.getSequenceAt(0).getSequenceAsString());
     assertEquals("GGGTTTCCC", cds.getSequenceAt(1).getSequenceAsString());
@@ -1009,19 +1033,45 @@ public class AlignmentUtilsTests
      * verify shared, extended alignment dataset
      */
     assertSame(dna.getDataset(), cds.getDataset());
-    assertTrue(dna.getDataset().getSequences()
-            .contains(cds.getSequenceAt(0).getDatasetSequence()));
-    assertTrue(dna.getDataset().getSequences()
-            .contains(cds.getSequenceAt(1).getDatasetSequence()));
+    SequenceI cds1Dss = cds.getSequenceAt(0).getDatasetSequence();
+    SequenceI cds2Dss = cds.getSequenceAt(1).getDatasetSequence();
+    assertTrue(dna.getDataset().getSequences().contains(cds1Dss));
+    assertTrue(dna.getDataset().getSequences().contains(cds2Dss));
+
+    /*
+     * verify CDS has a dbref with mapping to peptide
+     */
+    assertNotNull(cds1Dss.getDBRefs());
+    assertEquals(1, cds1Dss.getDBRefs().length);
+    dbref = cds1Dss.getDBRefs()[0];
+    assertEquals("UNIPROT", dbref.getSource());
+    assertEquals("0", dbref.getVersion());
+    assertEquals("pep1", dbref.getAccessionId());
+    assertNotNull(dbref.getMap());
+    assertSame(pep1.getDatasetSequence(), dbref.getMap().getTo());
+    MapList cdsMapping = new MapList(new int[] { 1, 6 },
+            new int[] { 1, 2 }, 3, 1);
+    assertEquals(cdsMapping, dbref.getMap().getMap());
 
     /*
-     * Verify mappings from CDS to peptide, cDNA to CDS, and cDNA to peptide
-     * the mappings are on the shared alignment dataset
+     * verify peptide has added a dbref with reverse mapping to CDS
      */
-    List<AlignedCodonFrame> cdsMappings = cds.getDataset().getCodonFrames();
+    assertNotNull(pep1.getDBRefs());
+    assertEquals(2, pep1.getDBRefs().length);
+    dbref = pep1.getDBRefs()[1];
+    assertEquals("ENSEMBL", dbref.getSource());
+    assertEquals("0", dbref.getVersion());
+    assertEquals("CDS|dna1", dbref.getAccessionId());
+    assertNotNull(dbref.getMap());
+    assertSame(cds1Dss, dbref.getMap().getTo());
+    assertEquals(cdsMapping.getInverse(), dbref.getMap().getMap());
+
     /*
+     * Verify mappings from CDS to peptide, cDNA to CDS, and cDNA to peptide
+     * the mappings are on the shared alignment dataset
      * 6 mappings, 2*(DNA->CDS), 2*(DNA->Pep), 2*(CDS->Pep) 
      */
+    List<AlignedCodonFrame> cdsMappings = cds.getDataset().getCodonFrames();
     assertEquals(6, cdsMappings.size());
 
     /*
@@ -1048,13 +1098,13 @@ public class AlignmentUtilsTests
     SearchResults sr = MappingUtils.buildSearchResults(pep1, 1, mappings);
     assertEquals(1, sr.getResults().size());
     Match m = sr.getResults().get(0);
-    assertSame(cds.getSequenceAt(0).getDatasetSequence(), m.getSequence());
+    assertSame(cds1Dss, m.getSequence());
     assertEquals(1, m.getStart());
     assertEquals(3, m.getEnd());
     // map F to TTT
     sr = MappingUtils.buildSearchResults(pep1, 2, mappings);
     m = sr.getResults().get(0);
-    assertSame(cds.getSequenceAt(0).getDatasetSequence(), m.getSequence());
+    assertSame(cds1Dss, m.getSequence());
     assertEquals(4, m.getStart());
     assertEquals(6, m.getEnd());
 
@@ -1072,19 +1122,19 @@ public class AlignmentUtilsTests
     sr = MappingUtils.buildSearchResults(pep2, 1, mappings);
     assertEquals(1, sr.getResults().size());
     m = sr.getResults().get(0);
-    assertSame(cds.getSequenceAt(1).getDatasetSequence(), m.getSequence());
+    assertSame(cds2Dss, m.getSequence());
     assertEquals(1, m.getStart());
     assertEquals(3, m.getEnd());
     // map F to TTT
     sr = MappingUtils.buildSearchResults(pep2, 2, mappings);
     m = sr.getResults().get(0);
-    assertSame(cds.getSequenceAt(1).getDatasetSequence(), m.getSequence());
+    assertSame(cds2Dss, m.getSequence());
     assertEquals(4, m.getStart());
     assertEquals(6, m.getEnd());
     // map P to CCC
     sr = MappingUtils.buildSearchResults(pep2, 3, mappings);
     m = sr.getResults().get(0);
-    assertSame(cds.getSequenceAt(1).getDatasetSequence(), m.getSequence());
+    assertSame(cds2Dss, m.getSequence());
     assertEquals(7, m.getStart());
     assertEquals(9, m.getEnd());
   }
@@ -1169,7 +1219,7 @@ public class AlignmentUtilsTests
     SequenceI cdsSeq = cds.get(0);
     assertEquals("GGGTTT", cdsSeq.getSequenceAsString());
     // assertEquals("dna1|A12345", cdsSeq.getName());
-    assertEquals("dna1|pep1", cdsSeq.getName());
+    assertEquals("CDS|dna1", cdsSeq.getName());
     // assertEquals(1, cdsSeq.getDBRefs().length);
     // DBRefEntry cdsRef = cdsSeq.getDBRefs()[0];
     // assertEquals("EMBLCDS", cdsRef.getSource());
@@ -1179,7 +1229,7 @@ public class AlignmentUtilsTests
     cdsSeq = cds.get(1);
     assertEquals("aaaccc", cdsSeq.getSequenceAsString());
     // assertEquals("dna1|A12346", cdsSeq.getName());
-    assertEquals("dna1|pep2", cdsSeq.getName());
+    assertEquals("CDS|dna1", cdsSeq.getName());
     // assertEquals(1, cdsSeq.getDBRefs().length);
     // cdsRef = cdsSeq.getDBRefs()[0];
     // assertEquals("EMBLCDS", cdsRef.getSource());
@@ -1189,7 +1239,7 @@ public class AlignmentUtilsTests
     cdsSeq = cds.get(2);
     assertEquals("aaaTTT", cdsSeq.getSequenceAsString());
     // assertEquals("dna1|A12347", cdsSeq.getName());
-    assertEquals("dna1|pep3", cdsSeq.getName());
+    assertEquals("CDS|dna1", cdsSeq.getName());
     // assertEquals(1, cdsSeq.getDBRefs().length);
     // cdsRef = cdsSeq.getDBRefs()[0];
     // assertEquals("EMBLCDS", cdsRef.getSource());
@@ -2255,7 +2305,7 @@ public class AlignmentUtilsTests
      * execute method under test to find CDS for EMBL peptides only
      */
     AlignmentI cds = AlignmentUtils.makeCdsAlignment(new SequenceI[] {
-        dna1, dna2 }, dna.getDataset(), emblPeptides);
+        dna1, dna2 }, dna.getDataset(), emblPeptides.getSequencesArray());
   
     assertEquals(2, cds.getSequences().size());
     assertEquals("GGGTTT", cds.getSequenceAt(0).getSequenceAsString());