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    * Constructor
  99    */
 100   public EnsemblSeqProxy()
 101   {
 102   }
 103
 104   /**
 105    * Makes the sequence queries to Ensembl's REST service and returns an
 106    * alignment consisting of the returned sequences
 107    */
 108   @Override
 109   public AlignmentI getSequenceRecords(String query) throws Exception
 110   {
 111     long now = System.currentTimeMillis();
 112     // TODO use a String... query vararg instead?
 113
 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;
 118     inProgress = true;
 119
 120     /*
 121      * execute queries, if necessary in batches of the
 122      * maximum allowed number of ids
 123      */
 124     int maxQueryCount = getMaximumQueryCount();
 125     for (int v = 0, vSize = allIds.size(); v < vSize; v += maxQueryCount)
 126     {
 127       int p = Math.min(vSize, v + maxQueryCount);
 128       List<String> ids = allIds.subList(v, p);
 129       try
 130       {
 131         alignment = fetchSequences(ids, alignment);
 132       } catch (Throwable r)
 133       {
 134         inProgress = false;
 135         String msg = "Aborting ID retrieval after " + v
 136                 + " chunks. Unexpected problem (" + r.getLocalizedMessage()
 137                 + ")";
 138         System.err.println(msg);
 139         if (alignment != null)
 140         {
 141           break; // return what we got
 142         }
 143         else
 144         {
 145           throw new JalviewException(msg, r);
 146         }
 147       }
 148     }
 149
 150     /*
 151      * fetch and transfer genomic sequence features
 152      */
 153     for (String accId : allIds)
 154     {
 155       addFeaturesAndProduct(accId, alignment);
 156     }
 157
 158     inProgress = false;
 159     System.out.println(getClass().getName() + " took "
 160             + (System.currentTimeMillis() - now) + "ms to fetch");
 161     return alignment;
 162   }
 163
 164   /**
 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
 168    * the dna sequence.
 169    *
 170    * @param accId
 171    * @param alignment
 172    */
 173   protected void addFeaturesAndProduct(String accId, AlignmentI alignment)
 174   {
 175     try
 176     {
 177       /*
 178        * get 'dummy' genomic sequence with exon, cds and variation features
 179        */
 180       EnsemblOverlap gffFetcher = new EnsemblOverlap();
 181       EnsemblFeatureType[] features = getFeaturesToFetch();
 182       AlignmentI geneFeatures = gffFetcher.getSequenceRecords(accId,
 183               features);
 184       if (geneFeatures.getHeight() > 0)
 185       {
 186         /*
 187          * transfer features to the query sequence
 188          */
 189         SequenceI genomicSequence = geneFeatures.getSequenceAt(0);
 190         SequenceI querySeq = alignment.findName(accId);
 191         transferFeatures(accId, genomicSequence, querySeq);
 192
 193         /*
 194          * fetch and map protein product, and add it as a cross-reference
 195          * of the retrieved sequence
 196          */
 197         addProteinProduct(querySeq);
 198       }
 199     } catch (IOException e)
 200     {
 201       System.err.println("Error transferring Ensembl features: "
 202               + e.getMessage());
 203     }
 204   }
 205
 206   /**
 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).
 211    *
 212    * @return
 213    */
 214   protected abstract EnsemblFeatureType[] getFeaturesToFetch();
 215
 216   /**
 217    * Fetches and maps the protein product, and adds it as a cross-reference of
 218    * the retrieved sequence
 219    */
 220   protected void addProteinProduct(SequenceI querySeq)
 221   {
 222     String accId = querySeq.getName();
 223     try
 224     {
 225       AlignmentI protein = new EnsemblProtein().getSequenceRecords(accId);
 226       if (protein == null || protein.getHeight() == 0)
 227       {
 228         System.out.println("Failed to retrieve protein for " + accId);
 229         return;
 230       }
 231       SequenceI proteinSeq = protein.getSequenceAt(0);
 232
 233       /*
 234        * need dataset sequences (to be the subject of mappings)
 235        */
 236       proteinSeq.createDatasetSequence();
 237       querySeq.createDatasetSequence();
 238
 239       MapList mapList = mapCdsToProtein(querySeq, proteinSeq);
 240       if (mapList != null)
 241       {
 242         Mapping map = new Mapping(proteinSeq.getDatasetSequence(), mapList);
 243         DBRefEntry dbr = new DBRefEntry(getDbSource(), getDbVersion(),
 244                 accId, map);
 245         querySeq.getDatasetSequence().addDBRef(dbr);
 246       }
 247     } catch (Exception e)
 248     {
 249       System.err
 250               .println(String.format("Error retrieving protein for %s: %s",
 251                       accId, e.getMessage()));
 252     }
 253   }
 254
 255   /**
 256    * Returns a mapping from dna to protein by inspecting sequence features of
 257    * type "CDS" on the dna.
 258    *
 259    * @param dnaSeq
 260    * @param proteinSeq
 261    * @return
 262    */
 263   protected MapList mapCdsToProtein(SequenceI dnaSeq, SequenceI proteinSeq)
 264   {
 265     SequenceFeature[] sfs = dnaSeq.getSequenceFeatures();
 266     if (sfs == null)
 267     {
 268       return null;
 269     }
 270
 271     List<int[]> ranges = new ArrayList<int[]>(50);
 272     SequenceOntology so = SequenceOntology.getInstance();
 273
 274     int mappedDnaLength = 0;
 275
 276     /*
 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.
 280      */
 281     for (SequenceFeature sf : sfs)
 282     {
 283       /*
 284        * process a CDS feature (or a sub-type of CDS)
 285        */
 286       if (so.isA(sf.getType(), SequenceOntology.CDS))
 287       {
 288         ranges.add(new int[] { sf.getBegin(), sf.getEnd() });
 289         mappedDnaLength += Math.abs(sf.getEnd() - sf.getBegin()) + 1;
 290       }
 291     }
 292     int proteinLength = proteinSeq.getLength();
 293     List<int[]> proteinRange = new ArrayList<int[]>();
 294     proteinRange.add(new int[] { 1, proteinLength });
 295
 296     /*
 297      * dna length should map to protein (or protein minus stop codon)
 298      */
 299     if (mappedDnaLength == 3 * proteinLength
 300             || mappedDnaLength == 3 * (proteinLength + 1))
 301     {
 302       return new MapList(ranges, proteinRange, 3, 1);
 303     }
 304     return null;
 305   }
 306
 307   /**
 308    * Fetches sequences for the list of accession ids and adds them to the
 309    * alignment. Returns the extended (or created) alignment.
 310    *
 311    * @param ids
 312    * @param alignment
 313    * @return
 314    * @throws JalviewException
 315    * @throws IOException
 316    */
 317   protected AlignmentI fetchSequences(List<String> ids, AlignmentI alignment)
 318           throws JalviewException, IOException
 319   {
 320     if (!isEnsemblAvailable())
 321     {
 322       inProgress = false;
 323       throw new JalviewException("ENSEMBL Rest API not available.");
 324     }
 325     FileParse fp = getSequenceReader(ids);
 326     FastaFile fr = new FastaFile(fp);
 327     if (fr.hasWarningMessage())
 328     {
 329       System.out.println(String.format(
 330               "Warning when retrieving %d ids %s\n%s", ids.size(),
 331               ids.toString(), fr.getWarningMessage()));
 332     }
 333     else if (fr.getSeqs().size() != ids.size())
 334     {
 335       System.out.println(String.format(
 336               "Only retrieved %d sequences for %d query strings", fr
 337                       .getSeqs().size(), ids.size()));
 338     }
 339     if (fr.getSeqs().size() > 0)
 340     {
 341       AlignmentI seqal = new Alignment(
 342               fr.getSeqsAsArray());
 343       for (SequenceI sq:seqal.getSequences())
 344       {
 345         if (sq.getDescription() == null)
 346         {
 347           sq.setDescription(getDbName());
 348         }
 349         String name = sq.getName();
 350         if (ids.contains(name)
 351                 || ids.contains(name.replace("ENSP", "ENST")))
 352         {
 353           DBRefUtils.parseToDbRef(sq, DBRefSource.ENSEMBL, "0", name);
 354         }
 355       }
 356       if (alignment == null)
 357       {
 358         alignment = seqal;
 359       }
 360       else
 361       {
 362         alignment.append(seqal);
 363       }
 364     }
 365     return alignment;
 366   }
 367
 368   /**
 369    * Returns the URL for the REST call
 370    *
 371    * @return
 372    * @throws MalformedURLException
 373    */
 374   @Override
 375   protected URL getUrl(List<String> ids) throws MalformedURLException
 376   {
 377     // ids are not used - they go in the POST body instead
 378     StringBuffer urlstring = new StringBuffer(128);
 379     urlstring.append(SEQUENCE_ID_URL);
 380
 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"));
 384
 385     URL url = new URL(urlstring.toString());
 386     return url;
 387   }
 388
 389   /**
 390    * A sequence/id POST request currently allows up to 50 queries
 391    *
 392    * @see http://rest.ensembl.org/documentation/info/sequence_id_post
 393    */
 394   @Override
 395   public int getMaximumQueryCount()
 396   {
 397     return 50;
 398   }
 399
 400   @Override
 401   public boolean useGetRequest()
 402   {
 403     return false;
 404   }
 405
 406   @Override
 407   public String getRequestMimeType()
 408   {
 409     return "application/json";
 410   }
 411
 412   @Override
 413   public String getResponseMimeType()
 414   {
 415     return "text/x-fasta";
 416   }
 417
 418   /**
 419    *
 420    * @return the configured sequence return type for this source
 421    */
 422   protected abstract EnsemblSeqType getSourceEnsemblType();
 423
 424   /**
 425    * Returns a list of [start, end] genomic ranges corresponding to the sequence
 426    * being retrieved.
 427    *
 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
 431    * <ul>
 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>
 438    * </ul>
 439    *
 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).
 443    *
 444    * @param sfs
 445    * @param accId
 446    * @return
 447    */
 448   protected MapList getGenomicRanges(SequenceFeature[] sfs, String accId)
 449   {
 450     /*
 451      * generously size for initial number of cds regions
 452      * (worst case titin Q8WZ42 has c. 313 exons)
 453      */
 454     List<int[]> regions = new ArrayList<int[]>(100);
 455     int mappedLength = 0;
 456     int direction = 1; // forward
 457     boolean directionSet = false;
 458
 459     for (SequenceFeature sf : sfs)
 460     {
 461       /*
 462        * accept the target feature type or a specialisation of it
 463        * (e.g. coding_exon for exon)
 464        */
 465       if (identifiesSequence(sf, accId))
 466       {
 467           int strand = sf.getStrand();
 468
 469           if (directionSet && strand != direction)
 470           {
 471             // abort - mix of forward and backward
 472           System.err.println("Error: forward and backward strand for "
 473                   + accId);
 474             return null;
 475           }
 476           direction = strand;
 477           directionSet = true;
 478
 479           /*
 480            * add to CDS ranges, semi-sorted forwards/backwards
 481            */
 482           if (strand < 0)
 483           {
 484             regions.add(0, new int[] { sf.getEnd(), sf.getBegin() });
 485           }
 486           else
 487           {
 488             regions.add(new int[] { sf.getBegin(), sf.getEnd() });
 489           }
 490           mappedLength += Math.abs(sf.getEnd() - sf.getBegin() + 1);
 491         }
 492     }
 493
 494     /*
 495      * a final sort is needed since Ensembl returns CDS sorted within source
 496      * (havana / ensembl_havana)
 497      */
 498     Collections.sort(regions, new RangeSorter(direction == 1));
 499
 500     List<int[]> to = new ArrayList<int[]>();
 501     to.add(new int[] { 1, mappedLength });
 502
 503     return new MapList(regions, to, 1, 1);
 504   }
 505
 506   /**
 507    * Returns true if the sequence feature identifies positions of the genomic
 508    * sequence feature which are within the sequence being retrieved.
 509    *
 510    * @param sf
 511    * @param accId
 512    * @return
 513    */
 514   protected abstract boolean identifiesSequence(SequenceFeature sf,
 515           String accId);
 516
 517   /**
 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.
 522    *
 523    * @param sf
 524    * @param targetSequence
 525    * @param overlap
 526    */
 527   protected void transferFeature(SequenceFeature sf,
 528           SequenceI targetSequence, MapList overlap)
 529   {
 530     String parent = (String) sf.getValue(PARENT);
 531     if (parent != null && !parent.contains(targetSequence.getName()))
 532     {
 533       // this genomic feature belongs to a different transcript
 534       return;
 535     }
 536
 537     int start = sf.getBegin();
 538     int end = sf.getEnd();
 539     int[] mappedRange = overlap.locateInTo(start, end);
 540
 541     if (mappedRange != null)
 542     {
 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);
 548
 549       /*
 550        * for sequence_variant, make an additional feature with consequence
 551        */
 552       if (SequenceOntology.getInstance().isSequenceVariant(sf.getType()))
 553       {
 554         String consequence = (String) sf.getValue("consequence_type");
 555         if (consequence != null)
 556         {
 557           SequenceFeature sf2 = new SequenceFeature("consequence",
 558                   consequence, copy.getBegin(), copy.getEnd(), 0f,
 559                   null);
 560           targetSequence.addSequenceFeature(sf2);
 561         }
 562       }
 563     }
 564
 565   }
 566
 567   /**
 568    * Transfers features from sourceSequence to targetSequence
 569    *
 570    * @param accessionId
 571    * @param sourceSequence
 572    * @param targetSequence
 573    */
 574   protected void transferFeatures(String accessionId,
 575           SequenceI sourceSequence, SequenceI targetSequence)
 576   {
 577     if (sourceSequence == null || targetSequence == null)
 578     {
 579       return;
 580     }
 581
 582     SequenceFeature[] sfs = sourceSequence.getSequenceFeatures();
 583     MapList overlap = getGenomicRanges(sfs, accessionId);
 584
 585     final boolean forwardStrand = overlap.isFromForwardStrand();
 586
 587     /*
 588      * sort features by start position (descending if reverse strand)
 589      * before transferring (in forwards order) to the target sequence
 590      */
 591     Arrays.sort(sfs, new Comparator<SequenceFeature>()
 592     {
 593       @Override
 594       public int compare(SequenceFeature o1, SequenceFeature o2)
 595       {
 596         int c = Integer.compare(o1.getBegin(), o2.getBegin());
 597         return forwardStrand ? c : -c;
 598       }
 599     });
 600
 601     for (SequenceFeature sf : sfs)
 602     {
 603       if (retainFeature(sf.getType()))
 604       {
 605         transferFeature(sf, targetSequence, overlap);
 606       }
 607     }
 608   }
 609
 610   /**
 611    * Answers true if the feature type is one to attach to the retrieved sequence
 612    *
 613    * @param type
 614    * @return
 615    */
 616   protected boolean retainFeature(@SuppressWarnings("unused") String type)
 617   {
 618     return true; // default is to keep all
 619   }
 620 }