JAL-2789 avoid using transcript (with no UTR) as CDS sequence
[jalview.git] / src / jalview / analysis / AlignmentUtils.java
index 2b9b9f9..343ebc7 100644 (file)
@@ -117,7 +117,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())
     {
@@ -247,7 +247,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();
@@ -256,7 +256,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);
@@ -283,8 +283,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
@@ -870,7 +870,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);
@@ -1081,7 +1081,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())
@@ -1127,9 +1127,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())
@@ -1308,7 +1308,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);
@@ -1445,7 +1445,7 @@ public class AlignmentUtils
       {
         continue;
       }
-      final List<AlignmentAnnotation> result = new ArrayList<AlignmentAnnotation>();
+      final List<AlignmentAnnotation> result = new ArrayList<>();
       for (AlignmentAnnotation dsann : datasetAnnotations)
       {
         /*
@@ -1627,13 +1627,13 @@ 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
@@ -1833,7 +1833,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 +1858,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 +1891,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())
@@ -2002,8 +2012,8 @@ public class AlignmentUtils
   {
 
     // gather direct refs from contig congrent with mapping
-    List<DBRefEntry> direct = new ArrayList<DBRefEntry>();
-    HashSet<String> directSources = new HashSet<String>();
+    List<DBRefEntry> direct = new ArrayList<>();
+    HashSet<String> directSources = new HashSet<>();
     if (contig.getDBRefs() != null)
     {
       for (DBRefEntry dbr : contig.getDBRefs())
@@ -2023,7 +2033,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)
@@ -2164,7 +2174,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
@@ -2176,6 +2189,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();
@@ -2190,7 +2213,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)
@@ -2199,8 +2222,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 });
@@ -2211,7 +2238,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.
@@ -2221,7 +2248,7 @@ public class AlignmentUtils
    */
   public 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);
@@ -2230,7 +2257,6 @@ public class AlignmentUtils
       return result;
     }
     SequenceFeatures.sortFeatures(sfs, true);
-    int startPhase = 0;
 
     for (SequenceFeature sf : sfs)
     {
@@ -2248,7 +2274,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)
@@ -2263,16 +2289,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
@@ -2518,7 +2534,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);
@@ -2553,9 +2569,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);
       }
 
@@ -2694,7 +2710,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();
@@ -2769,7 +2785,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();
@@ -2832,7 +2848,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
@@ -2937,7 +2953,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));