JAL-2738 use GeneLocus extends DBRefEntry to hold chromosomal mappings
[jalview.git] / src / jalview / ext / ensembl / EnsemblGene.java
index 3b32797..9c7ce63 100644 (file)
@@ -1,17 +1,38 @@
+/*
+ * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
+ * Copyright (C) $$Year-Rel$$ The Jalview Authors
+ * 
+ * This file is part of Jalview.
+ * 
+ * Jalview is free software: you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License 
+ * as published by the Free Software Foundation, either version 3
+ * of the License, or (at your option) any later version.
+ *  
+ * Jalview is distributed in the hope that it will be useful, but 
+ * WITHOUT ANY WARRANTY; without even the implied warranty 
+ * of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
+ * PURPOSE.  See the GNU General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
+ * The Jalview Authors are detailed in the 'AUTHORS' file.
+ */
 package jalview.ext.ensembl;
 
 import jalview.api.FeatureColourI;
-import jalview.api.FeatureSettingsI;
+import jalview.api.FeatureSettingsModelI;
 import jalview.datamodel.AlignmentI;
+import jalview.datamodel.GeneLociI;
 import jalview.datamodel.Sequence;
 import jalview.datamodel.SequenceFeature;
 import jalview.datamodel.SequenceI;
+import jalview.datamodel.features.SequenceFeatures;
 import jalview.io.gff.SequenceOntologyFactory;
 import jalview.io.gff.SequenceOntologyI;
-import jalview.schemes.FeatureColourAdapter;
+import jalview.schemes.FeatureColour;
 import jalview.schemes.FeatureSettingsAdapter;
 import jalview.util.MapList;
-import jalview.util.StringUtils;
 
 import java.awt.Color;
 import java.io.UnsupportedEncodingException;
@@ -29,11 +50,6 @@ import com.stevesoft.pat.Regex;
  */
 public class EnsemblGene extends EnsemblSeqProxy
 {
-  private static final List<String> CROSS_REFERENCES = Arrays
-          .asList(new String[] { "CCDS" });
-
-  private static final String GENE_PREFIX = "gene:";
-
   /*
    * accepts anything as we will attempt lookup of gene or 
    * transcript id or gene name
@@ -45,6 +61,8 @@ public class EnsemblGene extends EnsemblSeqProxy
       EnsemblFeatureType.exon, EnsemblFeatureType.cds,
       EnsemblFeatureType.variation };
 
+  private static final String CHROMOSOME = "chromosome";
+
   /**
    * Default constructor (to use rest.ensembl.org)
    */
@@ -81,6 +99,12 @@ public class EnsemblGene extends EnsemblSeqProxy
     return EnsemblSeqType.GENOMIC;
   }
 
+  @Override
+  protected String getObjectType()
+  {
+    return OBJECT_TYPE_GENE;
+  }
+
   /**
    * Returns an alignment containing the gene(s) for the given gene or
    * transcript identifier, or external identifier (e.g. Uniprot id). If given a
@@ -93,7 +117,8 @@ public class EnsemblGene extends EnsemblSeqProxy
    * <li>resolves an external identifier by looking up xref-ed gene ids</li>
    * <li>fetches the gene sequence</li>
    * <li>fetches features on the sequence</li>
-   * <li>identifies "transcript" features whose Parent is the requested gene</li>
+   * <li>identifies "transcript" features whose Parent is the requested
+   * gene</li>
    * <li>fetches the transcript sequence for each transcript</li>
    * <li>makes a mapping from the gene to each transcript</li>
    * <li>copies features from gene to transcript sequences</li>
@@ -104,82 +129,154 @@ public class EnsemblGene extends EnsemblSeqProxy
    * </ul>
    * 
    * @param query
-   *          one or more identifiers separated by a space
-   * @return an alignment containing one or more genes, and possibly
-   *         transcripts, or null
+   *          a single gene or transcript identifier or gene name
+   * @return an alignment containing a gene, and possibly transcripts, or null
    */
   @Override
   public AlignmentI getSequenceRecords(String query) throws Exception
   {
-    // todo: tidy up handling of one or multiple accession ids
-    String[] queries = query.split(getAccessionSeparator());
-
     /*
-     * if given a transcript id, look up its gene parent
+     * convert to a non-duplicated list of gene identifiers
      */
-    if (isTranscriptIdentifier(query))
+    List<String> geneIds = getGeneIds(query);
+
+    AlignmentI al = null;
+    for (String geneId : geneIds)
     {
-      // we are assuming all transcripts have the same gene parent here
-      query = new EnsemblLookup(getDomain()).getParent(queries[0]);
-      if (query == null)
+      /*
+       * fetch the gene sequence(s) with features and xrefs
+       */
+      AlignmentI geneAlignment = super.getSequenceRecords(geneId);
+      if (geneAlignment == null)
       {
-        return null;
+        continue;
       }
-    }
+      
+      if (geneAlignment.getHeight() == 1)
+      {
+        // ensure id has 'correct' case for the Ensembl identifier
+        geneId = geneAlignment.getSequenceAt(0).getName();
 
-    /*
-     * if given a gene or other external name, lookup and fetch 
-     * the corresponding gene for all model organisms 
-     */
-    if (!isGeneIdentifier(query))
-    {
-      List<String> geneIds = new EnsemblSymbol(getDomain()).getIds(query);
-      if (geneIds.isEmpty())
+        findGeneLoci(geneAlignment.getSequenceAt(0), geneId);
+
+        getTranscripts(geneAlignment, geneId);
+      }
+      if (al == null)
       {
-        return null;
+        al = geneAlignment;
+      }
+      else
+      {
+        al.append(geneAlignment);
       }
-      String theIds = StringUtils.listToDelimitedString(geneIds,
-              getAccessionSeparator());
-      return getSequenceRecords(theIds);
     }
+    return al;
+  }
 
