src/jalview/io/vcf/VCFLoader.java

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