JAL-1705 align CDS and peptide products to transcripts
[jalview.git] / src / jalview / ext / ensembl / EnsemblSeqProxy.java
index cbeaae9..869a702 100644 (file)
@@ -11,7 +11,8 @@ import jalview.datamodel.SequenceI;
 import jalview.exceptions.JalviewException;
 import jalview.io.FastaFile;
 import jalview.io.FileParse;
-import jalview.io.gff.SequenceOntology;
+import jalview.io.gff.SequenceOntologyFactory;
+import jalview.io.gff.SequenceOntologyI;
 import jalview.schemes.ResidueProperties;
 import jalview.util.DBRefUtils;
 import jalview.util.MapList;
@@ -36,39 +37,39 @@ import java.util.Map.Entry;
  */
 public abstract class EnsemblSeqProxy extends EnsemblRestClient
 {
+  private static final List<String> CROSS_REFERENCES = Arrays
+          .asList(new String[] { "CCDS", "Uniprot/SWISSPROT" });
+
   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";
-
   protected static final String NAME = "Name";
 
+  /*
+   * enum for 'type' parameter to the /sequence REST service
+   */
   public enum EnsemblSeqType
   {
     /**
-     * type=genomic for the full dna including introns
+     * type=genomic to fetch full dna including introns
      */
     GENOMIC("genomic"),
 
     /**
-     * type=cdna for transcribed dna including UTRs
+     * type=cdna to fetch dna including UTRs
      */
     CDNA("cdna"),
 
     /**
-     * type=cds for coding dna excluding UTRs
+     * type=cds to fetch coding dna excluding UTRs
      */
     CDS("cds"),
 
     /**
-     * type=protein for the peptide product sequence
+     * type=protein to fetch peptide product sequence
      */
     PROTEIN("protein");
 
@@ -124,7 +125,6 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
   @Override
   public AlignmentI getSequenceRecords(String query) throws Exception
   {
-    long now = System.currentTimeMillis();
     // TODO use a String... query vararg instead?
 
     // danger: accession separator used as a regex here, a string elsewhere
@@ -153,17 +153,15 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
                 + " chunks. Unexpected problem (" + r.getLocalizedMessage()
                 + ")";
         System.err.println(msg);
-        if (alignment != null)
-        {
-          break; // return what we got
-        }
-        else
-        {
-          throw new JalviewException(msg, r);
-        }
+        break;
       }
     }
 
+    if (alignment == null)
+    {
+      return null;
+    }
+
     /*
      * fetch and transfer genomic sequence features,
      * fetch protein product and add as cross-reference
@@ -173,9 +171,11 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
       addFeaturesAndProduct(accId, alignment);
     }
 
-    inProgress = false;
-    System.out.println(getClass().getName() + " took "
-            + (System.currentTimeMillis() - now) + "ms to fetch");
+    for (SequenceI seq : alignment.getSequences())
+    {
+      getCrossReferences(seq);
+    }
+
     return alignment;
   }
 
@@ -201,7 +201,7 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
        * get 'dummy' genomic sequence with exon, cds and variation features
        */
       SequenceI genomicSequence = null;
-      EnsemblOverlap gffFetcher = new EnsemblOverlap();
+      EnsemblFeatures gffFetcher = new EnsemblFeatures();
       EnsemblFeatureType[] features = getFeaturesToFetch();
       AlignmentI geneFeatures = gffFetcher.getSequenceRecords(accId,
               features);
@@ -268,7 +268,10 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
       MapList mapList = mapCdsToProtein(querySeq, proteinSeq);
       if (mapList != null)
       {
-        Mapping map = new Mapping(proteinSeq.getDatasetSequence(), mapList);
+        // clunky: ensure Uniprot xref if we have one is on mapped sequence
+        SequenceI ds = proteinSeq.getDatasetSequence();
+        ds.setSourceDBRef(proteinSeq.getSourceDBRef());
+        Mapping map = new Mapping(ds, mapList);
         DBRefEntry dbr = new DBRefEntry(getDbSource(), getDbVersion(),
                 accId, map);
         querySeq.getDatasetSequence().addDBRef(dbr);
@@ -288,6 +291,44 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
   }
 
   /**
+   * Get database xrefs from Ensembl, and attach them to the sequence
+   * 
+   * @param seq
+   */
+  protected void getCrossReferences(SequenceI seq)
+  {
+    while (seq.getDatasetSequence() != null)
+    {
+      seq = seq.getDatasetSequence();
+    }
+
+    EnsemblXref xrefFetcher = new EnsemblXref();
+    List<DBRefEntry> xrefs = xrefFetcher.getCrossReferences(seq.getName(),
+            getCrossReferenceDatabases());
+    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);
+      }
+    }
+  }
+
+  /**
+   * Returns a list of database names to be used when fetching cross-references.
+   * 
+   * @return
+   */
+  protected List<String> getCrossReferenceDatabases()
+  {
+    return CROSS_REFERENCES;
+  }
+
+  /**
    * Returns a mapping from dna to protein by inspecting sequence features of
    * type "CDS" on the dna.
    * 
@@ -302,7 +343,7 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
     int mappedDnaLength = getCdsRanges(dnaSeq, ranges);
 
     int proteinLength = proteinSeq.getLength();
-    List<int[]> proteinRange = new ArrayList<int[]>();
+    int proteinEnd = proteinLength;
     int proteinStart = 1;
 
     /*
@@ -314,26 +355,32 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
       proteinStart = 2;
       proteinLength--;
     }
-    proteinRange.add(new int[] { proteinStart, proteinLength });
+    List<int[]> proteinRange = new ArrayList<int[]>();
 
     /*
      * dna length should map to protein (or protein plus stop codon)
      */
     int codesForResidues = mappedDnaLength / 3;
