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