JAL-1705 derive transcript sequences from gene sequence/GFF; handle CDS
[jalview.git] / src / jalview / ext / ensembl / EnsemblSeqProxy.java
index 7153a9e..d7811b9 100644 (file)
@@ -30,10 +30,18 @@ import java.util.List;
  */
 public abstract class EnsemblSeqProxy extends EnsemblRestClient
 {
+  protected static final String CONSEQUENCE_TYPE = "consequence_type";
+
   protected static final String PARENT = "Parent";
 
   protected static final String ID = "ID";
 
+  /*
+   * this needs special handling, as it isA sequence_variant in the
+   * Sequence Ontology, but behaves in Ensembl as if it isA transcript
+   */
+  protected static final String NMD_VARIANT = "NMD_transcript_variant";
+
   public enum EnsemblSeqType
   {
     /**
@@ -94,11 +102,6 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
 
   }
 
-  /*
-   * genomic sequence, with features retrieved from the REST overlap service
-   */
-  private SequenceI genomicSequence;
-
   /**
    * Constructor
    */
@@ -108,17 +111,20 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
 
   /**
    * Makes the sequence queries to Ensembl's REST service and returns an
-   * alignment consisting of the returned sequences
+   * alignment consisting of the returned sequences. This overloaded method
+   * allows the genomic sequence (with features) to be passed in if it has
+   * already been retrieved, to avoid repeat calls to fetch it.
    */
-  @Override
-  public AlignmentI getSequenceRecords(String query) throws Exception
+  public AlignmentI getSequenceRecords(String query,
+          SequenceI genomicSequence) throws Exception
   {
     long now = System.currentTimeMillis();
     // TODO use a String... query vararg instead?
 
     // danger: accession separator used as a regex here, a string elsewhere
     // in this case it is ok (it is just a space), but (e.g.) '\' would not be
-    List<String> allIds = Arrays.asList(query.split(getAccessionSeparator()));
+    List<String> allIds = Arrays.asList(query
+            .split(getAccessionSeparator()));
     AlignmentI alignment = null;
     inProgress = true;
 
@@ -157,7 +163,7 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
      */
     for (String accId : allIds)
     {
-      addFeaturesAndProduct(accId, alignment);
+      addFeaturesAndProduct(accId, alignment, genomicSequence);
     }
 
     inProgress = false;
@@ -175,7 +181,8 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
    * @param accId
    * @param alignment
    */
-  protected void addFeaturesAndProduct(String accId, AlignmentI alignment)
+  protected void addFeaturesAndProduct(String accId, AlignmentI alignment,
+          SequenceI genomicSequence)
   {
     try
     {
@@ -190,14 +197,14 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
                 features);
         if (geneFeatures.getHeight() > 0)
         {
-          /*
-           * transfer features to the query sequence
-           */
           genomicSequence = geneFeatures.getSequenceAt(0);
         }
       }
       if (genomicSequence != null)
       {
+        /*
+         * transfer features to the query sequence
+         */
         SequenceI querySeq = alignment.findName(accId);
         if (transferFeatures(accId, genomicSequence, querySeq))
         {
@@ -291,6 +298,7 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
      * dna here since the retrieved sequence is as transcribed (reverse
      * complement for reverse strand), i.e in the same sense as the peptide. 
      */
+    boolean fivePrimeIncomplete = false;
     for (SequenceFeature sf : sfs)
     {
       /*
@@ -298,19 +306,48 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
        */
       if (so.isA(sf.getType(), SequenceOntology.CDS))
       {
-        ranges.add(new int[] { sf.getBegin(), sf.getEnd() });
-        mappedDnaLength += Math.abs(sf.getEnd() - sf.getBegin()) + 1;
+        int phase = 0;
+        try {
+          phase = Integer.parseInt(sf.getPhase());
+        } catch (NumberFormatException e)
+        {
+          // ignore
+        }
+        /*
+         * phase > 0 on first codon means 5' incomplete - skip to the start
+         * of the next codon; example ENST00000496384
+         */
+        int begin = sf.getBegin();
+        int end = sf.getEnd();
+        if (ranges.isEmpty() && phase > 0)
+        {
+          fivePrimeIncomplete = true;
+          begin += phase;
+          if (begin > end)
+          {
+            continue; // shouldn't happen?
+          }
+        }
+        ranges.add(new int[] { begin, end });
+        mappedDnaLength += Math.abs(end - begin) + 1;
       }
     }
     int proteinLength = proteinSeq.getLength();
     List<int[]> proteinRange = new ArrayList<int[]>();
-    proteinRange.add(new int[] { 1, proteinLength });
+    int proteinStart = 1;
+    if (fivePrimeIncomplete && proteinSeq.getCharAt(0) == 'X')
+    {
+      proteinStart = 2;
+      proteinLength--;
+    }
+    proteinRange.add(new int[] { proteinStart, proteinLength });
 
     /*
-     * dna length should map to protein (or protein minus stop codon)
+     * dna length should map to protein (or protein plus stop codon)
      */
-    if (mappedDnaLength == 3 * proteinLength
-            || mappedDnaLength == 3 * (proteinLength + 1))
+    int codesForResidues = mappedDnaLength / 3;
+    if (codesForResidues == proteinLength
+            || codesForResidues == (proteinLength + 1))
     {
       return new MapList(ranges, proteinRange, 3, 1);
     }
@@ -470,11 +507,10 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
     }
 
     /*
-     * generously size for initial number of cds regions
+     * generously initial size for number of cds regions
      * (worst case titin Q8WZ42 has c. 313 exons)
      */
     List<int[]> regions = new ArrayList<int[]>(100);
-    int sourceLength = sourceSequence.getLength();
     int mappedLength = 0;
     int direction = 1; // forward
     boolean directionSet = false;
@@ -512,11 +548,11 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
         }
         mappedLength += Math.abs(sf.getEnd() - sf.getBegin() + 1);
 
-        if (mappedLength >= sourceLength)
+        if (!isSpliceable())
         {
           /*
-           * break for the case of matching gene features v gene sequence
-           * - only need to locate the 'gene' feature for accId
+           * 'gene' sequence is contiguous so we can stop as soon as its
+           * identifying feature has been found
            */
           break;
         }
@@ -543,6 +579,15 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
   }
 
   /**
+   * Answers true if the sequence being retrieved may occupy discontiguous
+   * regions on the genomic sequence.
+   */
+  protected boolean isSpliceable()
+  {
+    return true;
+  }
+
+  /**
    * Returns true if the sequence feature marks positions of the genomic
    * sequence feature which are within the sequence being retrieved. For
    * example, an 'exon' feature whose parent is the target transcript marks the
@@ -556,21 +601,22 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
           String accId);
 
   /**
-   * Transfers the sequence feature to the target sequence, adjusting its start
-   * and end range based on the 'overlap' ranges. Features which do not overlap
-   * the target sequence are ignored, as are features with a parent other than
-   * the target sequence id.
+   * Transfers the sequence feature to the target sequence, locating its start
+   * and end range based on the mapping. Features which do not overlap the
+   * target sequence are ignored.
    * 
    * @param sf
    * @param targetSequence
-   * @param overlap
+   * @param mapping
+   *          mapping from the sequence feature's coordinates to the target
+   *          sequence
    */
   protected void transferFeature(SequenceFeature sf,
-          SequenceI targetSequence, MapList overlap)
+          SequenceI targetSequence, MapList mapping)
   {
     int start = sf.getBegin();
     int end = sf.getEnd();
-    int[] mappedRange = overlap.locateInTo(start, end);
+    int[] mappedRange = mapping.locateInTo(start, end);
   
     if (mappedRange != null)
     {
@@ -584,7 +630,7 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
        */
       if (SequenceOntology.getInstance().isSequenceVariant(sf.getType()))
       {
-        String consequence = (String) sf.getValue("consequence_type");
+        String consequence = (String) sf.getValue(CONSEQUENCE_TYPE);
         if (consequence != null)
         {
           SequenceFeature sf2 = new SequenceFeature("consequence",
@@ -594,7 +640,6 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
         }
       }
     }
-  
   }
 
   /**
@@ -614,20 +659,37 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
     }
 
     SequenceFeature[] sfs = sourceSequence.getSequenceFeatures();
-    MapList overlap = getGenomicRanges(sourceSequence, accessionId,
+    MapList mapping = getGenomicRanges(sourceSequence, accessionId,
             targetSequence.getStart());
-    if (overlap == null)
+    if (mapping == null)
     {
       return false;
     }
 
-    final boolean forwardStrand = overlap.isFromForwardStrand();
+    return transferFeatures(sfs, targetSequence, mapping, accessionId);
+  }
+
+  /**
+   * Transfer features to the target sequence. The start/end positions are
+   * converted using the mapping. Features which do not overlap are ignored.
+   * Features whose parent is not the specified identifier are also ignored.
+   * 
+   * @param features
+   * @param targetSequence
+   * @param mapping
+   * @param parentId
+   * @return
+   */
+  protected boolean transferFeatures(SequenceFeature[] features,
+          SequenceI targetSequence, MapList mapping, String parentId)
+  {
+    final boolean forwardStrand = mapping.isFromForwardStrand();
 
     /*
      * sort features by start position (descending if reverse strand) 
      * before transferring (in forwards order) to the target sequence
      */
-    Arrays.sort(sfs, new Comparator<SequenceFeature>()
+    Arrays.sort(features, new Comparator<SequenceFeature>()
     {
       @Override
       public int compare(SequenceFeature o1, SequenceFeature o2)
@@ -638,11 +700,11 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
     });
 
     boolean transferred = false;
-    for (SequenceFeature sf : sfs)
+    for (SequenceFeature sf : features)
     {
-      if (retainFeature(sf, accessionId))
+      if (retainFeature(sf, parentId))
       {
-        transferFeature(sf, targetSequence, overlap);
+        transferFeature(sf, targetSequence, mapping);
         transferred = true;
       }
     }
@@ -650,15 +712,31 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
   }
 
   /**
-   * Answers true if the feature is one to attach to the retrieved sequence
+   * Answers true if the feature type is one we want to keep for the sequence.
+   * Some features are only retrieved in order to identify the sequence range,
+   * and may then be discarded as redundant information (e.g. "CDS" feature for
+   * a CDS sequence).
+   */
+  @SuppressWarnings("unused")
+  protected boolean retainFeature(SequenceFeature sf, String accessionId)
+  {
+    return true; // override as required
+  }
+
+  /**
+   * Answers true if the feature has a Parent which refers to the given
+   * accession id, or if the feature has no parent. Answers false if the
+   * feature's Parent is for a different accession id.
    * 
-   * @param type
+   * @param sf
+   * @param identifier
    * @return
    */
-  protected boolean retainFeature(SequenceFeature sf, String accessionId)
+  protected boolean featureMayBelong(SequenceFeature sf, String identifier)
   {
     String parent = (String) sf.getValue(PARENT);
-    if (parent != null && !parent.contains(accessionId))
+    // using contains to allow for prefix "gene:", "transcript:" etc
+    if (parent != null && !parent.contains(identifier))
     {
       // this genomic feature belongs to a different transcript
       return false;
@@ -673,10 +751,57 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
             + " sequence with variant features";
   }
 
-  public AlignmentI getSequenceRecords(String transcriptId,
-          SequenceI geneSeq) throws Exception
+  @Override
+  public AlignmentI getSequenceRecords(String identifier) throws Exception
+  {
+    return getSequenceRecords(identifier, null);
+  }
+
+  /**
+   * Returns a (possibly empty) list of features on the sequence which have the
+   * specified sequence ontology type (or a sub-type of it), and the given
+   * identifier as parent
+   * 
+   * @param sequence
+   * @param type
+   * @param parentId
+   * @return
+   */
+  protected List<SequenceFeature> findFeatures(SequenceI sequence,
+          String type, String parentId)
+  {
+    List<SequenceFeature> result = new ArrayList<SequenceFeature>();
+    
+    SequenceFeature[] sfs = sequence.getSequenceFeatures();
+    if (sfs != null) {
+      SequenceOntology so = SequenceOntology.getInstance();
+      for (SequenceFeature sf :sfs) {
+        if (so.isA(sf.getType(), type))
+        {
+          String parent = (String) sf.getValue(PARENT);
+          if (parent.equals(parentId))
+          {
+            result.add(sf);
+          }
+        }
+      }
+    }
+    return result;
+  }
+
+  /**
+   * Answers true if the feature type is either 'NMD_transcript_variant' or
+   * 'transcript' or one of its sub-types in the Sequence Ontology. This is
+   * needed because NMD_transcript_variant behaves like 'transcript' in Ensembl
+   * although strictly speaking it is not (it is a sub-type of
+   * sequence_variant).
+   * 
+   * @param featureType
+   * @return
+   */
+  public static boolean isTranscript(String featureType)
   {
-    this.genomicSequence = geneSeq;
-    return getSequenceRecords(transcriptId);
+    return NMD_VARIANT.equals(featureType)
+            || SequenceOntology.getInstance().isA(featureType, SequenceOntology.TRANSCRIPT);
   }
 }