-    /*
-     * fetch the gene sequence(s) with features and xrefs
-     */
-    AlignmentI al = super.getSequenceRecords(query);
-
-    /*
-     * if we retrieved a single gene, get its transcripts as well
-     */
-    if (al.getHeight() == 1)
+  /**
+   * Calls the /lookup/id REST service, parses the response for gene
+   * coordinates, and if successful, adds these to the sequence. If this fails,
+   * fall back on trying to parse the sequence description in case it is in
+   * Ensembl-gene format e.g. chromosome:GRCh38:17:45051610:45109016:1.
+   * 
+   * @param seq
+   * @param geneId
+   */
+  void findGeneLoci(SequenceI seq, String geneId)
+  {
+    GeneLociI geneLoci = new EnsemblLookup(getDomain()).getGeneLoci(geneId);
+    if (geneLoci != null)
     {
-      getTranscripts(al, query);
+      seq.setGeneLoci(geneLoci.getSpeciesId(), geneLoci.getAssemblyId(),
+              geneLoci.getChromosomeId(), geneLoci.getMapping());
     }
+    else
+    {
+      parseChromosomeLocations(seq);
+    }
+  }
 
-    return al;
+  /**
+   * Parses and saves fields of an Ensembl-style description e.g.
+   * chromosome:GRCh38:17:45051610:45109016:1
+   * 
+   * @param seq
+   */
+  boolean parseChromosomeLocations(SequenceI seq)
+  {
+    String description = seq.getDescription();
+    if (description == null)
+    {
+      return false;
+    }
+    String[] tokens = description.split(":");
+    if (tokens.length == 6 && tokens[0].startsWith(CHROMOSOME))
+    {
+      String ref = tokens[1];
+      String chrom = tokens[2];
+      try
+      {
+        int chStart = Integer.parseInt(tokens[3]);
+        int chEnd = Integer.parseInt(tokens[4]);
+        boolean forwardStrand = "1".equals(tokens[5]);
+        String species = ""; // not known here
+        int[] from = new int[] { seq.getStart(), seq.getEnd() };
+        int[] to = new int[] { forwardStrand ? chStart : chEnd,
+            forwardStrand ? chEnd : chStart };
+        MapList map = new MapList(from, to, 1, 1);
+        seq.setGeneLoci(species, ref, chrom, map);
+        return true;
+      } catch (NumberFormatException e)
+      {
+        System.err.println("Bad integers in description " + description);
+      }
+    }
+    return false;
   }
 
   /**
-   * Attempts to get Ensembl stable identifiers for model organisms for a gene
-   * name by calling the xrefs symbol REST service to resolve the gene name.
+   * Converts a query, which may contain one or more gene, transcript, or
+   * external (to Ensembl) identifiers, into a non-redundant list of gene
+   * identifiers.
    * 
-   * @param query
+   * @param accessions
    * @return
    */