-    if (codesForResidues == proteinLength
-            || codesForResidues == (proteinLength + 1))
+    if (codesForResidues == (proteinLength + 1))
     {
+      MappingUtils.unmapStopCodon(ranges, mappedDnaLength);
+      codesForResidues--;
+    }
+    if (codesForResidues == proteinLength)
+    {
+      proteinRange.add(new int[] { proteinStart, proteinEnd });
       return new MapList(ranges, proteinRange, 3, 1);
     }
     return null;
   }
 
   /**
-   * Adds CDS ranges to the ranges list, and returns the total length mapped.
+   * Adds CDS ranges to the ranges list, and returns the total length mapped
+   * from.
    * 
-   * No need to worry about reverse strand dna here since the retrieved sequence
-   * is as transcribed (reverse complement for reverse strand), i.e in the same
-   * sense as the peptide.
+   * No need to worry about reverse strand dna, here since the retrieved
+   * sequence is as transcribed (reverse complement for reverse strand), i.e in
+   * the same sense as the peptide.
    * 
    * @param dnaSeq
    * @param ranges
@@ -346,13 +393,14 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
     {
       return 0;
     }
+    SequenceOntologyI so = SequenceOntologyFactory.getInstance();
     int mappedDnaLength = 0;
     for (SequenceFeature sf : sfs)
     {
       /*
        * process a CDS feature (or a sub-type of CDS)
        */
