src/jalview/io/vcf/VCFLoader.java

   1 package jalview.io.vcf;
   2
   3 import jalview.analysis.Dna;
   4 import jalview.api.AlignViewControllerGuiI;
   5 import jalview.bin.Cache;
   6 import jalview.datamodel.DBRefEntry;
   7 import jalview.datamodel.GeneLociI;
   8 import jalview.datamodel.Mapping;
   9 import jalview.datamodel.SequenceFeature;
  10 import jalview.datamodel.SequenceI;
  11 import jalview.datamodel.features.FeatureAttributeType;
  12 import jalview.datamodel.features.FeatureSource;
  13 import jalview.datamodel.features.FeatureSources;
  14 import jalview.ext.ensembl.EnsemblMap;
  15 import jalview.ext.htsjdk.HtsContigDb;
  16 import jalview.ext.htsjdk.VCFReader;
  17 import jalview.io.gff.Gff3Helper;
  18 import jalview.io.gff.SequenceOntologyI;
  19 import jalview.util.MapList;
  20 import jalview.util.MappingUtils;
  21 import jalview.util.MessageManager;
  22
  23 import java.io.File;
  24 import java.io.IOException;
  25 import java.io.UnsupportedEncodingException;
  26 import java.net.URLDecoder;
  27 import java.util.ArrayList;
  28 import java.util.HashMap;
  29 import java.util.List;
  30 import java.util.Map;
  31 import java.util.Map.Entry;
  32 import java.util.regex.Pattern;
  33 import java.util.regex.PatternSyntaxException;
  34
  35 import htsjdk.samtools.SAMException;
  36 import htsjdk.samtools.SAMSequenceDictionary;
  37 import htsjdk.samtools.SAMSequenceRecord;
  38 import htsjdk.samtools.util.CloseableIterator;
  39 import htsjdk.variant.variantcontext.Allele;
  40 import htsjdk.variant.variantcontext.VariantContext;
  41 import htsjdk.variant.vcf.VCFHeader;
  42 import htsjdk.variant.vcf.VCFHeaderLine;
  43 import htsjdk.variant.vcf.VCFHeaderLineCount;
  44 import htsjdk.variant.vcf.VCFHeaderLineType;
  45 import htsjdk.variant.vcf.VCFInfoHeaderLine;
  46
  47 /**
  48  * A class to read VCF data (using the htsjdk) and add variants as sequence
  49  * features on dna and any related protein product sequences
  50  *
  51  * @author gmcarstairs
  52  */
  53 public class VCFLoader
  54 {
  55   private static final String ENCODED_COMMA = "%2C";
  56
  57   private static final String ENCODED_PERCENT = "%25";
  58
  59   private static final String ENCODED_EQUALS = "%3D";
  60
  61   private static final String ENCODED_SEMICOLON = "%3B";
  62
  63   private static final String ENCODED_COLON = "%3A";
  64
  65   private static final String UTF_8 = "UTF-8";
  66
  67   private static final String DEFAULT_SPECIES = "homo_sapiens";
  68
  69   /**
  70    * A class to model the mapping from sequence to VCF coordinates. Cases include
  71    * <ul>
  72    * <li>a direct 1:1 mapping where the sequence is one of the VCF contigs</li>
  73    * <li>a mapping of sequence to chromosomal coordinates, where sequence and VCF
  74    * use the same reference assembly</li>
  75    * <li>a modified mapping of sequence to chromosomal coordinates, where sequence
  76    * and VCF use different reference assembles</li>
  77    * </ul>
  78    */
  79   class VCFMap
  80   {
  81     final String chromosome;
  82
  83     final MapList map;
  84
  85     VCFMap(String chr, MapList m)
  86     {
  87       chromosome = chr;
  88       map = m;
  89     }
  90
  91     @Override
  92     public String toString()
  93     {
  94       return chromosome + ":" + map.toString();
  95     }
  96   }
  97
  98   /*
  99    * Lookup keys, and default values, for Preference entries that describe
 100    * patterns for VCF and VEP fields to capture
 101    */
 102   private static final String VEP_FIELDS_PREF = "VEP_FIELDS";
 103
 104   private static final String VCF_FIELDS_PREF = "VCF_FIELDS";
 105
 106   private static final String DEFAULT_VCF_FIELDS = ".*";
 107
 108   private static final String DEFAULT_VEP_FIELDS = ".*";// "Allele,Consequence,IMPACT,SWISSPROT,SIFT,PolyPhen,CLIN_SIG";
 109
 110   /*
 111    * Lookup keys, and default values, for Preference entries that give
 112    * mappings from tokens in the 'reference' header to species or assembly
 113    */
 114   private static final String VCF_ASSEMBLY = "VCF_ASSEMBLY";
 115
 116   private static final String DEFAULT_VCF_ASSEMBLY = "assembly19=GRCh37,hs37=GRCh37,grch37=GRCh37,grch38=GRCh38";
 117
 118   private static final String VCF_SPECIES = "VCF_SPECIES"; // default is human
 119
 120   private static final String DEFAULT_REFERENCE = "grch37"; // fallback default is human GRCh37
 121
 122   /*
 123    * keys to fields of VEP CSQ consequence data
 124    * see https://www.ensembl.org/info/docs/tools/vep/vep_formats.html
 125    */
 126   private static final String CSQ_CONSEQUENCE_KEY = "Consequence";
 127   private static final String CSQ_ALLELE_KEY = "Allele";
 128   private static final String CSQ_ALLELE_NUM_KEY = "ALLELE_NUM"; // 0 (ref), 1...
 129   private static final String CSQ_FEATURE_KEY = "Feature"; // Ensembl stable id
 130
 131   /*
 132    * default VCF INFO key for VEP consequence data
 133    * NB this can be overridden running VEP with --vcf_info_field
 134    * - we don't handle this case (require identifier to be CSQ)
 135    */
 136   private static final String CSQ_FIELD = "CSQ";
 137
 138   /*
 139    * separator for fields in consequence data is '|'
 140    */
 141   private static final String PIPE_REGEX = "\\|";
 142
 143   /*
 144    * delimiter that separates multiple consequence data blocks
 145    */
 146   private static final String COMMA = ",";
 147
 148   /*
 149    * the feature group assigned to a VCF variant in Jalview
 150    */
 151   private static final String FEATURE_GROUP_VCF = "VCF";
 152
 153   /*
 154    * internal delimiter used to build keys for assemblyMappings
 155    *
 156    */
 157   private static final String EXCL = "!";
 158
 159   /*
 160    * the VCF file we are processing
 161    */
 162   protected String vcfFilePath;
 163
 164   /*
 165    * mappings between VCF and sequence reference assembly regions, as
 166    * key = "species!chromosome!fromAssembly!toAssembly
 167    * value = Map{fromRange, toRange}
 168    */
 169   private Map<String, Map<int[], int[]>> assemblyMappings;
 170
 171   private VCFReader reader;
 172
 173   /*
 174    * holds details of the VCF header lines (metadata)
 175    */
 176   private VCFHeader header;
 177
 178   /*
 179    * species (as a valid Ensembl term) the VCF is for
 180    */
 181   private String vcfSpecies;
 182
 183   /*
 184    * genome assembly version (as a valid Ensembl identifier) the VCF is for
 185    */
 186   private String vcfAssembly;
 187
 188   /*
 189    * a Dictionary of contigs (if present) referenced in the VCF file
 190    */
 191   private SAMSequenceDictionary dictionary;
 192
 193   /*
 194    * the position (0...) of field in each block of
 195    * CSQ (consequence) data (if declared in the VCF INFO header for CSQ)
 196    * see http://www.ensembl.org/info/docs/tools/vep/vep_formats.html
 197    */
 198   private int csqConsequenceFieldIndex = -1;
 199   private int csqAlleleFieldIndex = -1;
 200   private int csqAlleleNumberFieldIndex = -1;
 201   private int csqFeatureFieldIndex = -1;
 202
 203   // todo the same fields for SnpEff ANN data if wanted
 204   // see http://snpeff.sourceforge.net/SnpEff_manual.html#input
 205
 206   /*
 207    * a unique identifier under which to save metadata about feature
 208    * attributes (selected INFO field data)
 209    */
 210   private String sourceId;
 211
 212   /*
 213    * The INFO IDs of data that is both present in the VCF file, and
 214    * also matched by any filters for data of interest
 215    */
 216   List<String> vcfFieldsOfInterest;
 217
 218   /*
 219    * The field offsets and identifiers for VEP (CSQ) data that is both present
 220    * in the VCF file, and also matched by any filters for data of interest
 221    * for example 0 -> Allele, 1 -> Consequence, ..., 36 -> SIFT, ...
 222    */
 223   Map<Integer, String> vepFieldsOfInterest;
 224
 225   /**
 226    * Constructor given a VCF file
 227    *
 228    * @param alignment
 229    */
 230   public VCFLoader(String vcfFile)
 231   {
 232     try
 233     {
 234       initialise(vcfFile);
 235     } catch (IOException e)
 236     {
 237       System.err.println("Error opening VCF file: " + e.getMessage());
 238     }
 239
 240     // map of species!chromosome!fromAssembly!toAssembly to {fromRange, toRange}
 241     assemblyMappings = new HashMap<>();
 242   }
 243
 244   /**
 245    * Starts a new thread to query and load VCF variant data on to the given
 246    * sequences
 247    * <p>
 248    * This method is not thread safe - concurrent threads should use separate
 249    * instances of this class.
 250    *
 251    * @param seqs
 252    * @param gui
 253    */
 254   public void loadVCF(SequenceI[] seqs, final AlignViewControllerGuiI gui)
 255   {
 256     if (gui != null)
 257     {
 258       gui.setStatus(MessageManager.getString("label.searching_vcf"));
 259     }
 260
 261     new Thread()
 262     {
 263       @Override
 264       public void run()
 265       {
 266         VCFLoader.this.doLoad(seqs, gui);
 267       }
 268     }.start();
 269   }
 270
 271   /**
 272    * Reads the specified contig sequence and adds its VCF variants to it
 273    *
 274    * @param contig
 275    *          the id of a single sequence (contig) to load
 276    * @return
 277    */
 278   public SequenceI loadVCFContig(String contig)
 279   {
 280     VCFHeaderLine headerLine = header.getOtherHeaderLine(VCFHeader.REFERENCE_KEY);
 281     if (headerLine == null)
 282     {
 283       Cache.log.error("VCF reference header not found");
 284       return null;
 285     }
 286     String ref = headerLine.getValue();
 287     if (ref.startsWith("file://"))
 288     {
 289       ref = ref.substring(7);
 290     }
 291     setSpeciesAndAssembly(ref);
 292
 293     SequenceI seq = null;
 294     File dbFile = new File(ref);
 295
 296     if (dbFile.exists())
 297     {
 298       HtsContigDb db = new HtsContigDb("", dbFile);
 299       seq = db.getSequenceProxy(contig);
 300       loadSequenceVCF(seq);
 301       db.close();
 302     }
 303     else
 304     {
 305       Cache.log.error("VCF reference not found: " + ref);
 306     }
 307
 308     return seq;
 309   }
 310
 311   /**
 312    * Loads VCF on to one or more sequences
 313    *
 314    * @param seqs
 315    * @param gui
 316    *          optional callback handler for messages
 317    */
 318   protected void doLoad(SequenceI[] seqs, AlignViewControllerGuiI gui)
 319   {
 320     try
 321     {
 322       VCFHeaderLine ref = header
 323               .getOtherHeaderLine(VCFHeader.REFERENCE_KEY);
 324       String reference = ref == null ? null : ref.getValue();
 325
 326       setSpeciesAndAssembly(reference);
 327
 328       int varCount = 0;
 329       int seqCount = 0;
 330
 331       /*
 332        * query for VCF overlapping each sequence in turn
 333        */
 334       for (SequenceI seq : seqs)
 335       {
 336         int added = loadSequenceVCF(seq);
 337         if (added > 0)
 338         {
 339           seqCount++;
 340           varCount += added;
 341           transferAddedFeatures(seq);
 342         }
 343       }
 344       if (gui != null)
 345       {
 346         String msg = MessageManager.formatMessage("label.added_vcf",
 347                 varCount, seqCount);
 348         gui.setStatus(msg);
 349         if (gui.getFeatureSettingsUI() != null)
 350         {
 351           gui.getFeatureSettingsUI().discoverAllFeatureData();
 352         }
 353       }
 354     } catch (Throwable e)
 355     {
 356       System.err.println("Error processing VCF: " + e.getMessage());
 357       e.printStackTrace();
 358       if (gui != null)
 359       {
 360         gui.setStatus("Error occurred - see console for details");
 361       }
 362     } finally
 363     {
 364       if (reader != null)
 365       {
 366         try
 367         {
 368           reader.close();
 369         } catch (IOException e)
 370         {
 371           // ignore
 372         }
 373       }
 374       header = null;
 375       dictionary = null;
 376     }
 377   }
 378
 379   /**
 380    * Attempts to determine and save the species and genome assembly version to
 381    * which the VCF data applies. This may be done by parsing the {@code reference}
 382    * header line, configured in a property file, or (potentially) confirmed
 383    * interactively by the user.
 384    * <p>
 385    * The saved values should be identifiers valid for Ensembl's REST service
 386    * {@code map} endpoint, so they can be used (if necessary) to retrieve the
 387    * mapping between VCF coordinates and sequence coordinates.
 388    *
 389    * @param reference
 390    * @see https://rest.ensembl.org/documentation/info/assembly_map
 391    * @see https://rest.ensembl.org/info/assembly/human?content-type=text/xml
 392    * @see https://rest.ensembl.org/info/species?content-type=text/xml
 393    */
 394   protected void setSpeciesAndAssembly(String reference)
 395   {
 396     if (reference == null)
 397     {
 398       Cache.log.error("No VCF ##reference found, defaulting to "
 399               + DEFAULT_REFERENCE + ":" + DEFAULT_SPECIES);
 400       reference = DEFAULT_REFERENCE; // default to GRCh37 if not specified
 401     }
 402     reference = reference.toLowerCase();
 403
 404     /*
 405      * for a non-human species, or other assembly identifier,
 406      * specify as a Jalview property file entry e.g.
 407      * VCF_ASSEMBLY = hs37=GRCh37,assembly19=GRCh37
 408      * VCF_SPECIES = c_elegans=celegans
 409      * to map a token in the reference header to a value
 410      */
 411     String prop = Cache.getDefault(VCF_ASSEMBLY, DEFAULT_VCF_ASSEMBLY);
 412     for (String token : prop.split(","))
 413     {
 414       String[] tokens = token.split("=");
 415       if (tokens.length == 2)
 416       {
 417         if (reference.contains(tokens[0].trim().toLowerCase()))
 418         {
 419           vcfAssembly = tokens[1].trim();
 420           break;
 421         }
 422       }
 423     }
 424
 425     vcfSpecies = DEFAULT_SPECIES;
 426     prop = Cache.getProperty(VCF_SPECIES);
 427     if (prop != null)
 428     {
 429       for (String token : prop.split(","))
 430       {
 431         String[] tokens = token.split("=");
 432         if (tokens.length == 2)
 433         {
 434           if (reference.contains(tokens[0].trim().toLowerCase()))
 435           {
 436             vcfSpecies = tokens[1].trim();
 437             break;
 438           }
 439         }
 440       }
 441     }
 442   }
 443
 444   /**
 445    * Opens the VCF file and parses header data
 446    *
 447    * @param filePath
 448    * @throws IOException
 449    */
 450   private void initialise(String filePath) throws IOException
 451   {
 452     vcfFilePath = filePath;
 453
 454     reader = new VCFReader(filePath);
 455
 456     header = reader.getFileHeader();
 457
 458     try
 459     {
 460       dictionary = header.getSequenceDictionary();
 461     } catch (SAMException e)
 462     {
 463       // ignore - thrown if any contig line lacks length info
 464     }
 465
 466     sourceId = filePath;
 467
 468     saveMetadata(sourceId);
 469
 470     /*
 471      * get offset of CSQ ALLELE_NUM and Feature if declared
 472      */
 473     parseCsqHeader();
 474   }
 475
 476   /**
 477    * Reads metadata (such as INFO field descriptions and datatypes) and saves
 478    * them for future reference
 479    *
 480    * @param theSourceId
 481    */
 482   void saveMetadata(String theSourceId)
 483   {
 484     List<Pattern> vcfFieldPatterns = getFieldMatchers(VCF_FIELDS_PREF,
 485             DEFAULT_VCF_FIELDS);
 486     vcfFieldsOfInterest = new ArrayList<>();
 487
 488     FeatureSource metadata = new FeatureSource(theSourceId);
 489
 490     for (VCFInfoHeaderLine info : header.getInfoHeaderLines())
 491     {
 492       String attributeId = info.getID();
 493       String desc = info.getDescription();
 494       VCFHeaderLineType type = info.getType();
 495       FeatureAttributeType attType = null;
 496       switch (type)
 497       {
 498       case Character:
 499         attType = FeatureAttributeType.Character;
 500         break;
 501       case Flag:
 502         attType = FeatureAttributeType.Flag;
 503         break;
 504       case Float:
 505         attType = FeatureAttributeType.Float;
 506         break;
 507       case Integer:
 508         attType = FeatureAttributeType.Integer;
 509         break;
 510       case String:
 511         attType = FeatureAttributeType.String;
 512         break;
 513       }
 514       metadata.setAttributeName(attributeId, desc);
 515       metadata.setAttributeType(attributeId, attType);
 516
 517       if (isFieldWanted(attributeId, vcfFieldPatterns))
 518       {
 519         vcfFieldsOfInterest.add(attributeId);
 520       }
 521     }
 522
 523     FeatureSources.getInstance().addSource(theSourceId, metadata);
 524   }
 525
 526   /**
 527    * Answers true if the field id is matched by any of the filter patterns, else
 528    * false. Matching is against regular expression patterns, and is not
 529    * case-sensitive.
 530    *
 531    * @param id
 532    * @param filters
 533    * @return
 534    */
 535   private boolean isFieldWanted(String id, List<Pattern> filters)
 536   {
 537     for (Pattern p : filters)
 538     {
 539       if (p.matcher(id.toUpperCase()).matches())
 540       {
 541         return true;
 542       }
 543     }
 544     return false;
 545   }
 546
 547   /**
 548    * Records 'wanted' fields defined in the CSQ INFO header (if there is one).
 549    * Also records the position of selected fields (Allele, ALLELE_NUM, Feature)
 550    * required for processing.
 551    * <p>
 552    * CSQ fields are declared in the CSQ INFO Description e.g.
 553    * <p>
 554    * Description="Consequence ...from ... VEP. Format: Allele|Consequence|...
 555    */
 556   protected void parseCsqHeader()
 557   {
 558     List<Pattern> vepFieldFilters = getFieldMatchers(VEP_FIELDS_PREF,
 559             DEFAULT_VEP_FIELDS);
 560     vepFieldsOfInterest = new HashMap<>();
 561
 562     VCFInfoHeaderLine csqInfo = header.getInfoHeaderLine(CSQ_FIELD);
 563     if (csqInfo == null)
 564     {
 565       return;
 566     }
 567
 568     /*
 569      * parse out the pipe-separated list of CSQ fields; we assume here that
 570      * these form the last part of the description, and contain no spaces
 571      */
 572     String desc = csqInfo.getDescription();
 573     int spacePos = desc.lastIndexOf(" ");
 574     desc = desc.substring(spacePos + 1);
 575
 576     if (desc != null)
 577     {
 578       String[] format = desc.split(PIPE_REGEX);
 579       int index = 0;
 580       for (String field : format)
 581       {
 582         if (CSQ_CONSEQUENCE_KEY.equals(field))
 583         {
 584           csqConsequenceFieldIndex = index;
 585         }
 586         if (CSQ_ALLELE_NUM_KEY.equals(field))
 587         {
 588           csqAlleleNumberFieldIndex = index;
 589         }
 590         if (CSQ_ALLELE_KEY.equals(field))
 591         {
 592           csqAlleleFieldIndex = index;
 593         }
 594         if (CSQ_FEATURE_KEY.equals(field))
 595         {
 596           csqFeatureFieldIndex = index;
 597         }
 598
 599         if (isFieldWanted(field, vepFieldFilters))
 600         {
 601           vepFieldsOfInterest.put(index, field);
 602         }
 603
 604         index++;
 605       }
 606     }
 607   }
 608
 609   /**
 610    * Reads the Preference value for the given key, with default specified if no
 611    * preference set. The value is interpreted as a comma-separated list of
 612    * regular expressions, and converted into a list of compiled patterns ready
 613    * for matching. Patterns are forced to upper-case for non-case-sensitive
 614    * matching.
 615    * <p>
 616    * This supports user-defined filters for fields of interest to capture while
 617    * processing data. For example, VCF_FIELDS = AF,AC* would mean that VCF INFO
 618    * fields with an ID of AF, or starting with AC, would be matched.
 619    *
 620    * @param key
 621    * @param def
 622    * @return
 623    */
 624   private List<Pattern> getFieldMatchers(String key, String def)
 625   {
 626     String pref = Cache.getDefault(key, def);
 627     List<Pattern> patterns = new ArrayList<>();
 628     String[] tokens = pref.split(",");
 629     for (String token : tokens)
 630     {
 631       try
 632       {
 633       patterns.add(Pattern.compile(token.toUpperCase()));
 634       } catch (PatternSyntaxException e)
 635       {
 636         System.err.println("Invalid pattern ignored: " + token);
 637       }
 638     }
 639     return patterns;
 640   }
 641
 642   /**
 643    * Transfers VCF features to sequences to which this sequence has a mapping.
 644    * If the mapping is 3:1, computes peptide variants from nucleotide variants.
 645    *
 646    * @param seq
 647    */
 648   protected void transferAddedFeatures(SequenceI seq)
 649   {
 650     DBRefEntry[] dbrefs = seq.getDBRefs();
 651     if (dbrefs == null)
 652     {
 653       return;
 654     }
 655     for (DBRefEntry dbref : dbrefs)
 656     {
 657       Mapping mapping = dbref.getMap();
 658       if (mapping == null || mapping.getTo() == null)
 659       {
 660         continue;
 661       }
 662
 663       SequenceI mapTo = mapping.getTo();
 664       MapList map = mapping.getMap();
 665       if (map.getFromRatio() == 3)
 666       {
 667         /*
 668          * dna-to-peptide product mapping
 669          */
 670         // JAL-3187 render on the fly instead
 671         // AlignmentUtils.computeProteinFeatures(seq, mapTo, map);
 672       }
 673       else
 674       {
 675         /*
 676          * nucleotide-to-nucleotide mapping e.g. transcript to CDS
 677          */
 678         List<SequenceFeature> features = seq.getFeatures()
 679                 .getPositionalFeatures(SequenceOntologyI.SEQUENCE_VARIANT);
 680         for (SequenceFeature sf : features)
 681         {
 682           if (FEATURE_GROUP_VCF.equals(sf.getFeatureGroup()))
 683           {
 684             transferFeature(sf, mapTo, map);
 685           }
 686         }
 687       }
 688     }
 689   }
 690
 691   /**
 692    * Tries to add overlapping variants read from a VCF file to the given sequence,
 693    * and returns the number of variant features added
 694    *
 695    * @param seq
 696    * @return
 697    */
 698   protected int loadSequenceVCF(SequenceI seq)
 699   {
 700     VCFMap vcfMap = getVcfMap(seq);
 701     if (vcfMap == null)
 702     {
 703       return 0;
 704     }
 705
 706     /*
 707      * work with the dataset sequence here
 708      */
 709     SequenceI dss = seq.getDatasetSequence();
 710     if (dss == null)
 711     {
 712       dss = seq;
 713     }
 714     return addVcfVariants(dss, vcfMap);
 715   }
 716
 717   /**
 718    * Answers a map from sequence coordinates to VCF chromosome ranges
 719    *
 720    * @param seq
 721    * @return
 722    */
 723   private VCFMap getVcfMap(SequenceI seq)
 724   {
 725     /*
 726      * simplest case: sequence has id and length matching a VCF contig
 727      */
 728     VCFMap vcfMap = null;
 729     if (dictionary != null)
 730     {
 731       vcfMap = getContigMap(seq);
 732     }
 733     if (vcfMap != null)
 734     {
 735       return vcfMap;
 736     }
 737
 738     /*
 739      * otherwise, map to VCF from chromosomal coordinates
 740      * of the sequence (if known)
 741      */
 742     GeneLociI seqCoords = seq.getGeneLoci();
 743     if (seqCoords == null)
 744     {
 745       Cache.log.warn(String.format(
 746               "Can't query VCF for %s as chromosome coordinates not known",
 747               seq.getName()));
 748       return null;
 749     }
 750
 751     String species = seqCoords.getSpeciesId();
 752     String chromosome = seqCoords.getChromosomeId();
 753     String seqRef = seqCoords.getAssemblyId();
 754     MapList map = seqCoords.getMapping();
 755
 756     // note this requires the configured species to match that
 757     // returned with the Ensembl sequence; todo: support aliases?
 758     if (!vcfSpecies.equalsIgnoreCase(species))
 759     {
 760       Cache.log.warn("No VCF loaded to " + seq.getName()
 761               + " as species not matched");
 762       return null;
 763     }
 764
 765     if (seqRef.equalsIgnoreCase(vcfAssembly))
 766     {
 767       return new VCFMap(chromosome, map);
 768     }
 769
 770     /*
 771      * VCF data has a different reference assembly to the sequence:
 772      * query Ensembl to map chromosomal coordinates from sequence to VCF
 773      */
 774     List<int[]> toVcfRanges = new ArrayList<>();
 775     List<int[]> fromSequenceRanges = new ArrayList<>();
 776
 777     for (int[] range : map.getToRanges())
 778     {
 779       int[] fromRange = map.locateInFrom(range[0], range[1]);
 780       if (fromRange == null)
 781       {
 782         // corrupted map?!?
 783         continue;
 784       }
 785
 786       int[] newRange = mapReferenceRange(range, chromosome, "human", seqRef,
 787               vcfAssembly);
 788       if (newRange == null)
 789       {
 790         Cache.log.error(
 791                 String.format("Failed to map %s:%s:%s:%d:%d to %s", species,
 792                         chromosome, seqRef, range[0], range[1],
 793                         vcfAssembly));
 794         continue;
 795       }
 796       else
 797       {
 798         toVcfRanges.add(newRange);
 799         fromSequenceRanges.add(fromRange);
 800       }
 801     }
 802
 803     return new VCFMap(chromosome,
 804             new MapList(fromSequenceRanges, toVcfRanges, 1, 1));
 805   }
 806
 807   /**
 808    * If the sequence id matches a contig declared in the VCF file, and the
 809    * sequence length matches the contig length, then returns a 1:1 map of the
 810    * sequence to the contig, else returns null
 811    *
 812    * @param seq
 813    * @return
 814    */
 815   private VCFMap getContigMap(SequenceI seq)
 816   {
 817     String id = seq.getName();
 818     SAMSequenceRecord contig = dictionary.getSequence(id);
 819     if (contig != null)
 820     {
 821       int len = seq.getLength();
 822       if (len == contig.getSequenceLength())
 823       {
 824         MapList map = new MapList(new int[] { 1, len },
 825                 new int[]
 826                 { 1, len }, 1, 1);
 827         return new VCFMap(id, map);
 828       }
 829     }
 830     return null;
 831   }
 832
 833   /**
 834    * Queries the VCF reader for any variants that overlap the mapped chromosome
 835    * ranges of the sequence, and adds as variant features. Returns the number of
 836    * overlapping variants found.
 837    *
 838    * @param seq
 839    * @param map
 840    *          mapping from sequence to VCF coordinates
 841    * @return
 842    */
 843   protected int addVcfVariants(SequenceI seq, VCFMap map)
 844   {
 845     boolean forwardStrand = map.map.isToForwardStrand();
 846
 847     /*
 848      * query the VCF for overlaps of each contiguous chromosomal region
 849      */
 850     int count = 0;
 851
 852     for (int[] range : map.map.getToRanges())
 853     {
 854       int vcfStart = Math.min(range[0], range[1]);
 855       int vcfEnd = Math.max(range[0], range[1]);
 856       CloseableIterator<VariantContext> variants = reader
 857               .query(map.chromosome, vcfStart, vcfEnd);
 858       while (variants.hasNext())
 859       {
 860         VariantContext variant = variants.next();
 861
 862         int[] featureRange = map.map.locateInFrom(variant.getStart(),
 863                 variant.getEnd());
 864
 865         if (featureRange != null)
 866         {
 867           int featureStart = Math.min(featureRange[0], featureRange[1]);
 868           int featureEnd = Math.max(featureRange[0], featureRange[1]);
 869           count += addAlleleFeatures(seq, variant, featureStart, featureEnd,
 870                   forwardStrand);
 871         }
 872       }
 873       variants.close();
 874     }
 875
 876     return count;
 877   }
 878
 879   /**
 880    * A convenience method to get an attribute value for an alternate allele
 881    *
 882    * @param variant
 883    * @param attributeName
 884    * @param alleleIndex
 885    * @return
 886    */
 887   protected String getAttributeValue(VariantContext variant,
 888           String attributeName, int alleleIndex)
 889   {
 890     Object att = variant.getAttribute(attributeName);
 891
 892     if (att instanceof String)
 893     {
 894       return (String) att;
 895     }
 896     else if (att instanceof ArrayList)
 897     {
 898       return ((List<String>) att).get(alleleIndex);
 899     }
 900
 901     return null;
 902   }
 903
 904   /**
 905    * Adds one variant feature for each allele in the VCF variant record, and
 906    * returns the number of features added.
 907    *
 908    * @param seq
 909    * @param variant
 910    * @param featureStart
 911    * @param featureEnd
 912    * @param forwardStrand
 913    * @return
 914    */
 915   protected int addAlleleFeatures(SequenceI seq, VariantContext variant,
 916           int featureStart, int featureEnd, boolean forwardStrand)
 917   {
 918     int added = 0;
 919
 920     /*
 921      * Javadoc says getAlternateAlleles() imposes no order on the list returned
 922      * so we proceed defensively to get them in strict order
 923      */
 924     int altAlleleCount = variant.getAlternateAlleles().size();
 925     for (int i = 0; i < altAlleleCount; i++)
 926     {
 927       added += addAlleleFeature(seq, variant, i, featureStart, featureEnd,
 928               forwardStrand);
 929     }
 930     return added;
 931   }
 932
 933   /**
 934    * Inspects one allele and attempts to add a variant feature for it to the
 935    * sequence. The additional data associated with this allele is extracted to
 936    * store in the feature's key-value map. Answers the number of features added (0
 937    * or 1).
 938    *
 939    * @param seq
 940    * @param variant
 941    * @param altAlleleIndex
 942    *          (0, 1..)
 943    * @param featureStart
 944    * @param featureEnd
 945    * @param forwardStrand
 946    * @return
 947    */
 948   protected int addAlleleFeature(SequenceI seq, VariantContext variant,
 949           int altAlleleIndex, int featureStart, int featureEnd,
 950           boolean forwardStrand)
 951   {
 952     String reference = variant.getReference().getBaseString();
 953     Allele alt = variant.getAlternateAllele(altAlleleIndex);
 954     String allele = alt.getBaseString();
 955
 956     /*
 957      * insertion after a genomic base, if on reverse strand, has to be
 958      * converted to insertion of complement after the preceding position
 959      */
 960     int referenceLength = reference.length();
 961     if (!forwardStrand && allele.length() > referenceLength
 962             && allele.startsWith(reference))
 963     {
 964       featureStart -= referenceLength;
 965       featureEnd = featureStart;
 966       char insertAfter = seq.getCharAt(featureStart - seq.getStart());
 967       reference = Dna.reverseComplement(String.valueOf(insertAfter));
 968       allele = allele.substring(referenceLength) + reference;
 969     }
 970
 971     /*
 972      * build the ref,alt allele description e.g. "G,A", using the base
 973      * complement if the sequence is on the reverse strand
 974      */
 975     StringBuilder sb = new StringBuilder();
 976     sb.append(forwardStrand ? reference : Dna.reverseComplement(reference));
 977     sb.append(COMMA);
 978     sb.append(forwardStrand ? allele : Dna.reverseComplement(allele));
 979     String alleles = sb.toString(); // e.g. G,A
 980
 981     /*
 982      * pick out the consequence data (if any) that is for the current allele
 983      * and feature (transcript) that matches the current sequence
 984      */
 985     String consequence = getConsequenceForAlleleAndFeature(variant, CSQ_FIELD,
 986             altAlleleIndex, csqAlleleFieldIndex,
 987             csqAlleleNumberFieldIndex, seq.getName().toLowerCase(),
 988             csqFeatureFieldIndex);
 989
 990     /*
 991      * pick out the ontology term for the consequence type
 992      */
 993     String type = SequenceOntologyI.SEQUENCE_VARIANT;
 994     if (consequence != null)
 995     {
 996       type = getOntologyTerm(consequence);
 997     }
 998
 999     SequenceFeature sf = new SequenceFeature(type, alleles, featureStart,
1000             featureEnd, FEATURE_GROUP_VCF);
1001     sf.setSource(sourceId);
1002
1003     sf.setValue(Gff3Helper.ALLELES, alleles);
1004
1005     addAlleleProperties(variant, sf, altAlleleIndex, consequence);
1006
1007     seq.addSequenceFeature(sf);
1008
1009     return 1;
1010   }
1011
1012   /**
1013    * Determines the Sequence Ontology term to use for the variant feature type in
1014    * Jalview. The default is 'sequence_variant', but a more specific term is used
1015    * if:
1016    * <ul>
1017    * <li>VEP (or SnpEff) Consequence annotation is included in the VCF</li>
1018    * <li>sequence id can be matched to VEP Feature (or SnpEff Feature_ID)</li>
1019    * </ul>
1020    *
1021    * @param consequence
1022    * @return
1023    * @see http://www.sequenceontology.org/browser/current_svn/term/SO:0001060
1024    */
1025   String getOntologyTerm(String consequence)
1026   {
1027     String type = SequenceOntologyI.SEQUENCE_VARIANT;
1028
1029     /*
1030      * could we associate Consequence data with this allele and feature (transcript)?
1031      * if so, prefer the consequence term from that data
1032      */
1033     if (csqAlleleFieldIndex == -1) // && snpEffAlleleFieldIndex == -1
1034     {
1035       /*
1036        * no Consequence data so we can't refine the ontology term
1037        */
1038       return type;
1039     }
1040
1041     if (consequence != null)
1042     {
1043       String[] csqFields = consequence.split(PIPE_REGEX);
1044       if (csqFields.length > csqConsequenceFieldIndex)
1045       {
1046         type = csqFields[csqConsequenceFieldIndex];
1047       }
1048     }
1049     else
1050     {
1051       // todo the same for SnpEff consequence data matching if wanted
1052     }
1053
1054     /*
1055      * if of the form (e.g.) missense_variant&splice_region_variant,
1056      * just take the first ('most severe') consequence
1057      */
1058     if (type != null)
1059     {
1060       int pos = type.indexOf('&');
1061       if (pos > 0)
1062       {
1063         type = type.substring(0, pos);
1064       }
1065     }
1066     return type;
1067   }
1068
1069   /**
1070    * Returns matched consequence data if it can be found, else null.
1071    * <ul>
1072    * <li>inspects the VCF data for key 'vcfInfoId'</li>
1073    * <li>splits this on comma (to distinct consequences)</li>
1074    * <li>returns the first consequence (if any) where</li>
1075    * <ul>
1076    * <li>the allele matches the altAlleleIndex'th allele of variant</li>
1077    * <li>the feature matches the sequence name (e.g. transcript id)</li>
1078    * </ul>
1079    * </ul>
1080    * If matched, the consequence is returned (as pipe-delimited fields).
1081    *
1082    * @param variant
1083    * @param vcfInfoId
1084    * @param altAlleleIndex
1085    * @param alleleFieldIndex
1086    * @param alleleNumberFieldIndex
1087    * @param seqName
1088    * @param featureFieldIndex
1089    * @return
1090    */
1091   private String getConsequenceForAlleleAndFeature(VariantContext variant,
1092           String vcfInfoId, int altAlleleIndex, int alleleFieldIndex,
1093           int alleleNumberFieldIndex,
1094           String seqName, int featureFieldIndex)
1095   {
1096     if (alleleFieldIndex == -1 || featureFieldIndex == -1)
1097     {
1098       return null;
1099     }
1100     Object value = variant.getAttribute(vcfInfoId);
1101
1102     if (value == null || !(value instanceof List<?>))
1103     {
1104       return null;
1105     }
1106
1107     /*
1108      * inspect each consequence in turn (comma-separated blocks
1109      * extracted by htsjdk)
1110      */
1111     List<String> consequences = (List<String>) value;
1112
1113     for (String consequence : consequences)
1114     {
1115       String[] csqFields = consequence.split(PIPE_REGEX);
1116       if (csqFields.length > featureFieldIndex)
1117       {
1118         String featureIdentifier = csqFields[featureFieldIndex];
1119         if (featureIdentifier.length() > 4
1120                 && seqName.indexOf(featureIdentifier.toLowerCase()) > -1)
1121         {
1122           /*
1123            * feature (transcript) matched - now check for allele match
1124            */
1125           if (matchAllele(variant, altAlleleIndex, csqFields,
1126                   alleleFieldIndex, alleleNumberFieldIndex))
1127           {
1128             return consequence;
1129           }
1130         }
1131       }
1132     }
1133     return null;
1134   }
1135
1136   private boolean matchAllele(VariantContext variant, int altAlleleIndex,
1137           String[] csqFields, int alleleFieldIndex,
1138           int alleleNumberFieldIndex)
1139   {
1140     /*
1141      * if ALLELE_NUM is present, it must match altAlleleIndex
1142      * NB first alternate allele is 1 for ALLELE_NUM, 0 for altAlleleIndex
1143      */
1144     if (alleleNumberFieldIndex > -1)
1145     {
1146       if (csqFields.length <= alleleNumberFieldIndex)
1147       {
1148         return false;
1149       }
1150       String alleleNum = csqFields[alleleNumberFieldIndex];
1151       return String.valueOf(altAlleleIndex + 1).equals(alleleNum);
1152     }
1153
1154     /*
1155      * else consequence allele must match variant allele
1156      */
1157     if (alleleFieldIndex > -1 && csqFields.length > alleleFieldIndex)
1158     {
1159       String csqAllele = csqFields[alleleFieldIndex];
1160       String vcfAllele = variant.getAlternateAllele(altAlleleIndex)
1161               .getBaseString();
1162       return csqAllele.equals(vcfAllele);
1163     }
1164     return false;
1165   }
1166
1167   /**
1168    * Add any allele-specific VCF key-value data to the sequence feature
1169    *
1170    * @param variant
1171    * @param sf
1172    * @param altAlelleIndex
1173    *          (0, 1..)
1174    * @param consequence
1175    *          if not null, the consequence specific to this sequence (transcript
1176    *          feature) and allele
1177    */
1178   protected void addAlleleProperties(VariantContext variant,
1179           SequenceFeature sf, final int altAlelleIndex, String consequence)
1180   {
1181     Map<String, Object> atts = variant.getAttributes();
1182
1183     for (Entry<String, Object> att : atts.entrySet())
1184     {
1185       String key = att.getKey();
1186
1187       /*
1188        * extract Consequence data (if present) that we are able to
1189        * associated with the allele for this variant feature
1190        */
1191       if (CSQ_FIELD.equals(key))
1192       {
1193         addConsequences(variant, sf, consequence);
1194         continue;
1195       }
1196
1197       /*
1198        * filter out fields we don't want to capture
1199        */
1200       if (!vcfFieldsOfInterest.contains(key))
1201       {
1202         continue;
1203       }
1204
1205       /*
1206        * we extract values for other data which are allele-specific;
1207        * these may be per alternate allele (INFO[key].Number = 'A')
1208        * or per allele including reference (INFO[key].Number = 'R')
1209        */
1210       VCFInfoHeaderLine infoHeader = header.getInfoHeaderLine(key);
1211       if (infoHeader == null)
1212       {
1213         /*
1214          * can't be sure what data belongs to this allele, so
1215          * play safe and don't take any
1216          */
1217         continue;
1218       }
1219
1220       VCFHeaderLineCount number = infoHeader.getCountType();
1221       int index = altAlelleIndex;
1222       if (number == VCFHeaderLineCount.R)
1223       {
1224         /*
1225          * one value per allele including reference, so bump index
1226          * e.g. the 3rd value is for the  2nd alternate allele
1227          */
1228         index++;
1229       }
1230       else if (number != VCFHeaderLineCount.A)
1231       {
1232         /*
1233          * don't save other values as not allele-related
1234          */
1235         continue;
1236       }
1237
1238       /*
1239        * take the index'th value
1240        */
1241       String value = getAttributeValue(variant, key, index);
1242       if (value != null)
1243       {
1244         value = decodeSpecialCharacters(value);
1245         sf.setValue(key, value);
1246       }
1247     }
1248   }
1249
1250   /**
1251    * Decodes colon, semicolon, equals sign, percent sign, comma to their decoded
1252    * form. The VCF specification (para 1.2) requires these to be encoded where not
1253    * used with their special meaning in the VCF syntax. Note that general URL
1254    * decoding should not be applied, since this would incorrectly decode (for
1255    * example) a '+' sign.
1256    *
1257    * @param value
1258    * @return
1259    */
1260   protected static String decodeSpecialCharacters(String value)
1261   {
1262     /*
1263      * avoid regex compilation if it is not needed!
1264      */
1265     if (!value.contains(ENCODED_COLON) && !value.contains(ENCODED_SEMICOLON)
1266             && !value.contains(ENCODED_EQUALS)
1267             && !value.contains(ENCODED_PERCENT)
1268             && !value.contains(ENCODED_COMMA))
1269     {
1270       return value;
1271     }
1272
1273     value = value.replace(ENCODED_COLON, ":")
1274             .replace(ENCODED_SEMICOLON, ";").replace(ENCODED_EQUALS, "=")
1275             .replace(ENCODED_PERCENT, "%").replace(ENCODED_COMMA, ",");
1276     return value;
1277   }
1278
1279   /**
1280    * Inspects CSQ data blocks (consequences) and adds attributes on the sequence
1281    * feature.
1282    * <p>
1283    * If <code>myConsequence</code> is not null, then this is the specific
1284    * consequence data (pipe-delimited fields) that is for the current allele and
1285    * transcript (sequence) being processed)
1286    *
1287    * @param variant
1288    * @param sf
1289    * @param myConsequence
1290    */
1291   protected void addConsequences(VariantContext variant, SequenceFeature sf,
1292           String myConsequence)
1293   {
1294     Object value = variant.getAttribute(CSQ_FIELD);
1295
1296     if (value == null || !(value instanceof List<?>))
1297     {
1298       return;
1299     }
1300
1301     List<String> consequences = (List<String>) value;
1302
1303     /*
1304      * inspect CSQ consequences; restrict to the consequence
1305      * associated with the current transcript (Feature)
1306      */
1307     Map<String, String> csqValues = new HashMap<>();
1308
1309     for (String consequence : consequences)
1310     {
1311       if (myConsequence == null || myConsequence.equals(consequence))
1312       {
1313         String[] csqFields = consequence.split(PIPE_REGEX);
1314
1315         /*
1316          * inspect individual fields of this consequence, copying non-null
1317          * values which are 'fields of interest'
1318          */
1319         int i = 0;
1320         for (String field : csqFields)
1321         {
1322           if (field != null && field.length() > 0)
1323           {
1324             String id = vepFieldsOfInterest.get(i);
1325             if (id != null)
1326             {
1327               /*
1328                * VCF spec requires encoding of special characters e.g. '='
1329                * so decode them here before storing
1330                */
1331               try
1332               {
1333                 field = URLDecoder.decode(field, UTF_8);
1334               } catch (UnsupportedEncodingException e)
1335               {
1336               }
1337               csqValues.put(id, field);
1338             }
1339           }
1340           i++;
1341         }
1342       }
1343     }
1344
1345     if (!csqValues.isEmpty())
1346     {
1347       sf.setValue(CSQ_FIELD, csqValues);
1348     }
1349   }
1350
1351   /**
1352    * A convenience method to complement a dna base and return the string value
1353    * of its complement
1354    *
1355    * @param reference
1356    * @return
1357    */
1358   protected String complement(byte[] reference)
1359   {
1360     return String.valueOf(Dna.getComplement((char) reference[0]));
1361   }
1362
1363   /**
1364    * Determines the location of the query range (chromosome positions) in a
1365    * different reference assembly.
1366    * <p>
1367    * If the range is just a subregion of one for which we already have a mapping
1368    * (for example, an exon sub-region of a gene), then the mapping is just
1369    * computed arithmetically.
1370    * <p>
1371    * Otherwise, calls the Ensembl REST service that maps from one assembly
1372    * reference's coordinates to another's
1373    *
1374    * @param queryRange
1375    *          start-end chromosomal range in 'fromRef' coordinates
1376    * @param chromosome
1377    * @param species
1378    * @param fromRef
1379    *          assembly reference for the query coordinates
1380    * @param toRef
1381    *          assembly reference we wish to translate to
1382    * @return the start-end range in 'toRef' coordinates
1383    */
1384   protected int[] mapReferenceRange(int[] queryRange, String chromosome,
1385           String species, String fromRef, String toRef)
1386   {
1387     /*
1388      * first try shorcut of computing the mapping as a subregion of one
1389      * we already have (e.g. for an exon, if we have the gene mapping)
1390      */
1391     int[] mappedRange = findSubsumedRangeMapping(queryRange, chromosome,
1392             species, fromRef, toRef);
1393     if (mappedRange != null)
1394     {
1395       return mappedRange;
1396     }
1397
1398     /*
1399      * call (e.g.) http://rest.ensembl.org/map/human/GRCh38/17:45051610..45109016:1/GRCh37
1400      */
1401     EnsemblMap mapper = new EnsemblMap();
1402     int[] mapping = mapper.getAssemblyMapping(species, chromosome, fromRef,
1403             toRef, queryRange);
1404
1405     if (mapping == null)
1406     {
1407       // mapping service failure
1408       return null;
1409     }
1410
1411     /*
1412      * save mapping for possible future re-use
1413      */
1414     String key = makeRangesKey(chromosome, species, fromRef, toRef);
1415     if (!assemblyMappings.containsKey(key))
1416     {
1417       assemblyMappings.put(key, new HashMap<int[], int[]>());
1418     }
1419
1420     assemblyMappings.get(key).put(queryRange, mapping);
1421
1422     return mapping;
1423   }
1424
1425   /**
1426    * If we already have a 1:1 contiguous mapping which subsumes the given query
1427    * range, this method just calculates and returns the subset of that mapping,
1428    * else it returns null. In practical terms, if a gene has a contiguous
1429    * mapping between (for example) GRCh37 and GRCh38, then we assume that its
1430    * subsidiary exons occupy unchanged relative positions, and just compute
1431    * these as offsets, rather than do another lookup of the mapping.
1432    * <p>
1433    * If in future these assumptions prove invalid (e.g. for bacterial dna?!),
1434    * simply remove this method or let it always return null.
1435    * <p>
1436    * Warning: many rapid calls to the /map service map result in a 429 overload
1437    * error response
1438    *
1439    * @param queryRange
1440    * @param chromosome
1441    * @param species
1442    * @param fromRef
1443    * @param toRef
1444    * @return
1445    */
1446   protected int[] findSubsumedRangeMapping(int[] queryRange, String chromosome,
1447           String species, String fromRef, String toRef)
1448   {
1449     String key = makeRangesKey(chromosome, species, fromRef, toRef);
1450     if (assemblyMappings.containsKey(key))
1451     {
1452       Map<int[], int[]> mappedRanges = assemblyMappings.get(key);
1453       for (Entry<int[], int[]> mappedRange : mappedRanges.entrySet())
1454       {
1455         int[] fromRange = mappedRange.getKey();
1456         int[] toRange = mappedRange.getValue();
1457         if (fromRange[1] - fromRange[0] == toRange[1] - toRange[0])
1458         {
1459           /*
1460            * mapping is 1:1 in length, so we trust it to have no discontinuities
1461            */
1462           if (MappingUtils.rangeContains(fromRange, queryRange))
1463           {
1464             /*
1465              * fromRange subsumes our query range
1466              */
1467             int offset = queryRange[0] - fromRange[0];
1468             int mappedRangeFrom = toRange[0] + offset;
1469             int mappedRangeTo = mappedRangeFrom + (queryRange[1] - queryRange[0]);
1470             return new int[] { mappedRangeFrom, mappedRangeTo };
1471           }
1472         }
1473       }
1474     }
1475     return null;
1476   }
1477
1478   /**
1479    * Transfers the sequence feature to the target sequence, locating its start
1480    * and end range based on the mapping. Features which do not overlap the
1481    * target sequence are ignored.
1482    *
1483    * @param sf
1484    * @param targetSequence
1485    * @param mapping
1486    *          mapping from the feature's coordinates to the target sequence
1487    */
1488   protected void transferFeature(SequenceFeature sf,
1489           SequenceI targetSequence, MapList mapping)
1490   {
1491     int[] mappedRange = mapping.locateInTo(sf.getBegin(), sf.getEnd());
1492
1493     if (mappedRange != null)
1494     {
1495       String group = sf.getFeatureGroup();
1496       int newBegin = Math.min(mappedRange[0], mappedRange[1]);
1497       int newEnd = Math.max(mappedRange[0], mappedRange[1]);
1498       SequenceFeature copy = new SequenceFeature(sf, newBegin, newEnd,
1499               group, sf.getScore());
1500       targetSequence.addSequenceFeature(copy);
1501     }
1502   }
1503
1504   /**
1505    * Formats a ranges map lookup key
1506    *
1507    * @param chromosome
1508    * @param species
1509    * @param fromRef
1510    * @param toRef
1511    * @return
1512    */
1513   protected static String makeRangesKey(String chromosome, String species,
1514           String fromRef, String toRef)
1515   {
1516     return species + EXCL + chromosome + EXCL + fromRef + EXCL
1517             + toRef;
1518   }
1519 }