-  protected String getGeneIdentifiersForName(String query)
+  List<String> getGeneIds(String accessions)
   {
-    List<String> ids = new EnsemblSymbol(getDomain()).getIds(query);
-    if (ids != null)
+    List<String> geneIds = new ArrayList<>();
+
+    for (String acc : accessions.split(getAccessionSeparator()))
     {
-      for (String id : ids)
+      /*
+       * First try lookup as an Ensembl (gene or transcript) identifier
+       */
+      String geneId = new EnsemblLookup(getDomain()).getGeneId(acc);
+      if (geneId != null)
       {
-        if (isGeneIdentifier(id))
+        if (!geneIds.contains(geneId))
         {
-          return id;
+          geneIds.add(geneId);
+        }
+      }
+      else
+      {
+        /*
+         * if given a gene or other external name, lookup and fetch 
+         * the corresponding gene for all model organisms 
+         */
+        List<String> ids = new EnsemblSymbol(getDomain(), getDbSource(),
+                getDbVersion()).getGeneIds(acc);
+        for (String id : ids)
+        {
+          if (!geneIds.contains(id))
+          {
+            geneIds.add(id);
+          }
         }
       }
     }
-    return null;
+    return geneIds;
   }
 
   /**
@@ -214,23 +311,20 @@ public class EnsemblGene extends EnsemblSeqProxy
    */
   protected void clearGeneFeatures(SequenceI gene)
   {
-    SequenceFeature[] sfs = gene.getSequenceFeatures();
-    if (sfs != null)
+    /*
+     * Note we include NMD_transcript_variant here because it behaves like 
+     * 'transcript' in Ensembl, although strictly speaking it is not 
+     * (it is a sub-type of sequence_variant)    
+     */
+    String[] soTerms = new String[] {
+        SequenceOntologyI.NMD_TRANSCRIPT_VARIANT,
+        SequenceOntologyI.TRANSCRIPT, SequenceOntologyI.EXON,
+        SequenceOntologyI.CDS };
+    List<SequenceFeature> sfs = gene.getFeatures().getFeaturesByOntology(
+            soTerms);
+    for (SequenceFeature sf : sfs)
     {
-      SequenceOntologyI so = SequenceOntologyFactory.getInstance();
-      List<SequenceFeature> filtered = new ArrayList<SequenceFeature>();
-      for (SequenceFeature sf : sfs)
-      {
-        String type = sf.getType();
-        if (!isTranscript(type) && !so.isA(type, SequenceOntologyI.EXON)
-                && !so.isA(type, SequenceOntologyI.CDS))
-        {
-          filtered.add(sf);
-        }
-      }
-      gene.setSequenceFeatures(filtered
-              .toArray(new SequenceFeature[filtered
-              .size()]));
+      gene.deleteFeature(sf);
     }
   }
 
@@ -247,8 +341,8 @@ public class EnsemblGene extends EnsemblSeqProxy
    *          the parent gene sequence, with features
    * @return
    */
