src/jalview/ext/ensembl/EnsemblGene.java

   1 package jalview.ext.ensembl;
   2
   3 import jalview.api.FeatureColourI;
   4 import jalview.api.FeatureSettingsI;
   5 import jalview.datamodel.AlignmentI;
   6 import jalview.datamodel.Sequence;
   7 import jalview.datamodel.SequenceFeature;
   8 import jalview.datamodel.SequenceI;
   9 import jalview.io.gff.SequenceOntologyFactory;
  10 import jalview.io.gff.SequenceOntologyI;
  11 import jalview.schemes.FeatureColourAdapter;
  12 import jalview.schemes.FeatureSettingsAdapter;
  13 import jalview.util.MapList;
  14 import jalview.util.StringUtils;
  15
  16 import java.awt.Color;
  17 import java.io.UnsupportedEncodingException;
  18 import java.net.URLDecoder;
  19 import java.util.ArrayList;
  20 import java.util.Arrays;
  21 import java.util.List;
  22
  23 import com.stevesoft.pat.Regex;
  24
  25 /**
  26  * A class that fetches genomic sequence and all transcripts for an Ensembl gene
  27  *
  28  * @author gmcarstairs
  29  */
  30 public class EnsemblGene extends EnsemblSeqProxy
  31 {
  32   private static final List<String> CROSS_REFERENCES = Arrays
  33           .asList(new String[] { "CCDS" });
  34
  35   private static final String GENE_PREFIX = "gene:";
  36
  37   /*
  38    * accepts anything as we will attempt lookup of gene or
  39    * transcript id or gene name
  40    */
  41   private static final Regex ACCESSION_REGEX = new Regex(".*");
  42
  43   private static final EnsemblFeatureType[] FEATURES_TO_FETCH = {
  44       EnsemblFeatureType.gene, EnsemblFeatureType.transcript,
  45       EnsemblFeatureType.exon, EnsemblFeatureType.cds,
  46       EnsemblFeatureType.variation };
  47
  48   /**
  49    * Default constructor (to use rest.ensembl.org)
  50    */
  51   public EnsemblGene()
  52   {
  53     super();
  54   }
  55
  56   /**
  57    * Constructor given the target domain to fetch data from
  58    *
  59    * @param d
  60    */
  61   public EnsemblGene(String d)
  62   {
  63     super(d);
  64   }
  65
  66   @Override
  67   public String getDbName()
  68   {
  69     return "ENSEMBL";
  70   }
  71
  72   @Override
  73   protected EnsemblFeatureType[] getFeaturesToFetch()
  74   {
  75     return FEATURES_TO_FETCH;
  76   }
  77
  78   @Override
  79   protected EnsemblSeqType getSourceEnsemblType()
  80   {
  81     return EnsemblSeqType.GENOMIC;
  82   }
  83
  84   /**
  85    * Returns an alignment containing the gene(s) for the given gene or
  86    * transcript identifier, or external identifier (e.g. Uniprot id). If given a
  87    * gene name or external identifier, returns any related gene sequences found
  88    * for model organisms. If only a single gene is queried for, then its
  89    * transcripts are also retrieved and added to the alignment. <br>
  90    * Method:
  91    * <ul>
  92    * <li>resolves a transcript identifier by looking up its parent gene id</li>
  93    * <li>resolves an external identifier by looking up xref-ed gene ids</li>
  94    * <li>fetches the gene sequence</li>
  95    * <li>fetches features on the sequence</li>
  96    * <li>identifies "transcript" features whose Parent is the requested gene</li>
  97    * <li>fetches the transcript sequence for each transcript</li>
  98    * <li>makes a mapping from the gene to each transcript</li>
  99    * <li>copies features from gene to transcript sequences</li>
 100    * <li>fetches the protein sequence for each transcript, maps and saves it as
 101    * a cross-reference</li>
 102    * <li>aligns each transcript against the gene sequence based on the position
 103    * mappings</li>
 104    * </ul>
 105    *
 106    * @param query
 107    *          one or more identifiers separated by a space
 108    * @return an alignment containing one or more genes, and possibly
 109    *         transcripts, or null
 110    */
 111   @Override
 112   public AlignmentI getSequenceRecords(String query) throws Exception
 113   {
 114     // todo: tidy up handling of one or multiple accession ids
 115     String[] queries = query.split(getAccessionSeparator());
 116
 117     /*
 118      * if given a transcript id, look up its gene parent
 119      */
 120     if (isTranscriptIdentifier(query))
 121     {
 122       // we are assuming all transcripts have the same gene parent here
 123       query = new EnsemblLookup(getDomain()).getParent(queries[0]);
 124       if (query == null)
 125       {
 126         return null;
 127       }
 128     }
 129
 130     /*
 131      * if given a gene or other external name, lookup and fetch
 132      * the corresponding gene for all model organisms
 133      */
 134     if (!isGeneIdentifier(query))
 135     {
 136       List<String> geneIds = new EnsemblSymbol(getDomain()).getIds(query);
 137       if (geneIds.isEmpty())
 138       {
 139         return null;
 140       }
 141       String theIds = StringUtils.listToDelimitedString(geneIds,
 142               getAccessionSeparator());
 143       return getSequenceRecords(theIds);
 144     }
 145
 146     /*
 147      * fetch the gene sequence(s) with features and xrefs
 148      */
 149     AlignmentI al = super.getSequenceRecords(query);
 150
 151     /*
 152      * if we retrieved a single gene, get its transcripts as well
 153      */
 154     if (al.getHeight() == 1)
 155     {
 156       getTranscripts(al, query);
 157     }
 158
 159     return al;
 160   }
 161
 162   /**
 163    * Attempts to get Ensembl stable identifiers for model organisms for a gene
 164    * name by calling the xrefs symbol REST service to resolve the gene name.
 165    *
 166    * @param query
 167    * @return
 168    */
 169   protected String getGeneIdentifiersForName(String query)
 170   {
 171     List<String> ids = new EnsemblSymbol(getDomain()).getIds(query);
 172     if (ids != null)
 173     {
 174       for (String id : ids)
 175       {
 176         if (isGeneIdentifier(id))
 177         {
 178           return id;
 179         }
 180       }
 181     }
 182     return null;
 183   }
 184
 185   /**
 186    * Constructs all transcripts for the gene, as identified by "transcript"
 187    * features whose Parent is the requested gene. The coding transcript
 188    * sequences (i.e. with introns omitted) are added to the alignment.
 189    *
 190    * @param al
 191    * @param accId
 192    * @throws Exception
 193    */
 194   protected void getTranscripts(AlignmentI al, String accId)
 195           throws Exception
 196   {
 197     SequenceI gene = al.getSequenceAt(0);
 198     List<SequenceFeature> transcriptFeatures = getTranscriptFeatures(accId,
 199             gene);
 200
 201     for (SequenceFeature transcriptFeature : transcriptFeatures)
 202     {
 203       makeTranscript(transcriptFeature, al, gene);
 204     }
 205
 206     clearGeneFeatures(gene);
 207   }
 208
 209   /**
 210    * Remove unwanted features (transcript, exon, CDS) from the gene sequence
 211    * after we have used them to derive transcripts and transfer features
 212    *
 213    * @param gene
 214    */
 215   protected void clearGeneFeatures(SequenceI gene)
 216   {
 217     SequenceFeature[] sfs = gene.getSequenceFeatures();
 218     if (sfs != null)
 219     {
 220       SequenceOntologyI so = SequenceOntologyFactory.getInstance();
 221       List<SequenceFeature> filtered = new ArrayList<SequenceFeature>();
 222       for (SequenceFeature sf : sfs)
 223       {
 224         String type = sf.getType();
 225         if (!isTranscript(type) && !so.isA(type, SequenceOntologyI.EXON)
 226                 && !so.isA(type, SequenceOntologyI.CDS))
 227         {
 228           filtered.add(sf);
 229         }
 230       }
 231       gene.setSequenceFeatures(filtered
 232               .toArray(new SequenceFeature[filtered
 233               .size()]));
 234     }
 235   }
 236
 237   /**
 238    * Constructs a spliced transcript sequence by finding 'exon' features for the
 239    * given id (or failing that 'CDS'). Copies features on to the new sequence.
 240    * 'Aligns' the new sequence against the gene sequence by padding with gaps,
 241    * and adds it to the alignment.
 242    *
 243    * @param transcriptFeature
 244    * @param al
 245    *          the alignment to which to add the new sequence
 246    * @param gene
 247    *          the parent gene sequence, with features
 248    * @return
 249    */
 250   SequenceI makeTranscript(SequenceFeature transcriptFeature,
 251           AlignmentI al, SequenceI gene)
 252   {
 253     String accId = getTranscriptId(transcriptFeature);
 254     if (accId == null)
 255     {
 256       return null;
 257     }
 258
 259     /*
 260      * NB we are mapping from gene sequence (not genome), so do not
 261      * need to check for reverse strand (gene and transcript sequences
 262      * are in forward sense)
 263      */
 264
 265     /*
 266      * make a gene-length sequence filled with gaps
 267      * we will fill in the bases for transcript regions
 268      */
 269     char[] seqChars = new char[gene.getLength()];
 270     Arrays.fill(seqChars, al.getGapCharacter());
 271
 272     /*
 273      * look for exon features of the transcript, failing that for CDS
 274      * (for example ENSG00000124610 has 1 CDS but no exon features)
 275      */
 276     String parentId = "transcript:" + accId;
 277     List<SequenceFeature> splices = findFeatures(gene,
 278             SequenceOntologyI.EXON, parentId);
 279     if (splices.isEmpty())
 280     {
 281       splices = findFeatures(gene, SequenceOntologyI.CDS, parentId);
 282     }
 283
 284     int transcriptLength = 0;
 285     final char[] geneChars = gene.getSequence();
 286     int offset = gene.getStart(); // to convert to 0-based positions
 287     List<int[]> mappedFrom = new ArrayList<int[]>();
 288
 289     for (SequenceFeature sf : splices)
 290     {
 291       int start = sf.getBegin() - offset;
 292       int end = sf.getEnd() - offset;
 293       int spliceLength = end - start + 1;
 294       System.arraycopy(geneChars, start, seqChars, start, spliceLength);
 295       transcriptLength += spliceLength;
 296       mappedFrom.add(new int[] { sf.getBegin(), sf.getEnd() });
 297     }
 298
 299     Sequence transcript = new Sequence(accId, seqChars, 1, transcriptLength);
 300
 301     /*
 302      * Ensembl has gene name as transcript Name
 303      * EnsemblGenomes doesn't, but has a url-encoded description field
 304      */
 305     String description = (String) transcriptFeature.getValue(NAME);
 306     if (description == null)
 307     {
 308       description = (String) transcriptFeature.getValue(DESCRIPTION);
 309     }
 310     if (description != null)
 311     {
 312       try
 313       {
 314         transcript.setDescription(URLDecoder.decode(description, "UTF-8"));
 315       } catch (UnsupportedEncodingException e)
 316       {
 317         e.printStackTrace(); // as if
 318       }
 319     }
 320     transcript.createDatasetSequence();
 321
 322     al.addSequence(transcript);
 323
 324     /*
 325      * transfer features to the new sequence; we use EnsemblCdna to do this,
 326      * to filter out unwanted features types (see method retainFeature)
 327      */
 328     List<int[]> mapTo = new ArrayList<int[]>();
 329     mapTo.add(new int[] { 1, transcriptLength });
 330     MapList mapping = new MapList(mappedFrom, mapTo, 1, 1);
 331     new EnsemblCdna(getDomain()).transferFeatures(
 332             gene.getSequenceFeatures(), transcript.getDatasetSequence(),
 333             mapping, parentId);
 334
 335     /*
 336      * fetch and save cross-references
 337      */
 338     new EnsemblCdna(getDomain()).getCrossReferences(transcript);
 339
 340     /*
 341      * and finally fetch the protein product and save as a cross-reference
 342      */
 343     new EnsemblCdna(getDomain()).addProteinProduct(transcript);
 344
 345     return transcript;
 346   }
 347
 348   /**
 349    * Returns the 'transcript_id' property of the sequence feature (or null)
 350    *
 351    * @param feature
 352    * @return
 353    */
 354   protected String getTranscriptId(SequenceFeature feature)
 355   {
 356     return (String) feature.getValue("transcript_id");
 357   }
 358
 359   /**
 360    * Returns a list of the transcript features on the sequence whose Parent is
 361    * the gene for the accession id.
 362    *
 363    * @param accId
 364    * @param geneSequence
 365    * @return
 366    */
 367   protected List<SequenceFeature> getTranscriptFeatures(String accId,
 368           SequenceI geneSequence)
 369   {
 370     List<SequenceFeature> transcriptFeatures = new ArrayList<SequenceFeature>();
 371
 372     String parentIdentifier = GENE_PREFIX + accId;
 373     SequenceFeature[] sfs = geneSequence.getSequenceFeatures();
 374
 375     if (sfs != null)
 376     {
 377       for (SequenceFeature sf : sfs)
 378       {
 379         if (isTranscript(sf.getType()))
 380         {
 381           String parent = (String) sf.getValue(PARENT);
 382           if (parentIdentifier.equals(parent))
 383           {
 384             transcriptFeatures.add(sf);
 385           }
 386         }
 387       }
 388     }
 389
 390     return transcriptFeatures;
 391   }
 392
 393   @Override
 394   public String getDescription()
 395   {
 396     return "Fetches all transcripts and variant features for a gene or transcript";
 397   }
 398
 399   /**
 400    * Default test query is a gene id (can also enter a transcript id)
 401    */
 402   @Override
 403   public String getTestQuery()
 404   {
 405     return "ENSG00000157764"; // BRAF, 5 transcripts, reverse strand
 406     // ENSG00000090266 // NDUFB2, 15 transcripts, forward strand
 407     // ENSG00000101812 // H2BFM histone, 3 transcripts, forward strand
 408     // ENSG00000123569 // H2BFWT histone, 2 transcripts, reverse strand
 409   }
 410
 411   /**
 412    * Answers true for a feature of type 'gene' (or a sub-type of gene in the
 413    * Sequence Ontology), whose ID is the accession we are retrieving
 414    */
 415   @Override
 416   protected boolean identifiesSequence(SequenceFeature sf, String accId)
 417   {
 418     if (SequenceOntologyFactory.getInstance().isA(sf.getType(),
 419             SequenceOntologyI.GENE))
 420     {
 421       String id = (String) sf.getValue(ID);
 422       if ((GENE_PREFIX + accId).equals(id))
 423       {
 424         return true;
 425       }
 426     }
 427     return false;
 428   }
 429
 430   /**
 431    * Answers true unless feature type is 'gene', or 'transcript' with a parent
 432    * which is a different gene. We need the gene features to identify the range,
 433    * but it is redundant information on the gene sequence. Checking the parent
 434    * allows us to drop transcript features which belong to different
 435    * (overlapping) genes.
 436    */
 437   @Override
 438   protected boolean retainFeature(SequenceFeature sf, String accessionId)
 439   {
 440     SequenceOntologyI so = SequenceOntologyFactory.getInstance();
 441     String type = sf.getType();
 442     if (so.isA(type, SequenceOntologyI.GENE))
 443     {
 444       return false;
 445     }
 446     if (isTranscript(type))
 447     {
 448       String parent = (String) sf.getValue(PARENT);
 449       if (!(GENE_PREFIX + accessionId).equals(parent))
 450       {
 451         return false;
 452       }
 453     }
 454     return true;
 455   }
 456
 457   /**
 458    * Answers false. This allows an optimisation - a single 'gene' feature is all
 459    * that is needed to identify the positions of the gene on the genomic
 460    * sequence.
 461    */
 462   @Override
 463   protected boolean isSpliceable()
 464   {
 465     return false;
 466   }
 467
 468   @Override
 469   protected List<String> getCrossReferenceDatabases()
 470   {
 471     // found these for ENSG00000157764 on 30/01/2016:
 472     // return new String[] {"Vega_gene", "OTTG", "ENS_LRG_gene", "ArrayExpress",
 473     // "EntrezGene", "HGNC", "MIM_GENE", "MIM_MORBID", "WikiGene"};
 474     return CROSS_REFERENCES;
 475   }
 476
 477   /**
 478    * Override to do nothing as Ensembl doesn't return a protein sequence for a
 479    * gene identifier
 480    */
 481   @Override
 482   protected void addProteinProduct(SequenceI querySeq)
 483   {
 484   }
 485
 486   @Override
 487   public Regex getAccessionValidator()
 488   {
 489     return ACCESSION_REGEX;
 490   }
 491
 492   @Override
 493   public FeatureSettingsI getFeatureColourScheme()
 494   {
 495     return new FeatureSettingsAdapter()
 496     {
 497       SequenceOntologyI so = SequenceOntologyFactory.getInstance();
 498       @Override
 499       public boolean isFeatureDisplayed(String type)
 500       {
 501         return (so.isA(type, SequenceOntologyI.EXON) || so.isA(type,
 502                 SequenceOntologyI.SEQUENCE_VARIANT));
 503       }
 504
 505       @Override
 506       public FeatureColourI getFeatureColour(String type)
 507       {
 508         if (so.isA(type, SequenceOntologyI.EXON))
 509         {
 510           return new FeatureColourAdapter()
 511           {
 512             @Override
 513             public boolean isColourByLabel()
 514             {
 515               return true;
 516             }
 517           };
 518         }
 519         if (so.isA(type, SequenceOntologyI.SEQUENCE_VARIANT))
 520         {
 521           return new FeatureColourAdapter()
 522           {
 523
 524             @Override
 525             public Color getColour()
 526             {
 527               return Color.RED;
 528             }
 529           };
 530         }
 531         return null;
 532       }
 533
 534       /**
 535        * order to render sequence_variant after exon after the rest
 536        */
 537       @Override
 538       public int compare(String feature1, String feature2)
 539       {
 540         if (so.isA(feature1, SequenceOntologyI.SEQUENCE_VARIANT))
 541         {
 542           return +1;
 543         }
 544         if (so.isA(feature2, SequenceOntologyI.SEQUENCE_VARIANT))
 545         {
 546           return -1;
 547         }
 548         if (so.isA(feature1, SequenceOntologyI.EXON))
 549         {
 550           return +1;
 551         }
 552         if (so.isA(feature2, SequenceOntologyI.EXON))
 553         {
 554           return -1;
 555         }
 556         return 0;
 557       }
 558     };
 559   }
 560
 561 }