src/jalview/ext/ensembl/EnsemblSeqProxy.java

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