src/jalview/io/vcf/VCFLoader.java

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