-  SequenceI makeTranscript(SequenceFeature transcriptFeature,
-          AlignmentI al, SequenceI gene)
+  SequenceI makeTranscript(SequenceFeature transcriptFeature, AlignmentI al,
+          SequenceI gene)
   {
     String accId = getTranscriptId(transcriptFeature);
     if (accId == null)
@@ -273,18 +367,19 @@ public class EnsemblGene extends EnsemblSeqProxy
      * look for exon features of the transcript, failing that for CDS
      * (for example ENSG00000124610 has 1 CDS but no exon features)
      */
-    String parentId = "transcript:" + accId;
+    String parentId = accId;
     List<SequenceFeature> splices = findFeatures(gene,
             SequenceOntologyI.EXON, parentId);
     if (splices.isEmpty())
     {
       splices = findFeatures(gene, SequenceOntologyI.CDS, parentId);
     }
+    SequenceFeatures.sortFeatures(splices, true);
 
     int transcriptLength = 0;
     final char[] geneChars = gene.getSequence();
     int offset = gene.getStart(); // to convert to 0-based positions
-    List<int[]> mappedFrom = new ArrayList<int[]>();
+    List<int[]> mappedFrom = new ArrayList<>();
 
     for (SequenceFeature sf : splices)
     {
@@ -296,13 +391,14 @@ public class EnsemblGene extends EnsemblSeqProxy
       mappedFrom.add(new int[] { sf.getBegin(), sf.getEnd() });
     }
 
-    Sequence transcript = new Sequence(accId, seqChars, 1, transcriptLength);
+    Sequence transcript = new Sequence(accId, seqChars, 1,
+            transcriptLength);
 
     /*
      * Ensembl has gene name as transcript Name
      * EnsemblGenomes doesn't, but has a url-encoded description field
      */
-    String description = (String) transcriptFeature.getValue(NAME);
+    String description = transcriptFeature.getDescription();
     if (description == null)
     {
       description = (String) transcriptFeature.getValue(DESCRIPTION);
@@ -325,27 +421,65 @@ public class EnsemblGene extends EnsemblSeqProxy
      * transfer features to the new sequence; we use EnsemblCdna to do this,
      * to filter out unwanted features types (see method retainFeature)
      */
-    List<int[]> mapTo = new ArrayList<int[]>();
+    List<int[]> mapTo = new ArrayList<>();
     mapTo.add(new int[] { 1, transcriptLength });
     MapList mapping = new MapList(mappedFrom, mapTo, 1, 1);
-    new EnsemblCdna(getDomain()).transferFeatures(
-            gene.getSequenceFeatures(), transcript.getDatasetSequence(),
-            mapping, parentId);
+    EnsemblCdna cdna = new EnsemblCdna(getDomain());
+    cdna.transferFeatures(gene.getFeatures().getPositionalFeatures(),
+            transcript.getDatasetSequence(), mapping, parentId);
+
+    mapTranscriptToChromosome(transcript, gene, mapping);
 
     /*
      * fetch and save cross-references
      */
-    new EnsemblCdna(getDomain()).getCrossReferences(transcript);
+    cdna.getCrossReferences(transcript);
 
     /*
      * and finally fetch the protein product and save as a cross-reference
      */
-    new EnsemblCdna(getDomain()).addProteinProduct(transcript);
+    cdna.addProteinProduct(transcript);
 
     return transcript;
   }
 
   /**
+   * If the gene has a mapping to chromosome coordinates, derive the transcript
+   * chromosome regions and save on the transcript sequence
+   * 
+   * @param transcript
+   * @param gene
+   * @param mapping
+   *          the mapping from gene to transcript positions
+   */
+  protected void mapTranscriptToChromosome(SequenceI transcript,
+          SequenceI gene, MapList mapping)
+  {
+    GeneLociI loci = gene.getGeneLoci();
+    if (loci == null)
+    {
+      return;
+    }
+
+    MapList geneMapping = loci.getMapping();
+
+    List<int[]> exons = mapping.getFromRanges();
+    List<int[]> transcriptLoci = new ArrayList<>();
+
+    for (int[] exon : exons)
+    {
+      transcriptLoci.add(geneMapping.locateInTo(exon[0], exon[1]));
+    }
+
+    List<int[]> transcriptRange = Arrays.asList(new int[] {
+        transcript.getStart(), transcript.getEnd() });
+    MapList mapList = new MapList(transcriptRange, transcriptLoci, 1, 1);
+
+    transcript.setGeneLoci(loci.getSpeciesId(), loci.getAssemblyId(),
+            loci.getChromosomeId(), mapList);
+  }
+
+  /**
    * Returns the 'transcript_id' property of the sequence feature (or null)
    * 
    * @param feature
@@ -353,12 +487,18 @@ public class EnsemblGene extends EnsemblSeqProxy
    */
   protected String getTranscriptId(SequenceFeature feature)
   {
-    return (String) feature.getValue("transcript_id");
+    return (String) feature.getValue(JSON_ID);
   }
 
   /**
    * Returns a list of the transcript features on the sequence whose Parent is
    * the gene for the accession id.
+   * <p>
+   * Transcript features are those of type "transcript", or any of its sub-types
+   * in the Sequence Ontology e.g. "mRNA", "processed_transcript". We also
+   * include "NMD_transcript_variant", because this type behaves like a
+   * transcript identifier in Ensembl, although strictly speaking it is not in
+   * the SO.
    * 
    * @param accId
    * @param geneSequence
@@ -367,23 +507,21 @@ public class EnsemblGene extends EnsemblSeqProxy
   protected List<SequenceFeature> getTranscriptFeatures(String accId,
           SequenceI geneSequence)
   {
-    List<SequenceFeature> transcriptFeatures = new ArrayList<SequenceFeature>();
+    List<SequenceFeature> transcriptFeatures = new ArrayList<>();
 
-    String parentIdentifier = GENE_PREFIX + accId;
-    SequenceFeature[] sfs = geneSequence.getSequenceFeatures();
+    String parentIdentifier = accId;
 
-    if (sfs != null)
+    List<SequenceFeature> sfs = geneSequence.getFeatures()
+            .getFeaturesByOntology(SequenceOntologyI.TRANSCRIPT);
+    sfs.addAll(geneSequence.getFeatures().getPositionalFeatures(
+            SequenceOntologyI.NMD_TRANSCRIPT_VARIANT));
+
+    for (SequenceFeature sf : sfs)
     {
-      for (SequenceFeature sf : sfs)
+      String parent = (String) sf.getValue(PARENT);
+      if (parentIdentifier.equalsIgnoreCase(parent))
       {
-        if (isTranscript(sf.getType()))
-        {
-          String parent = (String) sf.getValue(PARENT);
-          if (parentIdentifier.equals(parent))
-          {
-            transcriptFeatures.add(sf);
-          }
-        }
+        transcriptFeatures.add(sf);
       }
     }
 
@@ -409,22 +547,26 @@ public class EnsemblGene extends EnsemblSeqProxy
   }
 
   /**
-   * Answers true for a feature of type 'gene' (or a sub-type of gene in the
-   * Sequence Ontology), whose ID is the accession we are retrieving
+   * Answers a list of sequence features (if any) whose type is 'gene' (or a
+   * subtype of gene in the Sequence Ontology), and whose ID is the accession we
+   * are retrieving
    */
   @Override
-  protected boolean identifiesSequence(SequenceFeature sf, String accId)
+  protected List<SequenceFeature> getIdentifyingFeatures(SequenceI seq,
+          String accId)
   {
-    if (SequenceOntologyFactory.getInstance().isA(sf.getType(),
-            SequenceOntologyI.GENE))
+    List<SequenceFeature> result = new ArrayList<>();
+    List<SequenceFeature> sfs = seq.getFeatures()
+            .getFeaturesByOntology(SequenceOntologyI.GENE);
+    for (SequenceFeature sf : sfs)
     {
-      String id = (String) sf.getValue(ID);
-      if ((GENE_PREFIX + accId).equals(id))
+      String id = (String) sf.getValue(JSON_ID);
+      if (accId.equalsIgnoreCase(id))
       {
-        return true;
+        result.add(sf);
       }
     }
-    return false;
+    return result;
   }
 
   /**
@@ -446,7 +588,7 @@ public class EnsemblGene extends EnsemblSeqProxy
     if (isTranscript(type))
     {
       String parent = (String) sf.getValue(PARENT);
-      if (!(GENE_PREFIX + accessionId).equals(parent))
+      if (!accessionId.equalsIgnoreCase(parent))
       {
         return false;
       }
@@ -455,26 +597,6 @@ public class EnsemblGene extends EnsemblSeqProxy
   }
 
   /**
-   * Answers false. This allows an optimisation - a single 'gene' feature is all
-   * that is needed to identify the positions of the gene on the genomic
-   * sequence.
-   */
-  @Override
-  protected boolean isSpliceable()
-  {
-    return false;
-  }
-
-  @Override
-  protected List<String> getCrossReferenceDatabases()
-  {
-    // found these for ENSG00000157764 on 30/01/2016:
-    // return new String[] {"Vega_gene", "OTTG", "ENS_LRG_gene", "ArrayExpress",
-    // "EntrezGene", "HGNC", "MIM_GENE", "MIM_MORBID", "WikiGene"};
-    return CROSS_REFERENCES;
-  }
-
-  /**
    * Override to do nothing as Ensembl doesn't return a protein sequence for a
    * gene identifier
    */
@@ -489,17 +611,28 @@ public class EnsemblGene extends EnsemblSeqProxy
     return ACCESSION_REGEX;
   }
 
