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