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