+  /**
+   * Returns a descriptor for suitable feature display settings with
+   * <ul>
+   * <li>only exon or sequence_variant features (or their subtypes in the
+   * Sequence Ontology) visible</li>
+   * <li>variant features coloured red</li>
+   * <li>exon features coloured by label (exon name)</li>
+   * <li>variants displayed above (on top of) exons</li>
+   * </ul>
+   */
   @Override
-  public FeatureSettingsI getFeatureColourScheme()
+  public FeatureSettingsModelI getFeatureColourScheme()
   {
     return new FeatureSettingsAdapter()
     {
       SequenceOntologyI so = SequenceOntologyFactory.getInstance();
+
       @Override
       public boolean isFeatureDisplayed(String type)
       {
-        return (so.isA(type, SequenceOntologyI.EXON) || so.isA(type,
-                SequenceOntologyI.SEQUENCE_VARIANT));
+        return (so.isA(type, SequenceOntologyI.EXON)
+                || so.isA(type, SequenceOntologyI.SEQUENCE_VARIANT));
       }
 
       @Override
@@ -507,7 +640,7 @@ public class EnsemblGene extends EnsemblSeqProxy
       {
         if (so.isA(type, SequenceOntologyI.EXON))
         {
-          return new FeatureColourAdapter()
+          return new FeatureColour()
           {
             @Override
             public boolean isColourByLabel()
@@ -518,7 +651,7 @@ public class EnsemblGene extends EnsemblSeqProxy
         }
         if (so.isA(type, SequenceOntologyI.SEQUENCE_VARIANT))
         {
-          return new FeatureColourAdapter()
+          return new FeatureColour()
           {
 
             @Override