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
  11 import java.util.ArrayList;
  12 import java.util.Arrays;
  13 import java.util.List;
  14
  15 /**
  16  * A class that fetches genomic sequence and all transcripts for an Ensembl gene
  17  *
  18  * @author gmcarstairs
  19  */
  20 public class EnsemblGene extends EnsemblSeqProxy
  21 {
  22   private static final EnsemblFeatureType[] FEATURES_TO_FETCH = {
  23       EnsemblFeatureType.gene, EnsemblFeatureType.transcript,
  24       EnsemblFeatureType.exon, EnsemblFeatureType.cds,
  25       EnsemblFeatureType.variation };
  26
  27   @Override
  28   public String getDbName()
  29   {
  30     return "ENSEMBL (GENE)";
  31   }
  32
  33   @Override
  34   protected EnsemblFeatureType[] getFeaturesToFetch()
  35   {
  36     return FEATURES_TO_FETCH;
  37   }
  38
  39   @Override
  40   protected EnsemblSeqType getSourceEnsemblType()
  41   {
  42     return EnsemblSeqType.GENOMIC;
  43   }
  44
  45   /**
  46    * Builds an alignment of all transcripts for the requested gene:
  47    * <ul>
  48    * <li>fetches the gene sequence</li>
  49    * <li>fetches features on the sequence</li>
  50    * <li>identifies "transcript" features whose Parent is the requested gene</li>
  51    * <li>fetches the transcript sequence for each transcript</li>
  52    * <li>makes a mapping from the gene to each transcript</li>
  53    * <li>copies features from gene to transcript sequences</li>
  54    * <li>fetches the protein sequence for each transcript, maps and saves it as
  55    * a cross-reference</li>
  56    * <li>aligns each transcript against the gene sequence based on the position
  57    * mappings</li>
  58    * </ul>
  59    */
  60   @Override
  61   public AlignmentI getSequenceRecords(String query) throws Exception
  62   {
  63     // TODO ? if an ENST identifier is supplied, convert to ENSG?
  64     AlignmentI al = super.getSequenceRecords(query);
  65     if (al.getHeight() > 0)
  66     {
  67       getTranscripts(al, query);
  68     }
  69
  70     return al;
  71   }
  72
  73   /**
  74    * Constructs all transcripts for the gene, as identified by "transcript"
  75    * features whose Parent is the requested gene. The coding transcript
  76    * sequences (i.e. with introns omitted) are added to the alignment.
  77    *
  78    * @param al
  79    * @param accId
  80    * @throws Exception
  81    */
  82   protected void getTranscripts(AlignmentI al, String accId)
  83           throws Exception
  84   {
  85     SequenceI gene = al.getSequenceAt(0);
  86     List<SequenceFeature> transcriptFeatures = getTranscriptFeatures(accId,
  87             gene);
  88
  89     for (SequenceFeature transcriptFeature : transcriptFeatures)
  90     {
  91       makeTranscript(transcriptFeature, al, gene);
  92     }
  93   }
  94
  95   /**
  96    * Constructs a spliced transcript sequence by finding 'exon' features for the
  97    * given id (or failing that 'CDS'). Copies features on to the new sequence.
  98    * 'Aligns' the new sequence against the gene sequence by padding with gaps,
  99    * and adds it to the alignment.
 100    *
 101    * @param transcriptFeature
 102    * @param al
 103    *          the alignment to which to add the new sequence
 104    * @param gene
 105    *          the parent gene sequence, with features
 106    * @return
 107    */
 108   SequenceI makeTranscript(SequenceFeature transcriptFeature,
 109           AlignmentI al, SequenceI gene)
 110   {
 111     String accId = (String) transcriptFeature.getValue("transcript_id");
 112     if (accId == null)
 113     {
 114       return null;
 115     }
 116
 117     /*
 118      * NB we are mapping from gene sequence (not genome), so do not
 119      * need to check for reverse strand (gene and transcript sequences
 120      * are in forward sense)
 121      */
 122
 123     /*
 124      * make a gene-length sequence filled with gaps
 125      * we will fill in the bases for transcript regions
 126      */
 127     char[] seqChars = new char[gene.getLength()];
 128     Arrays.fill(seqChars, al.getGapCharacter());
 129
 130     /*
 131      * look for exon features of the transcript, failing that for CDS
 132      * (for example ENSG00000124610 has 1 CDS but no exon features)
 133      */
 134     String parentId = "transcript:" + accId;
 135     List<SequenceFeature> splices = findFeatures(gene,
 136             SequenceOntologyI.EXON, parentId);
 137     if (splices.isEmpty())
 138     {
 139       splices = findFeatures(gene, SequenceOntologyI.CDS, parentId);
 140     }
 141
 142     int transcriptLength = 0;
 143     final char[] geneChars = gene.getSequence();
 144     int offset = gene.getStart(); // to convert to 0-based positions
 145     List<int[]> mappedFrom = new ArrayList<int[]>();
 146
 147     for (SequenceFeature sf : splices)
 148     {
 149       int start = sf.getBegin() - offset;
 150       int end = sf.getEnd() - offset;
 151       int spliceLength = end - start + 1;
 152       System.arraycopy(geneChars, start, seqChars, start, spliceLength);
 153       transcriptLength += spliceLength;
 154       mappedFrom.add(new int[] { sf.getBegin(), sf.getEnd() });
 155     }
 156
 157     Sequence transcript = new Sequence(accId, seqChars, 1, transcriptLength);
 158     String geneName = (String) transcriptFeature.getValue(NAME);
 159     if (geneName != null)
 160     {
 161       transcript.setDescription(geneName);
 162     }
 163     transcript.createDatasetSequence();
 164
 165     al.addSequence(transcript);
 166
 167     /*
 168      * transfer features to the new sequence; we use EnsemblCdna to do this,
 169      * to filter out unwanted features types (see method retainFeature)
 170      */
 171     List<int[]> mapTo = new ArrayList<int[]>();
 172     mapTo.add(new int[] { 1, transcriptLength });
 173     MapList mapping = new MapList(mappedFrom, mapTo, 1, 1);
 174     new EnsemblCdna().transferFeatures(gene.getSequenceFeatures(),
 175             transcript.getDatasetSequence(), mapping, parentId);
 176
 177     /*
 178      * and finally fetch the protein product and save as a cross-reference
 179      */
 180     new EnsemblCdna().addProteinProduct(transcript);
 181
 182     return transcript;
 183   }
 184
 185   /**
 186    * Returns a list of the transcript features on the sequence whose Parent is
 187    * the gene for the accession id.
 188    *
 189    * @param accId
 190    * @param geneSequence
 191    * @return
 192    */
 193   protected List<SequenceFeature> getTranscriptFeatures(String accId,
 194           SequenceI geneSequence)
 195   {
 196     List<SequenceFeature> transcriptFeatures = new ArrayList<SequenceFeature>();
 197
 198     String parentIdentifier = "gene:" + accId;
 199     SequenceFeature[] sfs = geneSequence.getSequenceFeatures();
 200
 201     if (sfs != null)
 202     {
 203       for (SequenceFeature sf : sfs)
 204       {
 205         if (isTranscript(sf.getType()))
 206         {
 207           String parent = (String) sf.getValue(PARENT);
 208           if (parentIdentifier.equals(parent))
 209           {
 210             transcriptFeatures.add(sf);
 211           }
 212         }
 213       }
 214     }
 215
 216     return transcriptFeatures;
 217   }
 218
 219   @Override
 220   public String getDescription()
 221   {
 222     return "Fetches all transcripts and variant features for a gene";
 223   }
 224
 225   /**
 226    * Default test query is a transcript
 227    */
 228   @Override
 229   public String getTestQuery()
 230   {
 231     return "ENSG00000157764"; // BRAF, 5 transcripts, reverse strand
 232     // ENSG00000090266 // NDUFB2, 15 transcripts, forward strand
 233     // ENSG00000101812 // H2BFM histone, 3 transcripts, forward strand
 234     // ENSG00000123569 // H2BFWT histone, 2 transcripts, reverse strand
 235   }
 236
 237   /**
 238    * Answers true for a feature of type 'gene' (or a sub-type of gene in the
 239    * Sequence Ontology), whose ID is the accession we are retrieving
 240    */
 241   @Override
 242   protected boolean identifiesSequence(SequenceFeature sf, String accId)
 243   {
 244     if (SequenceOntologyFactory.getInstance().isA(sf.getType(),
 245             SequenceOntologyI.GENE))
 246     {
 247       String id = (String) sf.getValue(ID);
 248       if (("gene:" + accId).equals(id))
 249       {
 250         return true;
 251       }
 252     }
 253     return false;
 254   }
 255
 256   /**
 257    * Answers true unless feature type is 'gene', or 'transcript' with a parent
 258    * which is a different gene. We need the gene features to identify the range,
 259    * but it is redundant information on the gene sequence. Checking the parent
 260    * allows us to drop transcript features which belong to different
 261    * (overlapping) genes.
 262    */
 263   @Override
 264   protected boolean retainFeature(SequenceFeature sf, String accessionId)
 265   {
 266     if (SequenceOntologyFactory.getInstance().isA(sf.getType(),
 267             SequenceOntologyI.GENE))
 268     {
 269       return false;
 270     }
 271
 272     if (isTranscript(sf.getType()))
 273     {
 274       String parent = (String) sf.getValue(PARENT);
 275       if (!("gene:" + accessionId).equals(parent))
 276       {
 277         return false;
 278       }
 279     }
 280     return true;
 281   }
 282
 283   /**
 284    * Answers false. This allows an optimisation - a single 'gene' feature is all
 285    * that is needed to identify the positions of the gene on the genomic
 286    * sequence.
 287    */
 288   @Override
 289   protected boolean isSpliceable()
 290   {
 291     return false;
 292   }
 293
 294   @Override
 295   protected List<String> getCrossReferenceDatabases()
 296   {
 297     // found these for ENSG00000157764 on 30/01/2016:
 298     // return new String[] {"Vega_gene", "OTTG", "ENS_LRG_gene", "ArrayExpress",
 299     // "EntrezGene", "HGNC", "MIM_GENE", "MIM_MORBID", "WikiGene"};
 300     return super.getCrossReferenceDatabases();
 301   }
 302
 303   /**
 304    * Override to do nothing as Ensembl doesn't return a protein sequence for a
 305    * gene identifier
 306    */
 307   @Override
 308   protected void addProteinProduct(SequenceI querySeq)
 309   {
 310   }
 311
 312 }