src/jalview/ext/ensembl/EnsemblGene.java

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