src/jalview/ext/ensembl/EnsemblGene.java

   1 package jalview.ext.ensembl;
   2
   3 import jalview.datamodel.AlignedCodonFrame;
   4 import jalview.datamodel.AlignmentI;
   5 import jalview.datamodel.SequenceFeature;
   6 import jalview.datamodel.SequenceI;
   7 import jalview.io.gff.SequenceOntology;
   8 import jalview.util.MapList;
   9
  10 import java.util.ArrayList;
  11 import java.util.List;
  12
  13 /**
  14  * A class that fetches genomic sequence and all transcripts for an Ensembl gene
  15  *
  16  * @author gmcarstairs
  17  */
  18 public class EnsemblGene extends EnsemblSeqProxy
  19 {
  20   private static final EnsemblFeatureType[] FEATURES_TO_FETCH = {
  21       EnsemblFeatureType.gene, EnsemblFeatureType.transcript,
  22       EnsemblFeatureType.exon, EnsemblFeatureType.cds,
  23       EnsemblFeatureType.variation };
  24
  25   @Override
  26   public String getDbName()
  27   {
  28     return "ENSEMBL (GENE)";
  29   }
  30
  31   @Override
  32   protected EnsemblFeatureType[] getFeaturesToFetch()
  33   {
  34     return FEATURES_TO_FETCH;
  35   }
  36
  37   @Override
  38   protected EnsemblSeqType getSourceEnsemblType()
  39   {
  40     return EnsemblSeqType.GENOMIC;
  41   }
  42
  43   /**
  44    * Builds an alignment of all transcripts for the requested gene:
  45    * <ul>
  46    * <li>fetches the gene sequence</li>
  47    * <li>fetches features on the sequence</li>
  48    * <li>identifies "transcript" features whose Parent is the requested gene</li>
  49    * <li>fetches the transcript sequence for each transcript</li>
  50    * <li>makes a mapping from the gene to each transcript</li>
  51    * <li>copies features from gene to transcript sequences</li>
  52    * <li>fetches the protein sequence for each transcript, maps and saves it as
  53    * a cross-reference</li>
  54    * <li>aligns each transcript against the gene sequence based on the position
  55    * mappings</li>
  56    * </ul>
  57    */
  58   @Override
  59   public AlignmentI getSequenceRecords(String query) throws Exception
  60   {
  61     AlignmentI al = super.getSequenceRecords(query);
  62     if (al.getHeight() > 0)
  63     {
  64       getTranscripts(al, query);
  65     }
  66
  67     return al;
  68   }
  69
  70   /**
  71    * Find and fetch all transcripts for the gene, as identified by "transcript"
  72    * features whose Parent is the requested gene
  73    *
  74    * @param al
  75    * @param accId
  76    * @throws Exception
  77    */
  78   protected void getTranscripts(AlignmentI al, String accId)
  79           throws Exception
  80   {
  81     SequenceI gene = al.getSequenceAt(0);
  82     List<String> transcriptIds = getTranscriptIds(accId, gene);
  83
  84     // TODO: could just use features and genomic sequence
  85     // to generate the transcript sequences - faster
  86     // could also grab "Name" as transcript description (gene name)
  87     for (String transcriptId : transcriptIds)
  88     {
  89       /*
  90        * fetch and map the transcript sequence; we can pass in the gene
  91        * sequence with features marked to save fetching it again
  92        */
  93       EnsemblCdna cdnaFetcher = new EnsemblCdna();
  94       AlignmentI al2 = cdnaFetcher.getSequenceRecords(transcriptId,
  95               gene);
  96       for (SequenceI seq : al2.getSequences())
  97       {
  98         /*
  99          * build mapping from gene sequence to transcript
 100          */
 101         MapList mapping = cdnaFetcher.getGenomicRanges(gene, transcriptId,
 102                 seq.getStart());
 103
 104         /*
 105          * align the transcript to the gene
 106          */
 107         AlignedCodonFrame acf = new AlignedCodonFrame();
 108         acf.addMap(gene, seq, mapping);
 109         char gap = al.getGapCharacter();
 110         // AlignmentUtils.alignSequenceAs(seq, gene, acf, String.valueOf(gap),
 111         // gap, false, false);
 112
 113         al.addSequence(seq);
 114       }
 115     }
 116   }
 117
 118   /**
 119    * Returns a list of the ids of transcript features on the sequence whose
 120    * Parent is the gene for the accession id
 121    *
 122    * @param accId
 123    * @param geneSequence
 124    * @return
 125    */
 126   protected List<String> getTranscriptIds(String accId, SequenceI geneSequence)
 127   {
 128     SequenceOntology so = SequenceOntology.getInstance();
 129     List<String> transcriptIds = new ArrayList<String>();
 130
 131     /*
 132      * scan for transcript features belonging to our gene;
 133      * also remove any which belong to other genes
 134      */
 135     SequenceFeature[] sfs = geneSequence.getSequenceFeatures();
 136     List<SequenceFeature> keptFeatures = new ArrayList<SequenceFeature>();
 137     boolean featureDropped = false;
 138     String parentIdentifier = "gene:" + accId;
 139     for (SequenceFeature sf : sfs)
 140     {
 141       if (so.isA(sf.getType(), SequenceOntology.TRANSCRIPT))
 142       {
 143         String parent = (String) sf.getValue(PARENT);
 144         if (parentIdentifier.equals(parent))
 145         {
 146           transcriptIds.add((String) sf.getValue("transcript_id"));
 147           keptFeatures.add(sf);
 148         }
 149         else
 150         {
 151           featureDropped = true;
 152         }
 153       }
 154       else
 155       {
 156         keptFeatures.add(sf);
 157       }
 158     }
 159     if (featureDropped)
 160     {
 161       geneSequence.getDatasetSequence().setSequenceFeatures(
 162               keptFeatures.toArray(new SequenceFeature[keptFeatures
 163                       .size()]));
 164     }
 165     return transcriptIds;
 166   }
 167
 168   @Override
 169   public String getDescription()
 170   {
 171     return "Fetches all transcripts and variant features for a gene";
 172   }
 173
 174   /**
 175    * Default test query is a transcript
 176    */
 177   @Override
 178   public String getTestQuery()
 179   {
 180     return "ENSG00000157764"; // reverse strand
 181     // ENSG00000090266 // forward strand
 182   }
 183
 184   /**
 185    * Answers true for a feature of type 'gene' (or a sub-type of gene in the
 186    * Sequence Ontology), whose ID is the accession we are retrieving
 187    */
 188   @Override
 189   protected boolean identifiesSequence(SequenceFeature sf, String accId)
 190   {
 191     if (SequenceOntology.getInstance().isA(sf.getType(),
 192             SequenceOntology.GENE))
 193     {
 194       String id = (String) sf.getValue(ID);
 195       if (("gene:" + accId).equals(id))
 196       {
 197         return true;
 198       }
 199     }
 200     return false;
 201   }
 202
 203   /**
 204    * Answers true unless feature type is 'gene'. We need the gene features to
 205    * identify the range, but it is redundant information on the gene sequence.
 206    */
 207   @Override
 208   protected boolean retainFeature(SequenceFeature sf, String accessionId)
 209   {
 210     return !SequenceOntology.getInstance().isA(sf.getType(),
 211             SequenceOntology.GENE);
 212   }
 213
 214 }