src/jalview/ext/ensembl/EnsemblSeqProxy.java

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