-      if (SequenceOntology.getInstance().isA(sf.getType(), SequenceOntology.CDS))
+      if (so.isA(sf.getType(), SequenceOntologyI.CDS))
       {
         int phase = 0;
         try {
@@ -367,7 +415,7 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
          */
         int begin = sf.getBegin();
         int end = sf.getEnd();
-        if (ranges.isEmpty() && phase > 0)
+        if (ranges.isEmpty())
         {
           begin += phase;
           if (begin > end)
@@ -540,7 +588,7 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
    *          the start position of the sequence we are mapping to
    * @return
    */
-  protected MapList getGenomicRanges(SequenceI sourceSequence,
+  protected MapList getGenomicRangesFromFeatures(SequenceI sourceSequence,
           String accId, int start)
   {
     SequenceFeature[] sfs = sourceSequence.getSequenceFeatures();
@@ -566,11 +614,12 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
        */
       if (identifiesSequence(sf, accId))
       {
-          int strand = sf.getStrand();
-  
-          if (directionSet && strand != direction)
-          {
-            // abort - mix of forward and backward
+        int strand = sf.getStrand();
+        strand = strand == 0 ? 1 : strand; // treat unknown as forward
+
+        if (directionSet && strand != direction)
+        {
+          // abort - mix of forward and backward
           System.err.println("Error: forward and backward strand for "
                   + accId);
             return null;
@@ -615,8 +664,8 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
      */
     Collections.sort(regions, new RangeSorter(direction == 1));
   
-    List<int[]> to = new ArrayList<int[]>();
-    to.add(new int[] { start, start + mappedLength - 1 });
+    List<int[]> to = Arrays.asList(new int[] { start,
+        start + mappedLength - 1 });
   
     return new MapList(regions, to, 1, 1);
   }
@@ -671,17 +720,18 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
       /*
        * for sequence_variant, make an additional feature with consequence
        */
-      if (SequenceOntology.getInstance().isSequenceVariant(sf.getType()))
-      {
-        String consequence = (String) sf.getValue(CONSEQUENCE_TYPE);
-        if (consequence != null)
-        {
-          SequenceFeature sf2 = new SequenceFeature("consequence",
-                  consequence, copy.getBegin(), copy.getEnd(), 0f,
-                  null);
-          targetSequence.addSequenceFeature(sf2);
-        }
-      }
+      // if (SequenceOntologyFactory.getInstance().isA(sf.getType(),
+      // SequenceOntologyI.SEQUENCE_VARIANT))
+      // {
+      // String consequence = (String) sf.getValue(CONSEQUENCE_TYPE);
+      // if (consequence != null)
+      // {
+      // SequenceFeature sf2 = new SequenceFeature("consequence",
+      // consequence, copy.getBegin(), copy.getEnd(), 0f,
+      // null);
+      // targetSequence.addSequenceFeature(sf2);
+      // }
+      // }
     }
   }
 
@@ -701,15 +751,22 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
       return false;
     }
 
+    // long start = System.currentTimeMillis();
     SequenceFeature[] sfs = sourceSequence.getSequenceFeatures();
-    MapList mapping = getGenomicRanges(sourceSequence, accessionId,
+    MapList mapping = getGenomicRangesFromFeatures(sourceSequence, accessionId,
             targetSequence.getStart());
     if (mapping == null)
     {
       return false;
     }
 
-    return transferFeatures(sfs, targetSequence, mapping, accessionId);
+    boolean result = transferFeatures(sfs, targetSequence, mapping,
+            accessionId);
+    // System.out.println("transferFeatures (" + (sfs.length) + " --> "
+    // + targetSequence.getSequenceFeatures().length + ") to "
+    // + targetSequence.getName()
+    // + " took " + (System.currentTimeMillis() - start) + "ms");
+    return result;
   }
 
   /**
@@ -811,7 +868,7 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
     
     SequenceFeature[] sfs = sequence.getSequenceFeatures();
     if (sfs != null) {
-      SequenceOntology so = SequenceOntology.getInstance();
+      SequenceOntologyI so = SequenceOntologyFactory.getInstance();
       for (SequenceFeature sf :sfs) {
         if (so.isA(sf.getType(), type))
         {
@@ -849,7 +906,7 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
     }
   
     AlignmentUtils.transferFeatures(dnaSeq, peptide, dnaToProtein,
-            SequenceOntology.EXON);
+            SequenceOntologyI.EXON);
 
     LinkedHashMap<Integer, String[][]> variants = buildDnaVariantsMap(
             dnaSeq, dnaToProtein);
@@ -870,12 +927,28 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
         String desc = StringUtils.listToDelimitedString(peptideVariants,
                 ", ");
         SequenceFeature sf = new SequenceFeature(
-                SequenceOntology.SEQUENCE_VARIANT, desc, peptidePos,
+                SequenceOntologyI.SEQUENCE_VARIANT, desc, peptidePos,
                 peptidePos, 0f, null);
         peptide.addSequenceFeature(sf);
         count++;
       }
     }
+
+    /*
+     * ugly sort to get sequence features in start position order
+     * - would be better to store in Sequence as a TreeSet instead?
+     */
+    Arrays.sort(peptide.getSequenceFeatures(),
+            new Comparator<SequenceFeature>()
+            {
+              @Override
+              public int compare(SequenceFeature o1, SequenceFeature o2)
+              {
+                int c = Integer.compare(o1.getBegin(), o2.getBegin());
+                return c == 0 ? Integer.compare(o1.getEnd(), o2.getEnd())
+                        : c;
+              }
+            });
     return count;
   }
 
@@ -895,7 +968,7 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
      * LinkedHashMap ensures we add the peptide features in sequence order
      */
     LinkedHashMap<Integer, String[][]> variants = new LinkedHashMap<Integer, String[][]>();
-    SequenceOntology so = SequenceOntology.getInstance();
+    SequenceOntologyI so = SequenceOntologyFactory.getInstance();
   
     SequenceFeature[] dnaFeatures = dnaSeq.getSequenceFeatures();
     if (dnaFeatures == null)
@@ -918,7 +991,7 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
         // not handling multi-locus variant features
         continue;
       }
-      if (so.isSequenceVariant(sf.getType()))
+      if (so.isA(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT))
       {
         int[] mapsTo = dnaToProtein.locateInTo(dnaCol, dnaCol);
         if (mapsTo == null)
@@ -1056,7 +1129,8 @@ public abstract class EnsemblSeqProxy extends EnsemblRestClient
    */
   public static boolean isTranscript(String featureType)
   {
-    return NMD_VARIANT.equals(featureType)
-            || SequenceOntology.getInstance().isA(featureType, SequenceOntology.TRANSCRIPT);
+    return SequenceOntologyI.NMD_TRANSCRIPT_VARIANT.equals(featureType)
+            || SequenceOntologyFactory.getInstance().isA(featureType,
+                    SequenceOntologyI.TRANSCRIPT);
   }
 }