1 package jalview.ext.ensembl;
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.SequenceOntology;
8 import jalview.util.MapList;
10 import java.util.ArrayList;
11 import java.util.Arrays;
12 import java.util.List;
15 * A class that fetches genomic sequence and all transcripts for an Ensembl gene
19 public class EnsemblGene extends EnsemblSeqProxy
21 private static final EnsemblFeatureType[] FEATURES_TO_FETCH = {
22 EnsemblFeatureType.gene, EnsemblFeatureType.transcript,
23 EnsemblFeatureType.exon, EnsemblFeatureType.cds,
24 EnsemblFeatureType.variation };
27 public String getDbName()
29 return "ENSEMBL (GENE)";
33 protected EnsemblFeatureType[] getFeaturesToFetch()
35 return FEATURES_TO_FETCH;
39 protected EnsemblSeqType getSourceEnsemblType()
41 return EnsemblSeqType.GENOMIC;
45 * Builds an alignment of all transcripts for the requested gene:
47 * <li>fetches the gene sequence</li>
48 * <li>fetches features on the sequence</li>
49 * <li>identifies "transcript" features whose Parent is the requested gene</li>
50 * <li>fetches the transcript sequence for each transcript</li>
51 * <li>makes a mapping from the gene to each transcript</li>
52 * <li>copies features from gene to transcript sequences</li>
53 * <li>fetches the protein sequence for each transcript, maps and saves it as
54 * a cross-reference</li>
55 * <li>aligns each transcript against the gene sequence based on the position
60 public AlignmentI getSequenceRecords(String query) throws Exception
62 // TODO ? if an ENST identifier is supplied, convert to ENSG?
63 AlignmentI al = super.getSequenceRecords(query);
64 if (al.getHeight() > 0)
66 getTranscripts(al, query);
73 * Constructs all transcripts for the gene, as identified by "transcript"
74 * features whose Parent is the requested gene. The coding transcript
75 * sequences (i.e. with introns omitted) are added to the alignment.
81 protected void getTranscripts(AlignmentI al, String accId)
84 SequenceI gene = al.getSequenceAt(0);
85 List<SequenceFeature> transcriptFeatures = getTranscriptFeatures(accId,
88 for (SequenceFeature transcriptFeature : transcriptFeatures)
90 makeTranscript(transcriptFeature, al, gene);
95 * Constructs a spliced transcript sequence by finding 'exon' features for the
96 * given id (or failing that 'CDS'). Copies features on to the new sequence.
97 * 'Aligns' the new sequence against the gene sequence by padding with gaps,
98 * and adds it to the alignment.
100 * @param transcriptFeature
102 * the alignment to which to add the new sequence
104 * the parent gene sequence, with features
107 SequenceI makeTranscript(SequenceFeature transcriptFeature,
108 AlignmentI al, SequenceI gene)
110 String accId = (String) transcriptFeature.getValue("transcript_id");
117 * NB we are mapping from gene sequence (not genome), so do not
118 * need to check for reverse strand (gene and transcript sequences
119 * are in forward sense)
123 * make a gene-length sequence filled with gaps
124 * we will fill in the bases for transcript regions
126 char[] seqChars = new char[gene.getLength()];
127 Arrays.fill(seqChars, al.getGapCharacter());
130 * look for exon features of the transcript, failing that for CDS
131 * (for example ENSG00000124610 has 1 CDS but no exon features)
133 String parentId = "transcript:" + accId;
134 List<SequenceFeature> splices = findFeatures(gene,
135 SequenceOntology.EXON, parentId);
136 if (splices.isEmpty())
138 splices = findFeatures(gene, SequenceOntology.CDS, parentId);
141 int transcriptLength = 0;
142 final char[] geneChars = gene.getSequence();
143 int offset = gene.getStart(); // to convert to 0-based positions
144 List<int[]> mappedFrom = new ArrayList<int[]>();
146 for (SequenceFeature sf : splices)
148 int start = sf.getBegin() - offset;
149 int end = sf.getEnd() - offset;
150 int spliceLength = end - start + 1;
151 System.arraycopy(geneChars, start, seqChars, start, spliceLength);
152 transcriptLength += spliceLength;
153 mappedFrom.add(new int[] { sf.getBegin(), sf.getEnd() });
156 Sequence transcript = new Sequence(accId, seqChars, 1, transcriptLength);
157 String geneName = (String) transcriptFeature.getValue(NAME);
158 if (geneName != null)
160 transcript.setDescription(geneName);
162 transcript.createDatasetSequence();
164 al.addSequence(transcript);
167 * transfer features to the new sequence; we use EnsemblCdna to do this,
168 * to filter out unwanted features types (see method retainFeature)
170 List<int[]> mapTo = new ArrayList<int[]>();
171 mapTo.add(new int[] { 1, transcriptLength });
172 MapList mapping = new MapList(mappedFrom, mapTo, 1, 1);
173 new EnsemblCdna().transferFeatures(gene.getSequenceFeatures(),
174 transcript.getDatasetSequence(), mapping, parentId);
177 * and finally fetch the protein product and save as a cross-reference
179 addProteinProduct(transcript);
185 * Returns a list of the transcript features on the sequence whose Parent is
186 * the gene for the accession id.
189 * @param geneSequence
192 protected List<SequenceFeature> getTranscriptFeatures(String accId,
193 SequenceI geneSequence)
195 List<SequenceFeature> transcriptFeatures = new ArrayList<SequenceFeature>();
197 String parentIdentifier = "gene:" + accId;
198 SequenceFeature[] sfs = geneSequence.getSequenceFeatures();
202 for (SequenceFeature sf : sfs)
204 if (isTranscript(sf.getType()))
206 String parent = (String) sf.getValue(PARENT);
207 if (parentIdentifier.equals(parent))
209 transcriptFeatures.add(sf);
215 return transcriptFeatures;
219 public String getDescription()
221 return "Fetches all transcripts and variant features for a gene";
225 * Default test query is a transcript
228 public String getTestQuery()
230 return "ENSG00000157764"; // BRAF, 5 transcripts, reverse strand
231 // ENSG00000090266 // NDUFB2, 15 transcripts, forward strand
232 // ENSG00000101812 // H2BFM histone, 3 transcripts, forward strand
233 // ENSG00000123569 // H2BFWT histone, 2 transcripts, reverse strand
237 * Answers true for a feature of type 'gene' (or a sub-type of gene in the
238 * Sequence Ontology), whose ID is the accession we are retrieving
241 protected boolean identifiesSequence(SequenceFeature sf, String accId)
243 if (SequenceOntology.getInstance().isA(sf.getType(),
244 SequenceOntology.GENE))
246 String id = (String) sf.getValue(ID);
247 if (("gene:" + accId).equals(id))
256 * Answers true unless feature type is 'gene', or 'transcript' with a parent
257 * which is a different gene. We need the gene features to identify the range,
258 * but it is redundant information on the gene sequence. Checking the parent
259 * allows us to drop transcript features which belong to different
260 * (overlapping) genes.
263 protected boolean retainFeature(SequenceFeature sf, String accessionId)
265 if (SequenceOntology.getInstance().isA(sf.getType(),
266 SequenceOntology.GENE))
271 if (isTranscript(sf.getType()))
273 String parent = (String) sf.getValue(PARENT);
274 if (!("gene:" + accessionId).equals(parent))
283 * Answers false. This allows an optimisation - a single 'gene' feature is all
284 * that is needed to identify the positions of the gene on the genomic
288 protected boolean isSpliceable()