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