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