src/jalview/ext/ensembl/EnsemblGene.java

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