JAL-3076 refactor for more efficient scan of 'gene' features
[jalview.git] / src / jalview / ext / ensembl / EnsemblGene.java
index df4e45a..36b19e2 100644 (file)
@@ -1,16 +1,43 @@
+/*
+ * 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.FeatureSettingsI;
+import jalview.api.FeatureColourI;
+import jalview.api.FeatureSettingsModelI;
 import jalview.datamodel.AlignmentI;
+import jalview.datamodel.DBRefEntry;
+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.FeatureColourScheme;
+import jalview.schemes.FeatureColour;
+import jalview.schemes.FeatureSettingsAdapter;
 import jalview.util.MapList;
-import jalview.util.StringUtils;
 
+import java.awt.Color;
+import java.io.UnsupportedEncodingException;
+import java.net.URLDecoder;
 import java.util.ArrayList;
 import java.util.Arrays;
 import java.util.List;
@@ -37,10 +64,28 @@ public class EnsemblGene extends EnsemblSeqProxy
       EnsemblFeatureType.exon, EnsemblFeatureType.cds,
       EnsemblFeatureType.variation };
 
+  /**
+   * Default constructor (to use rest.ensembl.org)
+   */
+  public EnsemblGene()
+  {
+    super();
+  }
+
+  /**
+   * Constructor given the target domain to fetch data from
+   * 
+   * @param d
+   */
+  public EnsemblGene(String d)
+  {
+    super(d);
+  }
+
   @Override
   public String getDbName()
   {
-    return "ENSEMBL (GENE)";
+    return "ENSEMBL";
   }
 
   @Override
@@ -55,6 +100,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
@@ -67,7 +118,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>
@@ -78,79 +130,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().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().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;
+  }
 
-    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.getMap());
     }
+    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(DBRefEntry.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().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 (!geneIds.contains(geneId))
+        {
+          geneIds.add(geneId);
+        }
+      }
+      else
       {
-        if (isGeneIdentifier(id))
+        /*
+         * 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)
         {
-          return id;
+          if (!geneIds.contains(id))
+          {
+            geneIds.add(id);
+          }
         }
       }
     }
-    return null;
+    return geneIds;
   }
 
   /**
@@ -173,6 +300,33 @@ public class EnsemblGene extends EnsemblSeqProxy
     {
       makeTranscript(transcriptFeature, al, gene);
     }
+
+    clearGeneFeatures(gene);
+  }
+
+  /**
+   * Remove unwanted features (transcript, exon, CDS) from the gene sequence
+   * after we have used them to derive transcripts and transfer features
+   * 
+   * @param gene
+   */
+  protected void clearGeneFeatures(SequenceI gene)
+  {
+    /*
+     * 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)
+    {
+      gene.deleteFeature(sf);
+    }
   }
 
   /**
@@ -188,8 +342,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)
@@ -221,11 +375,12 @@ public class EnsemblGene extends EnsemblSeqProxy
     {
       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)
     {
@@ -237,11 +392,27 @@ public class EnsemblGene extends EnsemblSeqProxy
       mappedFrom.add(new int[] { sf.getBegin(), sf.getEnd() });
     }
 
-    Sequence transcript = new Sequence(accId, seqChars, 1, transcriptLength);
-    String geneName = (String) transcriptFeature.getValue(NAME);
-    if (geneName != null)
+    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);
+    if (description == null)
     {
-      transcript.setDescription(geneName);
+      description = (String) transcriptFeature.getValue(DESCRIPTION);
+    }
+    if (description != null)
+    {
+      try
+      {
+        transcript.setDescription(URLDecoder.decode(description, "UTF-8"));
+      } catch (UnsupportedEncodingException e)
+      {
+        e.printStackTrace(); // as if
+      }
     }
     transcript.createDatasetSequence();
 
@@ -251,26 +422,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().transferFeatures(gene.getSequenceFeatures(),
+    EnsemblCdna cdna = new EnsemblCdna(getDomain());
+    cdna.transferFeatures(gene.getFeatures().getPositionalFeatures(),
             transcript.getDatasetSequence(), mapping, parentId);
 
+    mapTranscriptToChromosome(transcript, gene, mapping);
+
     /*
      * fetch and save cross-references
      */
-    super.getCrossReferences(transcript);
+    cdna.getCrossReferences(transcript);
 
     /*
      * and finally fetch the protein product and save as a cross-reference
      */
-    new EnsemblCdna().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.getMap();
+
+    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
@@ -284,6 +494,12 @@ public class EnsemblGene extends EnsemblSeqProxy
   /**
    * 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
@@ -292,23 +508,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();
 
-    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);
       }
     }
 
@@ -334,22 +548,27 @@ 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))
+      // NB features as gff use 'ID'; rest services return as 'id'
+      String id = (String) sf.getValue("ID");
+      if ((GENE_PREFIX + accId).equalsIgnoreCase(id))
       {
-        return true;
+        result.add(sf);
       }
     }
-    return false;
+    return result;
   }
 
   /**
@@ -362,16 +581,16 @@ public class EnsemblGene extends EnsemblSeqProxy
   @Override
   protected boolean retainFeature(SequenceFeature sf, String accessionId)
   {
-    if (SequenceOntologyFactory.getInstance().isA(sf.getType(),
-            SequenceOntologyI.GENE))
+    SequenceOntologyI so = SequenceOntologyFactory.getInstance();
+    String type = sf.getType();
+    if (so.isA(type, SequenceOntologyI.GENE))
     {
       return false;
     }
-
-    if (isTranscript(sf.getType()))
+    if (isTranscript(type))
     {
       String parent = (String) sf.getValue(PARENT);
-      if (!(GENE_PREFIX + accessionId).equals(parent))
+      if (!(GENE_PREFIX + accessionId).equalsIgnoreCase(parent))
       {
         return false;
       }
@@ -380,26 +599,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 super.getCrossReferenceDatabases();
-  }
-
-  /**
    * Override to do nothing as Ensembl doesn't return a protein sequence for a
    * gene identifier
    */
@@ -414,10 +613,84 @@ 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 FeatureColourScheme.EnsemblVariants;
+    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));
+      }
+
+      @Override
+      public FeatureColourI getFeatureColour(String type)
+      {
+        if (so.isA(type, SequenceOntologyI.EXON))
+        {
+          return new FeatureColour()
+          {
+            @Override
+            public boolean isColourByLabel()
+            {
+              return true;
+            }
+          };
+        }
+        if (so.isA(type, SequenceOntologyI.SEQUENCE_VARIANT))
+        {
+          return new FeatureColour()
+          {
+
+            @Override
+            public Color getColour()
+            {
+              return Color.RED;
+            }
+          };
+        }
+        return null;
+      }
+
+      /**
+       * order to render sequence_variant after exon after the rest
+       */
+      @Override
+      public int compare(String feature1, String feature2)
+      {
+        if (so.isA(feature1, SequenceOntologyI.SEQUENCE_VARIANT))
+        {
+          return +1;
+        }
+        if (so.isA(feature2, SequenceOntologyI.SEQUENCE_VARIANT))
+        {
+          return -1;
+        }
+        if (so.isA(feature1, SequenceOntologyI.EXON))
+        {
+          return +1;
+        }
+        if (so.isA(feature2, SequenceOntologyI.EXON))
+        {
+          return -1;
+        }
+        return 0;
+      }
+    };
   }
 
 }