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