JAL-2738 use GeneLocus extends DBRefEntry to hold chromosomal mappings
[jalview.git] / src / jalview / ext / ensembl / EnsemblGene.java
index dc28796..9c7ce63 100644 (file)
@@ -1,13 +1,42 @@
+/*
+ * 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.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.FeatureColour;
+import jalview.schemes.FeatureSettingsAdapter;
 import jalview.util.MapList;
 
+import java.awt.Color;
+import java.io.UnsupportedEncodingException;
+import java.net.URLDecoder;
 import java.util.ArrayList;
 import java.util.Arrays;
 import java.util.List;
@@ -21,19 +50,41 @@ import com.stevesoft.pat.Regex;
  */
 public class EnsemblGene extends EnsemblSeqProxy
 {
-  // TODO modify to accept other species e.g. ENSMUSGnnn
-  private static final Regex ACCESSION_REGEX = new Regex(
-          "((ENSG)[0-9]{11})");
+  /*
+   * accepts anything as we will attempt lookup of gene or 
+   * transcript id or gene name
+   */
+  private static final Regex ACCESSION_REGEX = new Regex(".*");
 
   private static final EnsemblFeatureType[] FEATURES_TO_FETCH = {
       EnsemblFeatureType.gene, EnsemblFeatureType.transcript,
       EnsemblFeatureType.exon, EnsemblFeatureType.cds,
       EnsemblFeatureType.variation };
 
+  private static final String CHROMOSOME = "chromosome";
+
+  /**
+   * 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
@@ -48,12 +99,26 @@ public class EnsemblGene extends EnsemblSeqProxy
     return EnsemblSeqType.GENOMIC;
   }
 
+  @Override
+  protected String getObjectType()
+  {
+    return OBJECT_TYPE_GENE;
+  }
+
   /**
-   * Builds an alignment of all transcripts for the requested 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
+   * gene name or external identifier, returns any related gene sequences found
+   * for model organisms. If only a single gene is queried for, then its
+   * transcripts are also retrieved and added to the alignment. <br>
+   * Method:
    * <ul>
+   * <li>resolves a transcript identifier by looking up its parent gene id</li>
+   * <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>
@@ -62,21 +127,159 @@ public class EnsemblGene extends EnsemblSeqProxy
    * <li>aligns each transcript against the gene sequence based on the position
    * mappings</li>
    * </ul>
+   * 
+   * @param query
+   *          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 ? if an ENST identifier is supplied, convert to ENSG?
-    AlignmentI al = super.getSequenceRecords(query);
-    if (al.getHeight() > 0)
+    /*
+     * convert to a non-duplicated list of gene identifiers
+     */
+    List<String> geneIds = getGeneIds(query);
+
+    AlignmentI al = null;
+    for (String geneId : geneIds)
     {
-      getTranscripts(al, query);
-    }
+      /*
+       * fetch the gene sequence(s) with features and xrefs
+       */
+      AlignmentI geneAlignment = super.getSequenceRecords(geneId);
+      if (geneAlignment == null)
+      {
+        continue;
+      }
+      
+      if (geneAlignment.getHeight() == 1)
+      {
+        // ensure id has 'correct' case for the Ensembl identifier
+        geneId = geneAlignment.getSequenceAt(0).getName();
 
+        findGeneLoci(geneAlignment.getSequenceAt(0), geneId);
+
+        getTranscripts(geneAlignment, geneId);
+      }
+      if (al == null)
+      {
+        al = geneAlignment;
+      }
+      else
+      {
+        al.append(geneAlignment);
+      }
+    }
     return al;
   }
 
   /**
+   * 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)
+    {
+      seq.setGeneLoci(geneLoci.getSpeciesId(), geneLoci.getAssemblyId(),
+              geneLoci.getChromosomeId(), geneLoci.getMapping());
+    }
+    else
+    {
+      parseChromosomeLocations(seq);
+    }
+  }
+
+  /**
+   * 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;
+  }
+
+  /**
+   * 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 accessions
+   * @return
+   */
+  List<String> getGeneIds(String accessions)
+  {
+    List<String> geneIds = new ArrayList<>();
+
+    for (String acc : accessions.split(getAccessionSeparator()))
+    {
+      /*
+       * 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 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 geneIds;
+  }
+
+  /**
    * Constructs all transcripts for the gene, as identified by "transcript"
    * features whose Parent is the requested gene. The coding transcript
    * sequences (i.e. with introns omitted) are added to the alignment.
@@ -96,6 +299,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);
+    }
   }
 
   /**
@@ -111,10 +341,10 @@ 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 = (String) transcriptFeature.getValue("transcript_id");
+    String accId = getTranscriptId(transcriptFeature);
     if (accId == null)
     {
       return null;
@@ -137,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)
     {
@@ -160,11 +391,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 = transcriptFeature.getDescription();
+    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();
 
@@ -174,23 +421,84 @@ 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
+     */
+    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.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
+   * @return
+   */
+  protected String getTranscriptId(SequenceFeature feature)
+  {
+    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
@@ -199,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:" + 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);
       }
     }
 
@@ -225,11 +531,11 @@ public class EnsemblGene extends EnsemblSeqProxy
   @Override
   public String getDescription()
   {
-    return "Fetches all transcripts and variant features for a gene";
+    return "Fetches all transcripts and variant features for a gene or transcript";
   }
 
   /**
-   * Default test query is a transcript
+   * Default test query is a gene id (can also enter a transcript id)
    */
   @Override
   public String getTestQuery()
@@ -241,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:" + accId).equals(id))
+      String id = (String) sf.getValue(JSON_ID);
+      if (accId.equalsIgnoreCase(id))
       {
-        return true;
+        result.add(sf);
       }
     }
-    return false;
+    return result;
   }
 
   /**
@@ -269,16 +579,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:" + accessionId).equals(parent))
+      if (!accessionId.equalsIgnoreCase(parent))
       {
         return false;
       }
@@ -287,38 +597,98 @@ 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 to do nothing as Ensembl doesn't return a protein sequence for a
+   * gene identifier
    */
   @Override
-  protected boolean isSpliceable()
+  protected void addProteinProduct(SequenceI querySeq)
   {
-    return false;
   }
 
   @Override
-  protected List<String> getCrossReferenceDatabases()
+  public Regex getAccessionValidator()
   {
-    // 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();
+    return ACCESSION_REGEX;
   }
 
   /**
-   * Override to do nothing as Ensembl doesn't return a protein sequence for a
-   * gene identifier
+   * 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
-  protected void addProteinProduct(SequenceI querySeq)
+  public FeatureSettingsModelI getFeatureColourScheme()
   {
-  }
+    return new FeatureSettingsAdapter()
+    {
+      SequenceOntologyI so = SequenceOntologyFactory.getInstance();
 
-  @Override
-  public Regex getAccessionValidator()
-  {
-    return ACCESSION_REGEX;
+      @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;
+      }
+    };
   }
 
 }