src/jalview/ext/ensembl/EnsemblGene.java

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