src/jalview/ext/ensembl/EnsemblSeqProxy.java

   1 package jalview.ext.ensembl;
   2
   3 import jalview.analysis.AlignmentUtils;
   4 import jalview.analysis.Dna;
   5 import jalview.datamodel.Alignment;
   6 import jalview.datamodel.AlignmentI;
   7 import jalview.datamodel.DBRefEntry;
   8 import jalview.datamodel.DBRefSource;
   9 import jalview.datamodel.Mapping;
  10 import jalview.datamodel.SequenceFeature;
  11 import jalview.datamodel.SequenceI;
  12 import jalview.exceptions.JalviewException;
  13 import jalview.io.FastaFile;
  14 import jalview.io.FileParse;
  15 import jalview.io.gff.SequenceOntologyFactory;
  16 import jalview.io.gff.SequenceOntologyI;
  17 import jalview.util.Comparison;
  18 import jalview.util.DBRefUtils;
  19 import jalview.util.MapList;
  20
  21 import java.io.IOException;
  22 import java.net.MalformedURLException;
  23 import java.net.URL;
  24 import java.util.ArrayList;
  25 import java.util.Arrays;
  26 import java.util.Collections;
  27 import java.util.Comparator;
  28 import java.util.List;
  29
  30 /**
  31  * Base class for Ensembl sequence fetchers
  32  *
  33  * @see http://rest.ensembl.org/documentation/info/sequence_id
  34  * @author gmcarstairs
  35  */
  36 public abstract class EnsemblSeqProxy extends EnsemblRestClient
  37 {
  38   private static final String ALLELES = "alleles";
  39
  40   protected static final String PARENT = "Parent";
  41
  42   protected static final String ID = "ID";
  43
  44   protected static final String NAME = "Name";
  45
  46   protected static final String DESCRIPTION = "description";
  47
  48   /*
  49    * enum for 'type' parameter to the /sequence REST service
  50    */
  51   public enum EnsemblSeqType
  52   {
  53     /**
  54      * type=genomic to fetch full dna including introns
  55      */
  56     GENOMIC("genomic"),
  57
  58     /**
  59      * type=cdna to fetch coding dna including UTRs
  60      */
  61     CDNA("cdna"),
  62
  63     /**
  64      * type=cds to fetch coding dna excluding UTRs
  65      */
  66     CDS("cds"),
  67
  68     /**
  69      * type=protein to fetch peptide product sequence
  70      */
  71     PROTEIN("protein");
  72
  73     /*
  74      * the value of the 'type' parameter to fetch this version of
  75      * an Ensembl sequence
  76      */
  77     private String type;
  78
  79     EnsemblSeqType(String t)
  80     {
  81       type = t;
  82     }
  83
  84     public String getType()
  85     {
  86       return type;
  87     }
  88
  89   }
  90
  91   /**
  92    * A comparator to sort ranges into ascending start position order
  93    */
  94   private class RangeSorter implements Comparator<int[]>
  95   {
  96     boolean forwards;
  97
  98     RangeSorter(boolean forward)
  99     {
 100       forwards = forward;
 101     }
 102
 103     @Override
 104     public int compare(int[] o1, int[] o2)
 105     {
 106       return (forwards ? 1 : -1) * Integer.compare(o1[0], o2[0]);
 107     }
 108
 109   }
 110
 111   /**
 112    * Default constructor (to use rest.ensembl.org)
 113    */
 114   public EnsemblSeqProxy()
 115   {
 116     super();
 117   }
 118
 119   /**
 120    * Constructor given the target domain to fetch data from
 121    */
 122   public EnsemblSeqProxy(String d)
 123   {
 124     super(d);
 125   }
 126
 127   /**
 128    * Makes the sequence queries to Ensembl's REST service and returns an
 129    * alignment consisting of the returned sequences.
 130    */
 131   @Override
 132   public AlignmentI getSequenceRecords(String query) throws Exception
 133   {
 134     // TODO use a String... query vararg instead?
 135
 136     // danger: accession separator used as a regex here, a string elsewhere
 137     // in this case it is ok (it is just a space), but (e.g.) '\' would not be
 138     List<String> allIds = Arrays.asList(query
 139             .split(getAccessionSeparator()));
 140     AlignmentI alignment = null;
 141     inProgress = true;
 142
 143     /*
 144      * execute queries, if necessary in batches of the
 145      * maximum allowed number of ids
 146      */
 147     int maxQueryCount = getMaximumQueryCount();
 148     for (int v = 0, vSize = allIds.size(); v < vSize; v += maxQueryCount)
 149     {
 150       int p = Math.min(vSize, v + maxQueryCount);
 151       List<String> ids = allIds.subList(v, p);
 152       try
 153       {
 154         alignment = fetchSequences(ids, alignment);
 155       } catch (Throwable r)
 156       {
 157         inProgress = false;
 158         String msg = "Aborting ID retrieval after " + v
 159                 + " chunks. Unexpected problem (" + r.getLocalizedMessage()
 160                 + ")";
 161         System.err.println(msg);
 162         break;
 163       }
 164     }
 165
 166     if (alignment == null)
 167     {
 168       return null;
 169     }
 170
 171     /*
 172      * fetch and transfer genomic sequence features,
 173      * fetch protein product and add as cross-reference
 174      */
 175     for (String accId : allIds)
 176     {
 177       addFeaturesAndProduct(accId, alignment);
 178     }
 179
 180     for (SequenceI seq : alignment.getSequences())
 181     {
 182       getCrossReferences(seq);
 183     }
 184
 185     return alignment;
 186   }
 187
 188   /**
 189    * Fetches Ensembl features using the /overlap REST endpoint, and adds them to
 190    * the sequence in the alignment. Also fetches the protein product, maps it
 191    * from the CDS features of the sequence, and saves it as a cross-reference of
 192    * the dna sequence.
 193    *
 194    * @param accId
 195    * @param alignment
 196    */
 197   protected void addFeaturesAndProduct(String accId, AlignmentI alignment)
 198   {
 199     if (alignment == null)
 200     {
 201       return;
 202     }
 203
 204     try
 205     {
 206       /*
 207        * get 'dummy' genomic sequence with exon, cds and variation features
 208        */
 209       SequenceI genomicSequence = null;
 210       EnsemblFeatures gffFetcher = new EnsemblFeatures(getDomain());
 211       EnsemblFeatureType[] features = getFeaturesToFetch();
 212       AlignmentI geneFeatures = gffFetcher.getSequenceRecords(accId,
 213               features);
 214       if (geneFeatures.getHeight() > 0)
 215       {
 216         genomicSequence = geneFeatures.getSequenceAt(0);
 217       }
 218       if (genomicSequence != null)
 219       {
 220         /*
 221          * transfer features to the query sequence
 222          */
 223         SequenceI querySeq = alignment.findName(accId);
 224         if (transferFeatures(accId, genomicSequence, querySeq))
 225         {
 226
 227           /*
 228            * fetch and map protein product, and add it as a cross-reference
 229            * of the retrieved sequence
 230            */
 231           addProteinProduct(querySeq);
 232         }
 233       }
 234     } catch (IOException e)
 235     {
 236       System.err.println("Error transferring Ensembl features: "
 237               + e.getMessage());
 238     }
 239   }
 240
 241   /**
 242    * Returns those sequence feature types to fetch from Ensembl. We may want
 243    * features either because they are of interest to the user, or as means to
 244    * identify the locations of the sequence on the genomic sequence (CDS
 245    * features identify CDS, exon features identify cDNA etc).
 246    *
 247    * @return
 248    */
 249   protected abstract EnsemblFeatureType[] getFeaturesToFetch();
 250
 251   /**
 252    * Fetches and maps the protein product, and adds it as a cross-reference of
 253    * the retrieved sequence
 254    */
 255   protected void addProteinProduct(SequenceI querySeq)
 256   {
 257     String accId = querySeq.getName();
 258     try
 259     {
 260       AlignmentI protein = new EnsemblProtein(getDomain())
 261               .getSequenceRecords(accId);
 262       if (protein == null || protein.getHeight() == 0)
 263       {
 264         System.out.println("No protein product found for " + accId);
 265         return;
 266       }
 267       SequenceI proteinSeq = protein.getSequenceAt(0);
 268
 269       /*
 270        * need dataset sequences (to be the subject of mappings)
 271        */
 272       proteinSeq.createDatasetSequence();
 273       querySeq.createDatasetSequence();
 274
 275       MapList mapList = AlignmentUtils.mapCdsToProtein(querySeq, proteinSeq);
 276       if (mapList != null)
 277       {
 278         // clunky: ensure Uniprot xref if we have one is on mapped sequence
 279         SequenceI ds = proteinSeq.getDatasetSequence();
 280         ds.setSourceDBRef(proteinSeq.getSourceDBRef());
 281
 282         Mapping map = new Mapping(ds, mapList);
 283         DBRefEntry dbr = new DBRefEntry(getDbSource(),
 284                 getEnsemblDataVersion(), proteinSeq.getName(), map);
 285         querySeq.getDatasetSequence().addDBRef(dbr);
 286
 287         /*
 288          * copy exon features to protein, compute peptide variants from dna
 289          * variants and add as features on the protein sequence ta-da
 290          */
 291         AlignmentUtils.computeProteinFeatures(querySeq, proteinSeq, mapList);
 292       }
 293     } catch (Exception e)
 294     {
 295       System.err
 296               .println(String.format("Error retrieving protein for %s: %s",
 297                       accId, e.getMessage()));
 298     }
 299   }
 300
 301   /**
 302    * Get database xrefs from Ensembl, and attach them to the sequence
 303    *
 304    * @param seq
 305    */
 306   protected void getCrossReferences(SequenceI seq)
 307   {
 308     while (seq.getDatasetSequence() != null)
 309     {
 310       seq = seq.getDatasetSequence();
 311     }
 312
 313     EnsemblXref xrefFetcher = new EnsemblXref(getDomain());
 314     List<DBRefEntry> xrefs = xrefFetcher.getCrossReferences(seq.getName());
 315     for (DBRefEntry xref : xrefs)
 316     {
 317       seq.addDBRef(xref);
 318       /*
 319        * Save any Uniprot xref to be the reference for SIFTS mapping
 320        */
 321       if (DBRefSource.UNIPROT.equals(xref.getSource()))
 322       {
 323         seq.setSourceDBRef(xref);
 324       }
 325     }
 326
 327     /*
 328      * and add a reference to itself
 329      */
 330     DBRefEntry self = new DBRefEntry(getDbSource(),
 331             getEnsemblDataVersion(), seq.getName());
 332     seq.addDBRef(self);
 333   }
 334
 335   /**
 336    * Fetches sequences for the list of accession ids and adds them to the
 337    * alignment. Returns the extended (or created) alignment.
 338    *
 339    * @param ids
 340    * @param alignment
 341    * @return
 342    * @throws JalviewException
 343    * @throws IOException
 344    */
 345   protected AlignmentI fetchSequences(List<String> ids, AlignmentI alignment)
 346           throws JalviewException, IOException
 347   {
 348     if (!isEnsemblAvailable())
 349     {
 350       inProgress = false;
 351       throw new JalviewException("ENSEMBL Rest API not available.");
 352     }
 353     FileParse fp = getSequenceReader(ids);
 354     FastaFile fr = new FastaFile(fp);
 355     if (fr.hasWarningMessage())
 356     {
 357       System.out.println(String.format(
 358               "Warning when retrieving %d ids %s\n%s", ids.size(),
 359               ids.toString(), fr.getWarningMessage()));
 360     }
 361     else if (fr.getSeqs().size() != ids.size())
 362     {
 363       System.out.println(String.format(
 364               "Only retrieved %d sequences for %d query strings", fr
 365                       .getSeqs().size(), ids.size()));
 366     }
 367
 368     if (fr.getSeqs().size() == 1 && fr.getSeqs().get(0).getLength() == 0)
 369     {
 370       /*
 371        * POST request has returned an empty FASTA file e.g. for invalid id
 372        */
 373       throw new IOException("No data returned for " + ids);
 374     }
 375
 376     if (fr.getSeqs().size() > 0)
 377     {
 378       AlignmentI seqal = new Alignment(
 379               fr.getSeqsAsArray());
 380       for (SequenceI sq:seqal.getSequences())
 381       {
 382         if (sq.getDescription() == null)
 383         {
 384           sq.setDescription(getDbName());
 385         }
 386         String name = sq.getName();
 387         if (ids.contains(name)
 388                 || ids.contains(name.replace("ENSP", "ENST")))
 389         {
 390           DBRefUtils.parseToDbRef(sq, getDbSource(),
 391                   getEnsemblDataVersion(), name);
 392         }
 393       }
 394       if (alignment == null)
 395       {
 396         alignment = seqal;
 397       }
 398       else
 399       {
 400         alignment.append(seqal);
 401       }
 402     }
 403     return alignment;
 404   }
 405
 406   /**
 407    * Returns the URL for the REST call
 408    *
 409    * @return
 410    * @throws MalformedURLException
 411    */
 412   @Override
 413   protected URL getUrl(List<String> ids) throws MalformedURLException
 414   {
 415     /*
 416      * a single id is included in the URL path
 417      * multiple ids go in the POST body instead
 418      */
 419     StringBuffer urlstring = new StringBuffer(128);
 420     urlstring.append(getDomain() + "/sequence/id");
 421     if (ids.size() == 1)
 422     {
 423       urlstring.append("/").append(ids.get(0));
 424     }
 425     // @see https://github.com/Ensembl/ensembl-rest/wiki/Output-formats
 426     urlstring.append("?type=").append(getSourceEnsemblType().getType());
 427     urlstring.append(("&Accept=text/x-fasta"));
 428
 429     URL url = new URL(urlstring.toString());
 430     return url;
 431   }
 432
 433   /**
 434    * A sequence/id POST request currently allows up to 50 queries
 435    *
 436    * @see http://rest.ensembl.org/documentation/info/sequence_id_post
 437    */
 438   @Override
 439   public int getMaximumQueryCount()
 440   {
 441     return 50;
 442   }
 443
 444   @Override
 445   protected boolean useGetRequest()
 446   {
 447     return false;
 448   }
 449
 450   @Override
 451   protected String getRequestMimeType(boolean multipleIds)
 452   {
 453     return multipleIds ? "application/json" : "text/x-fasta";
 454   }
 455
 456   @Override
 457   protected String getResponseMimeType()
 458   {
 459     return "text/x-fasta";
 460   }
 461
 462   /**
 463    *
 464    * @return the configured sequence return type for this source
 465    */
 466   protected abstract EnsemblSeqType getSourceEnsemblType();
 467
 468   /**
 469    * Returns a list of [start, end] genomic ranges corresponding to the sequence
 470    * being retrieved.
 471    *
 472    * The correspondence between the frames of reference is made by locating
 473    * those features on the genomic sequence which identify the retrieved
 474    * sequence. Specifically
 475    * <ul>
 476    * <li>genomic sequence is identified by "transcript" features with
 477    * ID=transcript:transcriptId</li>
 478    * <li>cdna sequence is identified by "exon" features with
 479    * Parent=transcript:transcriptId</li>
 480    * <li>cds sequence is identified by "CDS" features with
 481    * Parent=transcript:transcriptId</li>
 482    * </ul>
 483    *
 484    * The returned ranges are sorted to run forwards (for positive strand) or
 485    * backwards (for negative strand). Aborts and returns null if both positive
 486    * and negative strand are found (this should not normally happen).
 487    *
 488    * @param sourceSequence
 489    * @param accId
 490    * @param start
 491    *          the start position of the sequence we are mapping to
 492    * @return
 493    */
 494   protected MapList getGenomicRangesFromFeatures(SequenceI sourceSequence,
 495           String accId, int start)
 496   {
 497     SequenceFeature[] sfs = sourceSequence.getSequenceFeatures();
 498     if (sfs == null)
 499     {
 500       return null;
 501     }
 502
 503     /*
 504      * generously initial size for number of cds regions
 505      * (worst case titin Q8WZ42 has c. 313 exons)
 506      */
 507     List<int[]> regions = new ArrayList<int[]>(100);
 508     int mappedLength = 0;
 509     int direction = 1; // forward
 510     boolean directionSet = false;
 511
 512     for (SequenceFeature sf : sfs)
 513     {
 514       /*
 515        * accept the target feature type or a specialisation of it
 516        * (e.g. coding_exon for exon)
 517        */
 518       if (identifiesSequence(sf, accId))
 519       {
 520         int strand = sf.getStrand();
 521         strand = strand == 0 ? 1 : strand; // treat unknown as forward
 522
 523         if (directionSet && strand != direction)
 524         {
 525           // abort - mix of forward and backward
 526           System.err.println("Error: forward and backward strand for "
 527                   + accId);
 528             return null;
 529           }
 530           direction = strand;
 531           directionSet = true;
 532
 533           /*
 534            * add to CDS ranges, semi-sorted forwards/backwards
 535            */
 536           if (strand < 0)
 537           {
 538             regions.add(0, new int[] { sf.getEnd(), sf.getBegin() });
 539           }
 540           else
 541           {
 542           regions.add(new int[] { sf.getBegin(), sf.getEnd() });
 543         }
 544         mappedLength += Math.abs(sf.getEnd() - sf.getBegin() + 1);
 545
 546         if (!isSpliceable())
 547         {
 548           /*
 549            * 'gene' sequence is contiguous so we can stop as soon as its
 550            * identifying feature has been found
 551            */
 552           break;
 553         }
 554       }
 555     }
 556
 557     if (regions.isEmpty())
 558     {
 559       System.out.println("Failed to identify target sequence for " + accId
 560               + " from genomic features");
 561       return null;
 562     }
 563
 564     /*
 565      * a final sort is needed since Ensembl returns CDS sorted within source
 566      * (havana / ensembl_havana)
 567      */
 568     Collections.sort(regions, new RangeSorter(direction == 1));
 569
 570     List<int[]> to = Arrays.asList(new int[] { start,
 571         start + mappedLength - 1 });
 572
 573     return new MapList(regions, to, 1, 1);
 574   }
 575
 576   /**
 577    * Answers true if the sequence being retrieved may occupy discontiguous
 578    * regions on the genomic sequence.
 579    */
 580   protected boolean isSpliceable()
 581   {
 582     return true;
 583   }
 584
 585   /**
 586    * Returns true if the sequence feature marks positions of the genomic
 587    * sequence feature which are within the sequence being retrieved. For
 588    * example, an 'exon' feature whose parent is the target transcript marks the
 589    * cdna positions of the transcript.
 590    *
 591    * @param sf
 592    * @param accId
 593    * @return
 594    */
 595   protected abstract boolean identifiesSequence(SequenceFeature sf,
 596           String accId);
 597
 598   /**
 599    * Transfers the sequence feature to the target sequence, locating its start
 600    * and end range based on the mapping. Features which do not overlap the
 601    * target sequence are ignored.
 602    *
 603    * @param sf
 604    * @param targetSequence
 605    * @param mapping
 606    *          mapping from the sequence feature's coordinates to the target
 607    *          sequence
 608    * @param forwardStrand
 609    */
 610   protected void transferFeature(SequenceFeature sf,
 611           SequenceI targetSequence, MapList mapping, boolean forwardStrand)
 612   {
 613     int start = sf.getBegin();
 614     int end = sf.getEnd();
 615     int[] mappedRange = mapping.locateInTo(start, end);
 616
 617     if (mappedRange != null)
 618     {
 619       SequenceFeature copy = new SequenceFeature(sf);
 620       copy.setBegin(Math.min(mappedRange[0], mappedRange[1]));
 621       copy.setEnd(Math.max(mappedRange[0], mappedRange[1]));
 622       targetSequence.addSequenceFeature(copy);
 623
 624       /*
 625        * for sequence_variant on reverse strand, have to convert the allele
 626        * values to their complements
 627        */
 628       if (!forwardStrand
 629               && SequenceOntologyFactory.getInstance().isA(sf.getType(),
 630                       SequenceOntologyI.SEQUENCE_VARIANT))
 631       {
 632         reverseComplementAlleles(copy);
 633       }
 634     }
 635   }
 636
 637   /**
 638    * Change the 'alleles' value of a feature by converting to complementary
 639    * bases, and also update the feature description to match
 640    *
 641    * @param sf
 642    */
 643   static void reverseComplementAlleles(SequenceFeature sf)
 644   {
 645     final String alleles = (String) sf.getValue(ALLELES);
 646     if (alleles == null)
 647     {
 648       return;
 649     }
 650     StringBuilder complement = new StringBuilder(alleles.length());
 651     for (String allele : alleles.split(","))
 652     {
 653       reverseComplementAllele(complement, allele);
 654     }
 655     String comp = complement.toString();
 656     sf.setValue(ALLELES, comp);
 657     sf.setDescription(comp);
 658
 659     /*
 660      * replace value of "alleles=" in sf.ATTRIBUTES as well
 661      * so 'output as GFF' shows reverse complement alleles
 662      */
 663     String atts = sf.getAttributes();
 664     if (atts != null)
 665     {
 666       atts = atts.replace(ALLELES + "=" + alleles, ALLELES + "=" + comp);
 667       sf.setAttributes(atts);
 668     }
 669   }
 670
 671   /**
 672    * Makes the 'reverse complement' of the given allele and appends it to the
 673    * buffer, after a comma separator if not the first
 674    *
 675    * @param complement
 676    * @param allele
 677    */
 678   static void reverseComplementAllele(StringBuilder complement,
 679           String allele)
 680   {
 681     if (complement.length() > 0)
 682     {
 683       complement.append(",");
 684     }
 685
 686     /*
 687      * some 'alleles' are actually descriptive terms
 688      * e.g. HGMD_MUTATION, PhenCode_variation
 689      * - we don't want to 'reverse complement' these
 690      */
 691     if (!Comparison.isNucleotideSequence(allele, true))
 692     {
 693       complement.append(allele);
 694     }
 695     else
 696     {
 697       for (int i = allele.length() - 1; i >= 0; i--)
 698       {
 699         complement.append(Dna.getComplement(allele.charAt(i)));
 700       }
 701     }
 702   }
 703
 704   /**
 705    * Transfers features from sourceSequence to targetSequence
 706    *
 707    * @param accessionId
 708    * @param sourceSequence
 709    * @param targetSequence
 710    * @return true if any features were transferred, else false
 711    */
 712   protected boolean transferFeatures(String accessionId,
 713           SequenceI sourceSequence, SequenceI targetSequence)
 714   {
 715     if (sourceSequence == null || targetSequence == null)
 716     {
 717       return false;
 718     }
 719
 720     // long start = System.currentTimeMillis();
 721     SequenceFeature[] sfs = sourceSequence.getSequenceFeatures();
 722     MapList mapping = getGenomicRangesFromFeatures(sourceSequence, accessionId,
 723             targetSequence.getStart());
 724     if (mapping == null)
 725     {
 726       return false;
 727     }
 728
 729     boolean result = transferFeatures(sfs, targetSequence, mapping,
 730             accessionId);
 731     // System.out.println("transferFeatures (" + (sfs.length) + " --> "
 732     // + targetSequence.getSequenceFeatures().length + ") to "
 733     // + targetSequence.getName()
 734     // + " took " + (System.currentTimeMillis() - start) + "ms");
 735     return result;
 736   }
 737
 738   /**
 739    * Transfer features to the target sequence. The start/end positions are
 740    * converted using the mapping. Features which do not overlap are ignored.
 741    * Features whose parent is not the specified identifier are also ignored.
 742    *
 743    * @param features
 744    * @param targetSequence
 745    * @param mapping
 746    * @param parentId
 747    * @return
 748    */
 749   protected boolean transferFeatures(SequenceFeature[] features,
 750           SequenceI targetSequence, MapList mapping, String parentId)
 751   {
 752     final boolean forwardStrand = mapping.isFromForwardStrand();
 753
 754     /*
 755      * sort features by start position (which corresponds to end
 756      * position descending if reverse strand) so as to add them in
 757      * 'forwards' order to the target sequence
 758      */
 759     sortFeatures(features, forwardStrand);
 760
 761     boolean transferred = false;
 762     for (SequenceFeature sf : features)
 763     {
 764       if (retainFeature(sf, parentId))
 765       {
 766         transferFeature(sf, targetSequence, mapping, forwardStrand);
 767         transferred = true;
 768       }
 769     }
 770     return transferred;
 771   }
 772
 773   /**
 774    * Sort features by start position ascending (if on forward strand), or end
 775    * position descending (if on reverse strand)
 776    *
 777    * @param features
 778    * @param forwardStrand
 779    */
 780   protected static void sortFeatures(SequenceFeature[] features,
 781           final boolean forwardStrand)
 782   {
 783     Arrays.sort(features, new Comparator<SequenceFeature>()
 784     {
 785       @Override
 786       public int compare(SequenceFeature o1, SequenceFeature o2)
 787       {
 788         if (forwardStrand)
 789         {
 790           return Integer.compare(o1.getBegin(), o2.getBegin());
 791         }
 792         else
 793         {
 794           return Integer.compare(o2.getEnd(), o1.getEnd());
 795         }
 796       }
 797     });
 798   }
 799
 800   /**
 801    * Answers true if the feature type is one we want to keep for the sequence.
 802    * Some features are only retrieved in order to identify the sequence range,
 803    * and may then be discarded as redundant information (e.g. "CDS" feature for
 804    * a CDS sequence).
 805    */
 806   @SuppressWarnings("unused")
 807   protected boolean retainFeature(SequenceFeature sf, String accessionId)
 808   {
 809     return true; // override as required
 810   }
 811
 812   /**
 813    * Answers true if the feature has a Parent which refers to the given
 814    * accession id, or if the feature has no parent. Answers false if the
 815    * feature's Parent is for a different accession id.
 816    *
 817    * @param sf
 818    * @param identifier
 819    * @return
 820    */
 821   protected boolean featureMayBelong(SequenceFeature sf, String identifier)
 822   {
 823     String parent = (String) sf.getValue(PARENT);
 824     // using contains to allow for prefix "gene:", "transcript:" etc
 825     if (parent != null && !parent.contains(identifier))
 826     {
 827       // this genomic feature belongs to a different transcript
 828       return false;
 829     }
 830     return true;
 831   }
 832
 833   @Override
 834   public String getDescription()
 835   {
 836     return "Ensembl " + getSourceEnsemblType().getType()
 837             + " sequence with variant features";
 838   }
 839
 840   /**
 841    * Returns a (possibly empty) list of features on the sequence which have the
 842    * specified sequence ontology type (or a sub-type of it), and the given
 843    * identifier as parent
 844    *
 845    * @param sequence
 846    * @param type
 847    * @param parentId
 848    * @return
 849    */
 850   protected List<SequenceFeature> findFeatures(SequenceI sequence,
 851           String type, String parentId)
 852   {
 853     List<SequenceFeature> result = new ArrayList<SequenceFeature>();
 854
 855     SequenceFeature[] sfs = sequence.getSequenceFeatures();
 856     if (sfs != null) {
 857       SequenceOntologyI so = SequenceOntologyFactory.getInstance();
 858       for (SequenceFeature sf :sfs) {
 859         if (so.isA(sf.getType(), type))
 860         {
 861           String parent = (String) sf.getValue(PARENT);
 862           if (parent.equals(parentId))
 863           {
 864             result.add(sf);
 865           }
 866         }
 867       }
 868     }
 869     return result;
 870   }
 871
 872   /**
 873    * Answers true if the feature type is either 'NMD_transcript_variant' or
 874    * 'transcript' or one of its sub-types in the Sequence Ontology. This is
 875    * needed because NMD_transcript_variant behaves like 'transcript' in Ensembl
 876    * although strictly speaking it is not (it is a sub-type of
 877    * sequence_variant).
 878    *
 879    * @param featureType
 880    * @return
 881    */
 882   public static boolean isTranscript(String featureType)
 883   {
 884     return SequenceOntologyI.NMD_TRANSCRIPT_VARIANT.equals(featureType)
 885             || SequenceOntologyFactory.getInstance().isA(featureType,
 886                     SequenceOntologyI.TRANSCRIPT);
 887   }
 888 }