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