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