1 package jalview.ext.ensembl;
3 import jalview.datamodel.Alignment;
4 import jalview.datamodel.AlignmentI;
5 import jalview.datamodel.DBRefEntry;
6 import jalview.datamodel.DBRefSource;
7 import jalview.datamodel.Mapping;
8 import jalview.datamodel.SequenceFeature;
9 import jalview.datamodel.SequenceI;
10 import jalview.exceptions.JalviewException;
11 import jalview.io.FastaFile;
12 import jalview.io.FileParse;
13 import jalview.io.gff.SequenceOntology;
14 import jalview.util.DBRefUtils;
15 import jalview.util.MapList;
17 import java.io.IOException;
18 import java.net.MalformedURLException;
20 import java.util.ArrayList;
21 import java.util.Arrays;
22 import java.util.Collections;
23 import java.util.Comparator;
24 import java.util.List;
27 * Base class for Ensembl sequence fetchers
31 public abstract class EnsemblSeqProxy extends EnsemblRestClient
33 protected static final String PARENT = "Parent";
35 protected static final String ID = "ID";
37 public enum EnsemblSeqType
40 * type=genomic for the full dna including introns
45 * type=cdna for transcribed dna including UTRs
50 * type=cds for coding dna excluding UTRs
55 * type=protein for the peptide product sequence
60 * the value of the 'type' parameter to fetch this version of
65 EnsemblSeqType(String t)
70 public String getType()
78 * A comparator to sort ranges into ascending start position order
80 private class RangeSorter implements Comparator<int[]>
84 RangeSorter(boolean forward)
90 public int compare(int[] o1, int[] o2)
92 return (forwards ? 1 : -1) * Integer.compare(o1[0], o2[0]);
100 public EnsemblSeqProxy()
105 * Makes the sequence queries to Ensembl's REST service and returns an
106 * alignment consisting of the returned sequences
109 public AlignmentI getSequenceRecords(String query) throws Exception
111 long now = System.currentTimeMillis();
112 // TODO use a String... query vararg instead?
114 // danger: accession separator used as a regex here, a string elsewhere
115 // in this case it is ok (it is just a space), but (e.g.) '\' would not be
116 List<String> allIds = Arrays.asList(query.split(getAccessionSeparator()));
117 AlignmentI alignment = null;
121 * execute queries, if necessary in batches of the
122 * maximum allowed number of ids
124 int maxQueryCount = getMaximumQueryCount();
125 for (int v = 0, vSize = allIds.size(); v < vSize; v += maxQueryCount)
127 int p = Math.min(vSize, v + maxQueryCount);
128 List<String> ids = allIds.subList(v, p);
131 alignment = fetchSequences(ids, alignment);
132 } catch (Throwable r)
135 String msg = "Aborting ID retrieval after " + v
136 + " chunks. Unexpected problem (" + r.getLocalizedMessage()
138 System.err.println(msg);
139 if (alignment != null)
141 break; // return what we got
145 throw new JalviewException(msg, r);
151 * fetch and transfer genomic sequence features
153 for (String accId : allIds)
155 addFeaturesAndProduct(accId, alignment);
159 System.out.println(getClass().getName() + " took "
160 + (System.currentTimeMillis() - now) + "ms to fetch");
165 * Fetches Ensembl features using the /overlap REST endpoint, and adds them to
166 * the sequence in the alignment. Also fetches the protein product, maps it
167 * from the CDS features of the sequence, and saves it as a cross-reference of
173 protected void addFeaturesAndProduct(String accId, AlignmentI alignment)
178 * get 'dummy' genomic sequence with exon, cds and variation features
180 EnsemblOverlap gffFetcher = new EnsemblOverlap();
181 EnsemblFeatureType[] features = getFeaturesToFetch();
182 AlignmentI geneFeatures = gffFetcher.getSequenceRecords(accId,
184 if (geneFeatures.getHeight() > 0)
187 * transfer features to the query sequence
189 SequenceI genomicSequence = geneFeatures.getSequenceAt(0);
190 SequenceI querySeq = alignment.findName(accId);
191 transferFeatures(accId, genomicSequence, querySeq);
194 * fetch and map protein product, and add it as a cross-reference
195 * of the retrieved sequence
197 addProteinProduct(querySeq);
199 } catch (IOException e)
201 System.err.println("Error transferring Ensembl features: "
207 * Returns those sequence feature types to fetch from Ensembl. We may want
208 * features either because they are of interest to the user, or as means to
209 * identify the locations of the sequence on the genomic sequence (CDS
210 * features identify CDS, exon features identify cDNA etc).
214 protected abstract EnsemblFeatureType[] getFeaturesToFetch();
217 * Fetches and maps the protein product, and adds it as a cross-reference of
218 * the retrieved sequence
220 protected void addProteinProduct(SequenceI querySeq)
222 String accId = querySeq.getName();
225 AlignmentI protein = new EnsemblProtein().getSequenceRecords(accId);
226 if (protein == null || protein.getHeight() == 0)
228 System.out.println("Failed to retrieve protein for " + accId);
231 SequenceI proteinSeq = protein.getSequenceAt(0);
234 * need dataset sequences (to be the subject of mappings)
236 proteinSeq.createDatasetSequence();
237 querySeq.createDatasetSequence();
239 MapList mapList = mapCdsToProtein(querySeq, proteinSeq);
242 Mapping map = new Mapping(proteinSeq.getDatasetSequence(), mapList);
243 DBRefEntry dbr = new DBRefEntry(getDbSource(), getDbVersion(),
245 querySeq.getDatasetSequence().addDBRef(dbr);
247 } catch (Exception e)
250 .println(String.format("Error retrieving protein for %s: %s",
251 accId, e.getMessage()));
256 * Returns a mapping from dna to protein by inspecting sequence features of
257 * type "CDS" on the dna.
263 protected MapList mapCdsToProtein(SequenceI dnaSeq, SequenceI proteinSeq)
265 SequenceFeature[] sfs = dnaSeq.getSequenceFeatures();
271 List<int[]> ranges = new ArrayList<int[]>(50);
272 SequenceOntology so = SequenceOntology.getInstance();
274 int mappedDnaLength = 0;
277 * Map CDS columns of dna to peptide. No need to worry about reverse strand
278 * dna here since the retrieved sequence is as transcribed (reverse
279 * complement for reverse strand), i.e in the same sense as the peptide.
281 for (SequenceFeature sf : sfs)
284 * process a CDS feature (or a sub-type of CDS)
286 if (so.isA(sf.getType(), SequenceOntology.CDS))
288 ranges.add(new int[] { sf.getBegin(), sf.getEnd() });
289 mappedDnaLength += Math.abs(sf.getEnd() - sf.getBegin()) + 1;
292 int proteinLength = proteinSeq.getLength();
293 List<int[]> proteinRange = new ArrayList<int[]>();
294 proteinRange.add(new int[] { 1, proteinLength });
297 * dna length should map to protein (or protein minus stop codon)
299 if (mappedDnaLength == 3 * proteinLength
300 || mappedDnaLength == 3 * (proteinLength + 1))
302 return new MapList(ranges, proteinRange, 3, 1);
308 * Fetches sequences for the list of accession ids and adds them to the
309 * alignment. Returns the extended (or created) alignment.
314 * @throws JalviewException
315 * @throws IOException
317 protected AlignmentI fetchSequences(List<String> ids, AlignmentI alignment)
318 throws JalviewException, IOException
320 if (!isEnsemblAvailable())
323 throw new JalviewException("ENSEMBL Rest API not available.");
325 FileParse fp = getSequenceReader(ids);
326 FastaFile fr = new FastaFile(fp);
327 if (fr.hasWarningMessage())
329 System.out.println(String.format(
330 "Warning when retrieving %d ids %s\n%s", ids.size(),
331 ids.toString(), fr.getWarningMessage()));
333 else if (fr.getSeqs().size() != ids.size())
335 System.out.println(String.format(
336 "Only retrieved %d sequences for %d query strings", fr
337 .getSeqs().size(), ids.size()));
339 if (fr.getSeqs().size() > 0)
341 AlignmentI seqal = new Alignment(
342 fr.getSeqsAsArray());
343 for (SequenceI sq:seqal.getSequences())
345 if (sq.getDescription() == null)
347 sq.setDescription(getDbName());
349 String name = sq.getName();
350 if (ids.contains(name)
351 || ids.contains(name.replace("ENSP", "ENST")))
353 DBRefUtils.parseToDbRef(sq, DBRefSource.ENSEMBL, "0", name);
356 if (alignment == null)
362 alignment.append(seqal);
369 * Returns the URL for the REST call
372 * @throws MalformedURLException
375 protected URL getUrl(List<String> ids) throws MalformedURLException
377 // ids are not used - they go in the POST body instead
378 StringBuffer urlstring = new StringBuffer(128);
379 urlstring.append(SEQUENCE_ID_URL);
381 // @see https://github.com/Ensembl/ensembl-rest/wiki/Output-formats
382 urlstring.append("?type=").append(getSourceEnsemblType().getType());
383 urlstring.append(("&Accept=text/x-fasta"));
385 URL url = new URL(urlstring.toString());
390 * A sequence/id POST request currently allows up to 50 queries
392 * @see http://rest.ensembl.org/documentation/info/sequence_id_post
395 public int getMaximumQueryCount()
401 public boolean useGetRequest()
407 public String getRequestMimeType()
409 return "application/json";
413 public String getResponseMimeType()
415 return "text/x-fasta";
420 * @return the configured sequence return type for this source
422 protected abstract EnsemblSeqType getSourceEnsemblType();
425 * Returns a list of [start, end] genomic ranges corresponding to the sequence
428 * The correspondence between the frames of reference is made by locating
429 * those features on the genomic sequence which identify the retrieved
430 * sequence. Specifically
432 * <li>genomic sequence is identified by "transcript" features with
433 * ID=transcript:transcriptId</li>
434 * <li>cdna sequence is identified by "exon" features with
435 * Parent=transcript:transcriptId</li>
436 * <li>cds sequence is identified by "CDS" features with
437 * Parent=transcript:transcriptId</li>
440 * The returned ranges are sorted to run forwards (for positive strand) or
441 * backwards (for negative strand). Aborts and returns null if both positive
442 * and negative strand are found (this should not normally happen).
448 protected MapList getGenomicRanges(SequenceFeature[] sfs, String accId)
451 * generously size for initial number of cds regions
452 * (worst case titin Q8WZ42 has c. 313 exons)
454 List<int[]> regions = new ArrayList<int[]>(100);
455 int mappedLength = 0;
456 int direction = 1; // forward
457 boolean directionSet = false;
459 for (SequenceFeature sf : sfs)
462 * accept the target feature type or a specialisation of it
463 * (e.g. coding_exon for exon)
465 if (identifiesSequence(sf, accId))
467 int strand = sf.getStrand();
469 if (directionSet && strand != direction)
471 // abort - mix of forward and backward
472 System.err.println("Error: forward and backward strand for "
480 * add to CDS ranges, semi-sorted forwards/backwards
484 regions.add(0, new int[] { sf.getEnd(), sf.getBegin() });
488 regions.add(new int[] { sf.getBegin(), sf.getEnd() });
490 mappedLength += Math.abs(sf.getEnd() - sf.getBegin() + 1);
495 * a final sort is needed since Ensembl returns CDS sorted within source
496 * (havana / ensembl_havana)
498 Collections.sort(regions, new RangeSorter(direction == 1));
500 List<int[]> to = new ArrayList<int[]>();
501 to.add(new int[] { 1, mappedLength });
503 return new MapList(regions, to, 1, 1);
507 * Returns true if the sequence feature identifies positions of the genomic
508 * sequence feature which are within the sequence being retrieved.
514 protected abstract boolean identifiesSequence(SequenceFeature sf,
518 * Transfers the sequence feature to the target sequence, adjusting its start
519 * and end range based on the 'overlap' ranges. Features which do not overlap
520 * the target sequence are ignored, as are features with a parent other than
521 * the target sequence id.
524 * @param targetSequence
527 protected void transferFeature(SequenceFeature sf,
528 SequenceI targetSequence, MapList overlap)
530 String parent = (String) sf.getValue(PARENT);
531 if (parent != null && !parent.contains(targetSequence.getName()))
533 // this genomic feature belongs to a different transcript
537 int start = sf.getBegin();
538 int end = sf.getEnd();
539 int[] mappedRange = overlap.locateInTo(start, end);
541 if (mappedRange != null)
543 SequenceFeature copy = new SequenceFeature(sf);
544 int offset = targetSequence.getStart() - 1;
545 copy.setBegin(offset + Math.min(mappedRange[0], mappedRange[1]));
546 copy.setEnd(offset + Math.max(mappedRange[0], mappedRange[1]));
547 targetSequence.addSequenceFeature(copy);
550 * for sequence_variant, make an additional feature with consequence
552 if (SequenceOntology.getInstance().isSequenceVariant(sf.getType()))
554 String consequence = (String) sf.getValue("consequence_type");
555 if (consequence != null)
557 SequenceFeature sf2 = new SequenceFeature("consequence",
558 consequence, copy.getBegin(), copy.getEnd(), 0f,
560 targetSequence.addSequenceFeature(sf2);
568 * Transfers features from sourceSequence to targetSequence
571 * @param sourceSequence
572 * @param targetSequence
574 protected void transferFeatures(String accessionId,
575 SequenceI sourceSequence, SequenceI targetSequence)
577 if (sourceSequence == null || targetSequence == null)
582 SequenceFeature[] sfs = sourceSequence.getSequenceFeatures();
583 MapList overlap = getGenomicRanges(sfs, accessionId);
585 final boolean forwardStrand = overlap.isFromForwardStrand();
588 * sort features by start position (descending if reverse strand)
589 * before transferring (in forwards order) to the target sequence
591 Arrays.sort(sfs, new Comparator<SequenceFeature>()
594 public int compare(SequenceFeature o1, SequenceFeature o2)
596 int c = Integer.compare(o1.getBegin(), o2.getBegin());
597 return forwardStrand ? c : -c;
601 for (SequenceFeature sf : sfs)
603 if (retainFeature(sf.getType()))
605 transferFeature(sf, targetSequence, overlap);
611 * Answers true if the feature type is one to attach to the retrieved sequence
616 protected boolean retainFeature(@SuppressWarnings("unused") String type)
618 return true; // default is to keep all