JAL-3187 CDS has same dataset sequence as transcript
[jalview.git] / src / jalview / analysis / AlignmentUtils.java
index 5e11446..77cc4d6 100644 (file)
@@ -74,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
@@ -1743,8 +1746,10 @@ 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[]
+                  { cdsSeq.getStart(),
+                      cdsSeq.getLength() + cdsSeq.getStart() - 1 });
           MapList cdsToProteinMap = new MapList(cdsRange,
                   mapList.getToRanges(), mapList.getFromRatio(),
                   mapList.getToRatio());
@@ -1867,7 +1872,7 @@ public class AlignmentUtils
       return;
     }
 
-    MapList newMap = targetToFrom.traverse(fromLoci.getMap());
+    MapList newMap = targetToFrom.traverse(fromLoci.getMapping());
 
     if (newMap != null)
     {
@@ -1888,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,
@@ -1913,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;
+      }
     }
 
     /*
@@ -1938,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())
@@ -1971,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());
@@ -2056,9 +2093,11 @@ public class AlignmentUtils
   protected static List<DBRefEntry> propagateDBRefsToCDS(SequenceI cdsSeq,
           SequenceI contig, SequenceI proteinProduct, Mapping 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())
@@ -2144,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
@@ -2235,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();
@@ -2412,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;
@@ -2562,15 +2607,15 @@ public class AlignmentUtils
             peptidePos, var.getSource());
 
     StringBuilder attributes = new StringBuilder(32);
-    String id = (String) var.variant.getValue(ID);
+    String id = (String) var.variant.getValue(VARIANT_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);
+      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
@@ -2792,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);
@@ -2863,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)
   {
@@ -2890,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);
     }
 
     /*
@@ -2906,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)