src/jalview/io/vcf/VCFLoader.java

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