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