src/jalview/io/vcf/VCFLoader.java

   1 package jalview.io.vcf;
   2
   3 import htsjdk.samtools.util.CloseableIterator;
   4 import htsjdk.variant.variantcontext.Allele;
   5 import htsjdk.variant.variantcontext.VariantContext;
   6 import htsjdk.variant.vcf.VCFHeader;
   7 import htsjdk.variant.vcf.VCFHeaderLine;
   8 import htsjdk.variant.vcf.VCFHeaderLineCount;
   9 import htsjdk.variant.vcf.VCFHeaderLineType;
  10 import htsjdk.variant.vcf.VCFInfoHeaderLine;
  11
  12 import jalview.analysis.AlignmentUtils;
  13 import jalview.analysis.Dna;
  14 import jalview.api.AlignViewControllerGuiI;
  15 import jalview.datamodel.AlignmentI;
  16 import jalview.datamodel.DBRefEntry;
  17 import jalview.datamodel.GeneLociI;
  18 import jalview.datamodel.Mapping;
  19 import jalview.datamodel.SequenceFeature;
  20 import jalview.datamodel.SequenceI;
  21 import jalview.datamodel.features.FeatureAttributeType;
  22 import jalview.datamodel.features.FeatureSource;
  23 import jalview.datamodel.features.FeatureSources;
  24 import jalview.ext.ensembl.EnsemblMap;
  25 import jalview.ext.htsjdk.VCFReader;
  26 import jalview.io.gff.Gff3Helper;
  27 import jalview.io.gff.SequenceOntologyI;
  28 import jalview.util.MapList;
  29 import jalview.util.MappingUtils;
  30 import jalview.util.MessageManager;
  31
  32 import java.io.IOException;
  33 import java.util.ArrayList;
  34 import java.util.HashMap;
  35 import java.util.List;
  36 import java.util.Map;
  37 import java.util.Map.Entry;
  38
  39 /**
  40  * A class to read VCF data (using the htsjdk) and add variants as sequence
  41  * features on dna and any related protein product sequences
  42  *
  43  * @author gmcarstairs
  44  */
  45 public class VCFLoader
  46 {
  47   /*
  48    * keys to fields of VEP CSQ consequence data
  49    * see https://www.ensembl.org/info/docs/tools/vep/vep_formats.html
  50    */
  51   private static final String ALLELE_KEY = "Allele";
  52
  53   private static final String ALLELE_NUM_KEY = "ALLELE_NUM"; // 0 (ref), 1...
  54   private static final String FEATURE_KEY = "Feature"; // Ensembl stable id
  55
  56   /*
  57    * what comes before column headings in CSQ Description field
  58    */
  59   private static final String FORMAT = "Format: ";
  60
  61   /*
  62    * default VCF INFO key for VEP consequence data
  63    * NB this can be overridden running VEP with --vcf_info_field
  64    * - we don't handle this case (require identifier to be CSQ)
  65    */
  66   private static final String CSQ_FIELD = "CSQ";
  67
  68   /*
  69    * separator for fields in consequence data is '|'
  70    */
  71   private static final String PIPE_REGEX = "\\|";
  72
  73   /*
  74    * key for Allele Frequency output by VEP
  75    * see http://www.ensembl.org/info/docs/tools/vep/vep_formats.html
  76    */
  77   private static final String ALLELE_FREQUENCY_KEY = "AF";
  78
  79   /*
  80    * delimiter that separates multiple consequence data blocks
  81    */
  82   private static final String COMMA = ",";
  83
  84   /*
  85    * the feature group assigned to a VCF variant in Jalview
  86    */
  87   private static final String FEATURE_GROUP_VCF = "VCF";
  88
  89   /*
  90    * internal delimiter used to build keys for assemblyMappings
  91    *
  92    */
  93   private static final String EXCL = "!";
  94
  95   /*
  96    * the alignment we are associating VCF data with
  97    */
  98   private AlignmentI al;
  99
 100   /*
 101    * mappings between VCF and sequence reference assembly regions, as
 102    * key = "species!chromosome!fromAssembly!toAssembly
 103    * value = Map{fromRange, toRange}
 104    */
 105   private Map<String, Map<int[], int[]>> assemblyMappings;
 106
 107   /*
 108    * holds details of the VCF header lines (metadata)
 109    */
 110   private VCFHeader header;
 111
 112   /*
 113    * the position (0...) of field in each block of
 114    * CSQ (consequence) data (if declared in the VCF INFO header for CSQ)
 115    * see http://www.ensembl.org/info/docs/tools/vep/vep_formats.html
 116    */
 117   private int csqAlleleFieldIndex = -1;
 118   private int csqAlleleNumberFieldIndex = -1;
 119   private int csqFeatureFieldIndex = -1;
 120
 121   /*
 122    * a unique identifier under which to save metadata about feature
 123    * attributes (selected INFO field data)
 124    */
 125   private String sourceId;
 126
 127   List<String> vcfFieldsOfInterest;
 128
 129   List<String> vepFieldsOfInterest;
 130
 131   /**
 132    * Constructor given an alignment context
 133    *
 134    * @param alignment
 135    */
 136   public VCFLoader(AlignmentI alignment)
 137   {
 138     al = alignment;
 139
 140     // map of species!chromosome!fromAssembly!toAssembly to {fromRange, toRange}
 141     assemblyMappings = new HashMap<String, Map<int[], int[]>>();
 142   }
 143
 144   /**
 145    * Starts a new thread to query and load VCF variant data on to the alignment
 146    * <p>
 147    * This method is not thread safe - concurrent threads should use separate
 148    * instances of this class.
 149    *
 150    * @param filePath
 151    * @param gui
 152    */
 153   public void loadVCF(final String filePath,
 154           final AlignViewControllerGuiI gui)
 155   {
 156     if (gui != null)
 157     {
 158       gui.setStatus(MessageManager.getString("label.searching_vcf"));
 159     }
 160
 161     new Thread()
 162     {
 163
 164       @Override
 165       public void run()
 166       {
 167         VCFLoader.this.doLoad(filePath, gui);
 168       }
 169
 170     }.start();
 171   }
 172
 173   /**
 174    * Loads VCF on to an alignment - provided it can be related to one or more
 175    * sequence's chromosomal coordinates
 176    *
 177    * @param filePath
 178    * @param gui
 179    *          optional callback handler for messages
 180    */
 181   protected void doLoad(String filePath, AlignViewControllerGuiI gui)
 182   {
 183     VCFReader reader = null;
 184     try
 185     {
 186       // long start = System.currentTimeMillis();
 187       reader = new VCFReader(filePath);
 188
 189       header = reader.getFileHeader();
 190
 191       sourceId = filePath;
 192
 193       saveMetadata(sourceId);
 194
 195       /*
 196        * get offset of CSQ ALLELE_NUM and Feature if declared
 197        */
 198       locateCsqFields();
 199
 200       VCFHeaderLine ref = header
 201               .getOtherHeaderLine(VCFHeader.REFERENCE_KEY);
 202       String vcfAssembly = ref.getValue();
 203
 204       int varCount = 0;
 205       int seqCount = 0;
 206
 207       /*
 208        * query for VCF overlapping each sequence in turn
 209        */
 210       for (SequenceI seq : al.getSequences())
 211       {
 212         int added = loadSequenceVCF(seq, reader, vcfAssembly);
 213         if (added > 0)
 214         {
 215           seqCount++;
 216           varCount += added;
 217           transferAddedFeatures(seq);
 218         }
 219       }
 220       if (gui != null)
 221       {
 222         // long elapsed = System.currentTimeMillis() - start;
 223         String msg = MessageManager.formatMessage("label.added_vcf",
 224                 varCount, seqCount);
 225         gui.setStatus(msg);
 226         if (gui.getFeatureSettingsUI() != null)
 227         {
 228           gui.getFeatureSettingsUI().discoverAllFeatureData();
 229         }
 230       }
 231     } catch (Throwable e)
 232     {
 233       System.err.println("Error processing VCF: " + e.getMessage());
 234       e.printStackTrace();
 235       if (gui != null)
 236       {
 237         gui.setStatus("Error occurred - see console for details");
 238       }
 239     } finally
 240     {
 241       if (reader != null)
 242       {
 243         try
 244         {
 245           reader.close();
 246         } catch (IOException e)
 247         {
 248           // ignore
 249         }
 250       }
 251     }
 252   }
 253
 254   /**
 255    * Reads metadata (such as INFO field descriptions and datatypes) and saves
 256    * them for future reference
 257    *
 258    * @param theSourceId
 259    */
 260   void saveMetadata(String theSourceId)
 261   {
 262     vcfFieldsOfInterest = new ArrayList<>();
 263
 264     FeatureSource metadata = new FeatureSource(theSourceId);
 265
 266     for (VCFInfoHeaderLine info : header.getInfoHeaderLines())
 267     {
 268       String attributeId = info.getID();
 269       String desc = info.getDescription();
 270       VCFHeaderLineType type = info.getType();
 271       FeatureAttributeType attType = null;
 272       switch (type)
 273       {
 274       case Character:
 275         attType = FeatureAttributeType.Character;
 276         break;
 277       case Flag:
 278         attType = FeatureAttributeType.Flag;
 279         break;
 280       case Float:
 281         attType = FeatureAttributeType.Float;
 282         break;
 283       case Integer:
 284         attType = FeatureAttributeType.Integer;
 285         break;
 286       case String:
 287         attType = FeatureAttributeType.String;
 288         break;
 289       }
 290       metadata.setAttributeName(attributeId, desc);
 291       metadata.setAttributeType(attributeId, attType);
 292
 293       if (isVcfFieldWanted(attributeId))
 294       {
 295         vcfFieldsOfInterest.add(attributeId);
 296       }
 297     }
 298
 299     FeatureSources.getInstance().addSource(theSourceId, metadata);
 300   }
 301
 302   /**
 303    * Answers true if the VCF id is one we wish to capture in Jalview, else false
 304    *
 305    * @param id
 306    * @return
 307    */
 308   private boolean isVcfFieldWanted(String id)
 309   {
 310     // TODO option to match patterns in a Preferences entry?
 311     return true;
 312   }
 313
 314   /**
 315    * Answers true if the VEP (CSQ) id is one we wish to capture in Jalview, else
 316    * false
 317    *
 318    * @param id
 319    * @return
 320    */
 321   private boolean isVepFieldWanted(String id)
 322   {
 323     // TODO option to match patterns in a Preferences entry?
 324     return true;
 325   }
 326
 327   /**
 328    * Records the position of selected fields defined in the CSQ INFO header (if
 329    * there is one). CSQ fields are declared in the CSQ INFO Description e.g.
 330    * <p>
 331    * Description="Consequence ...from ... VEP. Format: Allele|Consequence|...
 332    */
 333   protected void locateCsqFields()
 334   {
 335     vepFieldsOfInterest = new ArrayList<>();
 336
 337     VCFInfoHeaderLine csqInfo = header.getInfoHeaderLine(CSQ_FIELD);
 338     if (csqInfo == null)
 339     {
 340       return;
 341     }
 342
 343     String desc = csqInfo.getDescription();
 344     int formatPos = desc.indexOf(FORMAT);
 345     if (formatPos == -1)
 346     {
 347       System.err.println("Parse error, failed to find " + FORMAT
 348               + " in " + desc);
 349       return;
 350     }
 351     desc = desc.substring(formatPos + FORMAT.length());
 352
 353     if (desc != null)
 354     {
 355       String[] format = desc.split(PIPE_REGEX);
 356       int index = 0;
 357       for (String field : format)
 358       {
 359         if (ALLELE_NUM_KEY.equals(field))
 360         {
 361           csqAlleleNumberFieldIndex = index;
 362         }
 363         if (ALLELE_KEY.equals(field))
 364         {
 365           csqAlleleFieldIndex = index;
 366         }
 367         if (FEATURE_KEY.equals(field))
 368         {
 369           csqFeatureFieldIndex = index;
 370         }
 371
 372         if (isVepFieldWanted(field))
 373         {
 374           vepFieldsOfInterest.add(field);
 375         }
 376
 377         index++;
 378       }
 379     }
 380   }
 381
 382   /**
 383    * Transfers VCF features to sequences to which this sequence has a mapping.
 384    * If the mapping is 3:1, computes peptide variants from nucleotide variants.
 385    *
 386    * @param seq
 387    */
 388   protected void transferAddedFeatures(SequenceI seq)
 389   {
 390     DBRefEntry[] dbrefs = seq.getDBRefs();
 391     if (dbrefs == null)
 392     {
 393       return;
 394     }
 395     for (DBRefEntry dbref : dbrefs)
 396     {
 397       Mapping mapping = dbref.getMap();
 398       if (mapping == null || mapping.getTo() == null)
 399       {
 400         continue;
 401       }
 402
 403       SequenceI mapTo = mapping.getTo();
 404       MapList map = mapping.getMap();
 405       if (map.getFromRatio() == 3)
 406       {
 407         /*
 408          * dna-to-peptide product mapping
 409          */
 410         AlignmentUtils.computeProteinFeatures(seq, mapTo, map);
 411       }
 412       else
 413       {
 414         /*
 415          * nucleotide-to-nucleotide mapping e.g. transcript to CDS
 416          */
 417         List<SequenceFeature> features = seq.getFeatures()
 418                 .getPositionalFeatures(SequenceOntologyI.SEQUENCE_VARIANT);
 419         for (SequenceFeature sf : features)
 420         {
 421           if (FEATURE_GROUP_VCF.equals(sf.getFeatureGroup()))
 422           {
 423             transferFeature(sf, mapTo, map);
 424           }
 425         }
 426       }
 427     }
 428   }
 429
 430   /**
 431    * Tries to add overlapping variants read from a VCF file to the given
 432    * sequence, and returns the number of variant features added. Note that this
 433    * requires the sequence to hold information as to its species, chromosomal
 434    * positions and reference assembly, in order to be able to map the VCF
 435    * variants to the sequence (or not)
 436    *
 437    * @param seq
 438    * @param reader
 439    * @param vcfAssembly
 440    * @return
 441    */
 442   protected int loadSequenceVCF(SequenceI seq, VCFReader reader,
 443           String vcfAssembly)
 444   {
 445     int count = 0;
 446     GeneLociI seqCoords = seq.getGeneLoci();
 447     if (seqCoords == null)
 448     {
 449       System.out.println(String.format(
 450               "Can't query VCF for %s as chromosome coordinates not known",
 451               seq.getName()));
 452       return 0;
 453     }
 454
 455     if (!vcfSpeciesMatchesSequence(vcfAssembly, seqCoords.getSpeciesId()))
 456     {
 457       return 0;
 458     }
 459
 460     List<int[]> seqChromosomalContigs = seqCoords.getMap().getToRanges();
 461     for (int[] range : seqChromosomalContigs)
 462     {
 463       count += addVcfVariants(seq, reader, range, vcfAssembly);
 464     }
 465
 466     return count;
 467   }
 468
 469   /**
 470    * Answers true if the species inferred from the VCF reference identifier
 471    * matches that for the sequence
 472    *
 473    * @param vcfAssembly
 474    * @param speciesId
 475    * @return
 476    */
 477   boolean vcfSpeciesMatchesSequence(String vcfAssembly, String speciesId)
 478   {
 479     // PROBLEM 1
 480     // there are many aliases for species - how to equate one with another?
 481     // PROBLEM 2
 482     // VCF ##reference header is an unstructured URI - how to extract species?
 483     // perhaps check if ref includes any (Ensembl) alias of speciesId??
 484     // TODO ask the user to confirm this??
 485
 486     if (vcfAssembly.contains("Homo_sapiens") // gnomAD exome data example
 487             && "HOMO_SAPIENS".equals(speciesId)) // Ensembl species id
 488     {
 489       return true;
 490     }
 491
 492     if (vcfAssembly.contains("c_elegans") // VEP VCF response example
 493             && "CAENORHABDITIS_ELEGANS".equals(speciesId)) // Ensembl
 494     {
 495       return true;
 496     }
 497
 498     // this is not a sustainable solution...
 499
 500     return false;
 501   }
 502
 503   /**
 504    * Queries the VCF reader for any variants that overlap the given chromosome
 505    * region of the sequence, and adds as variant features. Returns the number of
 506    * overlapping variants found.
 507    *
 508    * @param seq
 509    * @param reader
 510    * @param range
 511    *          start-end range of a sequence region in its chromosomal
 512    *          coordinates
 513    * @param vcfAssembly
 514    *          the '##reference' identifier for the VCF reference assembly
 515    * @return
 516    */
 517   protected int addVcfVariants(SequenceI seq, VCFReader reader,
 518           int[] range, String vcfAssembly)
 519   {
 520     GeneLociI seqCoords = seq.getGeneLoci();
 521
 522     String chromosome = seqCoords.getChromosomeId();
 523     String seqRef = seqCoords.getAssemblyId();
 524     String species = seqCoords.getSpeciesId();
 525
 526     /*
 527      * map chromosomal coordinates from sequence to VCF if the VCF
 528      * data has a different reference assembly to the sequence
 529      */
 530     // TODO generalise for non-human species
 531     // - or get the user to choose in a dialog
 532
 533     int offset = 0;
 534     if ("GRCh38".equalsIgnoreCase(seqRef) // Ensembl
 535             && vcfAssembly.contains("Homo_sapiens_assembly19")) // gnomAD
 536     {
 537       String toRef = "GRCh37";
 538       int[] newRange = mapReferenceRange(range, chromosome, "human",
 539               seqRef, toRef);
 540       if (newRange == null)
 541       {
 542         System.err.println(String.format(
 543                 "Failed to map %s:%s:%s:%d:%d to %s", species, chromosome,
 544                 seqRef, range[0], range[1], toRef));
 545         return 0;
 546       }
 547       offset = newRange[0] - range[0];
 548       range = newRange;
 549     }
 550
 551     boolean forwardStrand = range[0] <= range[1];
 552
 553     /*
 554      * query the VCF for overlaps
 555      * (convert a reverse strand range to forwards)
 556      */
 557     int count = 0;
 558     MapList mapping = seqCoords.getMap();
 559
 560     int fromLocus = Math.min(range[0], range[1]);
 561     int toLocus = Math.max(range[0], range[1]);
 562     CloseableIterator<VariantContext> variants = reader.query(chromosome,
 563             fromLocus, toLocus);
 564     while (variants.hasNext())
 565     {
 566       /*
 567        * get variant location in sequence chromosomal coordinates
 568        */
 569       VariantContext variant = variants.next();
 570
 571       int start = variant.getStart() - offset;
 572       int end = variant.getEnd() - offset;
 573
 574       /*
 575        * convert chromosomal location to sequence coordinates
 576        * - may be reverse strand (convert to forward for sequence feature)
 577        * - null if a partially overlapping feature
 578        */
 579       int[] seqLocation = mapping.locateInFrom(start, end);
 580       if (seqLocation != null)
 581       {
 582         int featureStart = Math.min(seqLocation[0], seqLocation[1]);
 583         int featureEnd = Math.max(seqLocation[0], seqLocation[1]);
 584         count += addAlleleFeatures(seq, variant, featureStart, featureEnd,
 585                 forwardStrand);
 586       }
 587     }
 588
 589     variants.close();
 590
 591     return count;
 592   }
 593
 594   /**
 595    * A convenience method to get the AF value for the given alternate allele
 596    * index
 597    *
 598    * @param variant
 599    * @param alleleIndex
 600    * @return
 601    */
 602   protected float getAlleleFrequency(VariantContext variant, int alleleIndex)
 603   {
 604     float score = 0f;
 605     String attributeValue = getAttributeValue(variant,
 606             ALLELE_FREQUENCY_KEY, alleleIndex);
 607     if (attributeValue != null)
 608     {
 609       try
 610       {
 611         score = Float.parseFloat(attributeValue);
 612       } catch (NumberFormatException e)
 613       {
 614         // leave as 0
 615       }
 616     }
 617
 618     return score;
 619   }
 620
 621   /**
 622    * A convenience method to get an attribute value for an alternate allele
 623    *
 624    * @param variant
 625    * @param attributeName
 626    * @param alleleIndex
 627    * @return
 628    */
 629   protected String getAttributeValue(VariantContext variant,
 630           String attributeName, int alleleIndex)
 631   {
 632     Object att = variant.getAttribute(attributeName);
 633
 634     if (att instanceof String)
 635     {
 636       return (String) att;
 637     }
 638     else if (att instanceof ArrayList)
 639     {
 640       return ((List<String>) att).get(alleleIndex);
 641     }
 642
 643     return null;
 644   }
 645
 646   /**
 647    * Adds one variant feature for each allele in the VCF variant record, and
 648    * returns the number of features added.
 649    *
 650    * @param seq
 651    * @param variant
 652    * @param featureStart
 653    * @param featureEnd
 654    * @param forwardStrand
 655    * @return
 656    */
 657   protected int addAlleleFeatures(SequenceI seq, VariantContext variant,
 658           int featureStart, int featureEnd, boolean forwardStrand)
 659   {
 660     int added = 0;
 661
 662     /*
 663      * Javadoc says getAlternateAlleles() imposes no order on the list returned
 664      * so we proceed defensively to get them in strict order
 665      */
 666     int altAlleleCount = variant.getAlternateAlleles().size();
 667     for (int i = 0; i < altAlleleCount; i++)
 668     {
 669       added += addAlleleFeature(seq, variant, i, featureStart, featureEnd,
 670               forwardStrand);
 671     }
 672     return added;
 673   }
 674
 675   /**
 676    * Inspects one allele and attempts to add a variant feature for it to the
 677    * sequence. We extract as much as possible of the additional data associated
 678    * with this allele to store in the feature's key-value map. Answers the
 679    * number of features added (0 or 1).
 680    *
 681    * @param seq
 682    * @param variant
 683    * @param altAlleleIndex
 684    *          (0, 1..)
 685    * @param featureStart
 686    * @param featureEnd
 687    * @param forwardStrand
 688    * @return
 689    */
 690   protected int addAlleleFeature(SequenceI seq, VariantContext variant,
 691           int altAlleleIndex, int featureStart, int featureEnd,
 692           boolean forwardStrand)
 693   {
 694     String reference = variant.getReference().getBaseString();
 695     Allele alt = variant.getAlternateAllele(altAlleleIndex);
 696     String allele = alt.getBaseString();
 697
 698     /*
 699      * build the ref,alt allele description e.g. "G,A", using the base
 700      * complement if the sequence is on the reverse strand
 701      */
 702     // TODO check how structural variants are shown on reverse strand
 703     StringBuilder sb = new StringBuilder();
 704     sb.append(forwardStrand ? reference : Dna.reverseComplement(reference));
 705     sb.append(COMMA);
 706     sb.append(forwardStrand ? allele : Dna.reverseComplement(allele));
 707     String alleles = sb.toString(); // e.g. G,A
 708
 709     String type = SequenceOntologyI.SEQUENCE_VARIANT;
 710     float score = getAlleleFrequency(variant, altAlleleIndex);
 711
 712     SequenceFeature sf = new SequenceFeature(type, alleles, featureStart,
 713             featureEnd, score, FEATURE_GROUP_VCF);
 714     sf.setSource(sourceId);
 715
 716     sf.setValue(Gff3Helper.ALLELES, alleles);
 717
 718     addAlleleProperties(variant, seq, sf, altAlleleIndex);
 719
 720     seq.addSequenceFeature(sf);
 721
 722     return 1;
 723   }
 724
 725   /**
 726    * Add any allele-specific VCF key-value data to the sequence feature
 727    *
 728    * @param variant
 729    * @param seq
 730    * @param sf
 731    * @param altAlelleIndex
 732    *          (0, 1..)
 733    */
 734   protected void addAlleleProperties(VariantContext variant, SequenceI seq,
 735           SequenceFeature sf, final int altAlelleIndex)
 736   {
 737     Map<String, Object> atts = variant.getAttributes();
 738
 739     for (Entry<String, Object> att : atts.entrySet())
 740     {
 741       String key = att.getKey();
 742
 743       /*
 744        * extract Consequence data (if present) that we are able to
 745        * associated with the allele for this variant feature
 746        */
 747       if (CSQ_FIELD.equals(key))
 748       {
 749         addConsequences(variant, seq, sf, altAlelleIndex);
 750         continue;
 751       }
 752
 753       /*
 754        * filter out fields we don't want to capture
 755        */
 756       if (!vcfFieldsOfInterest.contains(key))
 757       {
 758         continue;
 759       }
 760
 761       /*
 762        * we extract values for other data which are allele-specific;
 763        * these may be per alternate allele (INFO[key].Number = 'A')
 764        * or per allele including reference (INFO[key].Number = 'R')
 765        */
 766       VCFInfoHeaderLine infoHeader = header.getInfoHeaderLine(key);
 767       if (infoHeader == null)
 768       {
 769         /*
 770          * can't be sure what data belongs to this allele, so
 771          * play safe and don't take any
 772          */
 773         continue;
 774       }
 775
 776       VCFHeaderLineCount number = infoHeader.getCountType();
 777       int index = altAlelleIndex;
 778       if (number == VCFHeaderLineCount.R)
 779       {
 780         /*
 781          * one value per allele including reference, so bump index
 782          * e.g. the 3rd value is for the  2nd alternate allele
 783          */
 784         index++;
 785       }
 786       else if (number != VCFHeaderLineCount.A)
 787       {
 788         /*
 789          * don't save other values as not allele-related
 790          */
 791         continue;
 792       }
 793
 794       /*
 795        * take the index'th value
 796        */
 797       String value = getAttributeValue(variant, key, index);
 798       if (value != null)
 799       {
 800         sf.setValue(key, value);
 801       }
 802     }
 803   }
 804
 805   /**
 806    * Inspects CSQ data blocks (consequences) and adds attributes on the sequence
 807    * feature for the current allele (and transcript if applicable)
 808    * <p>
 809    * Allele matching: if field ALLELE_NUM is present, it must match
 810    * altAlleleIndex. If not present, then field Allele value must match the VCF
 811    * Allele.
 812    * <p>
 813    * Transcript matching: if sequence name can be identified to at least one of
 814    * the consequences' Feature values, then select only consequences that match
 815    * the value (i.e. consequences for the current transcript sequence). If not,
 816    * take all consequences (this is the case when adding features to the gene
 817    * sequence).
 818    *
 819    * @param variant
 820    * @param seq
 821    * @param sf
 822    * @param altAlelleIndex
 823    *          (0, 1..)
 824    */
 825   protected void addConsequences(VariantContext variant, SequenceI seq,
 826           SequenceFeature sf, int altAlelleIndex)
 827   {
 828     Object value = variant.getAttribute(CSQ_FIELD);
 829
 830     if (value == null || !(value instanceof ArrayList<?>))
 831     {
 832       return;
 833     }
 834
 835     List<String> consequences = (List<String>) value;
 836
 837     /*
 838      * if CSQ data includes 'Feature', and any value matches the sequence name,
 839      * then restrict consequence data to only the matching value (transcript)
 840      * i.e. just pick out consequences for the transcript the variant feature is on
 841      */
 842     String seqName = seq.getName()== null ? "" : seq.getName().toLowerCase();
 843     String matchFeature = null;
 844     if (csqFeatureFieldIndex > -1)
 845     {
 846       for (String consequence : consequences)
 847       {
 848         String[] csqFields = consequence.split(PIPE_REGEX);
 849         if (csqFields.length > csqFeatureFieldIndex)
 850         {
 851           String featureIdentifier = csqFields[csqFeatureFieldIndex];
 852           if (featureIdentifier.length() > 4
 853                   && seqName.indexOf(featureIdentifier.toLowerCase()) > -1)
 854           {
 855             matchFeature = featureIdentifier;
 856           }
 857         }
 858       }
 859     }
 860
 861     StringBuilder sb = new StringBuilder(128);
 862     boolean found = false;
 863
 864     // todo check against vepFieldsOfInterest as well somewhere
 865
 866     for (String consequence : consequences)
 867     {
 868       String[] csqFields = consequence.split(PIPE_REGEX);
 869
 870       if (includeConsequence(csqFields, matchFeature, variant,
 871               altAlelleIndex))
 872       {
 873         if (found)
 874         {
 875           sb.append(COMMA);
 876         }
 877         found = true;
 878         sb.append(consequence);
 879       }
 880     }
 881
 882     if (found)
 883     {
 884       sf.setValue(CSQ_FIELD, sb.toString());
 885     }
 886   }
 887
 888   /**
 889    * Answers true if we want to associate this block of consequence data with
 890    * the specified alternate allele of the VCF variant.
 891    * <p>
 892    * If consequence data includes the ALLELE_NUM field, then this has to match
 893    * altAlleleIndex. Otherwise the Allele field of the consequence data has to
 894    * match the allele value.
 895    * <p>
 896    * Optionally (if matchFeature is not null), restrict to only include
 897    * consequences whose Feature value matches. This allows us to attach
 898    * consequences to their respective transcripts.
 899    *
 900    * @param csqFields
 901    * @param matchFeature
 902    * @param variant
 903    * @param altAlelleIndex
 904    *          (0, 1..)
 905    * @return
 906    */
 907   protected boolean includeConsequence(String[] csqFields,
 908           String matchFeature, VariantContext variant, int altAlelleIndex)
 909   {
 910     /*
 911      * check consequence is for the current transcript
 912      */
 913     if (matchFeature != null)
 914     {
 915       if (csqFields.length <= csqFeatureFieldIndex)
 916       {
 917         return false;
 918       }
 919       String featureIdentifier = csqFields[csqFeatureFieldIndex];
 920       if (!featureIdentifier.equals(matchFeature))
 921       {
 922         return false; // consequence is for a different transcript
 923       }
 924     }
 925
 926     /*
 927      * if ALLELE_NUM is present, it must match altAlleleIndex
 928      * NB first alternate allele is 1 for ALLELE_NUM, 0 for altAlleleIndex
 929      */
 930     if (csqAlleleNumberFieldIndex > -1)
 931     {
 932       if (csqFields.length <= csqAlleleNumberFieldIndex)
 933       {
 934         return false;
 935       }
 936       String alleleNum = csqFields[csqAlleleNumberFieldIndex];
 937       return String.valueOf(altAlelleIndex + 1).equals(alleleNum);
 938     }
 939
 940     /*
 941      * else consequence allele must match variant allele
 942      */
 943     if (csqAlleleFieldIndex > -1 && csqFields.length > csqAlleleFieldIndex)
 944     {
 945       String csqAllele = csqFields[csqAlleleFieldIndex];
 946       String vcfAllele = variant.getAlternateAllele(altAlelleIndex)
 947               .getBaseString();
 948       return csqAllele.equals(vcfAllele);
 949     }
 950
 951     return false;
 952   }
 953
 954   /**
 955    * A convenience method to complement a dna base and return the string value
 956    * of its complement
 957    *
 958    * @param reference
 959    * @return
 960    */
 961   protected String complement(byte[] reference)
 962   {
 963     return String.valueOf(Dna.getComplement((char) reference[0]));
 964   }
 965
 966   /**
 967    * Determines the location of the query range (chromosome positions) in a
 968    * different reference assembly.
 969    * <p>
 970    * If the range is just a subregion of one for which we already have a mapping
 971    * (for example, an exon sub-region of a gene), then the mapping is just
 972    * computed arithmetically.
 973    * <p>
 974    * Otherwise, calls the Ensembl REST service that maps from one assembly
 975    * reference's coordinates to another's
 976    *
 977    * @param queryRange
 978    *          start-end chromosomal range in 'fromRef' coordinates
 979    * @param chromosome
 980    * @param species
 981    * @param fromRef
 982    *          assembly reference for the query coordinates
 983    * @param toRef
 984    *          assembly reference we wish to translate to
 985    * @return the start-end range in 'toRef' coordinates
 986    */
 987   protected int[] mapReferenceRange(int[] queryRange, String chromosome,
 988           String species, String fromRef, String toRef)
 989   {
 990     /*
 991      * first try shorcut of computing the mapping as a subregion of one
 992      * we already have (e.g. for an exon, if we have the gene mapping)
 993      */
 994     int[] mappedRange = findSubsumedRangeMapping(queryRange, chromosome,
 995             species, fromRef, toRef);
 996     if (mappedRange != null)
 997     {
 998       return mappedRange;
 999     }
1000
1001     /*
1002      * call (e.g.) http://rest.ensembl.org/map/human/GRCh38/17:45051610..45109016:1/GRCh37
1003      */
1004     EnsemblMap mapper = new EnsemblMap();
1005     int[] mapping = mapper.getAssemblyMapping(species, chromosome, fromRef,
1006             toRef, queryRange);
1007
1008     if (mapping == null)
1009     {
1010       // mapping service failure
1011       return null;
1012     }
1013
1014     /*
1015      * save mapping for possible future re-use
1016      */
1017     String key = makeRangesKey(chromosome, species, fromRef, toRef);
1018     if (!assemblyMappings.containsKey(key))
1019     {
1020       assemblyMappings.put(key, new HashMap<int[], int[]>());
1021     }
1022
1023     assemblyMappings.get(key).put(queryRange, mapping);
1024
1025     return mapping;
1026   }
1027
1028   /**
1029    * If we already have a 1:1 contiguous mapping which subsumes the given query
1030    * range, this method just calculates and returns the subset of that mapping,
1031    * else it returns null. In practical terms, if a gene has a contiguous
1032    * mapping between (for example) GRCh37 and GRCh38, then we assume that its
1033    * subsidiary exons occupy unchanged relative positions, and just compute
1034    * these as offsets, rather than do another lookup of the mapping.
1035    * <p>
1036    * If in future these assumptions prove invalid (e.g. for bacterial dna?!),
1037    * simply remove this method or let it always return null.
1038    * <p>
1039    * Warning: many rapid calls to the /map service map result in a 429 overload
1040    * error response
1041    *
1042    * @param queryRange
1043    * @param chromosome
1044    * @param species
1045    * @param fromRef
1046    * @param toRef
1047    * @return
1048    */
1049   protected int[] findSubsumedRangeMapping(int[] queryRange, String chromosome,
1050           String species, String fromRef, String toRef)
1051   {
1052     String key = makeRangesKey(chromosome, species, fromRef, toRef);
1053     if (assemblyMappings.containsKey(key))
1054     {
1055       Map<int[], int[]> mappedRanges = assemblyMappings.get(key);
1056       for (Entry<int[], int[]> mappedRange : mappedRanges.entrySet())
1057       {
1058         int[] fromRange = mappedRange.getKey();
1059         int[] toRange = mappedRange.getValue();
1060         if (fromRange[1] - fromRange[0] == toRange[1] - toRange[0])
1061         {
1062           /*
1063            * mapping is 1:1 in length, so we trust it to have no discontinuities
1064            */
1065           if (MappingUtils.rangeContains(fromRange, queryRange))
1066           {
1067             /*
1068              * fromRange subsumes our query range
1069              */
1070             int offset = queryRange[0] - fromRange[0];
1071             int mappedRangeFrom = toRange[0] + offset;
1072             int mappedRangeTo = mappedRangeFrom + (queryRange[1] - queryRange[0]);
1073             return new int[] { mappedRangeFrom, mappedRangeTo };
1074           }
1075         }
1076       }
1077     }
1078     return null;
1079   }
1080
1081   /**
1082    * Transfers the sequence feature to the target sequence, locating its start
1083    * and end range based on the mapping. Features which do not overlap the
1084    * target sequence are ignored.
1085    *
1086    * @param sf
1087    * @param targetSequence
1088    * @param mapping
1089    *          mapping from the feature's coordinates to the target sequence
1090    */
1091   protected void transferFeature(SequenceFeature sf,
1092           SequenceI targetSequence, MapList mapping)
1093   {
1094     int[] mappedRange = mapping.locateInTo(sf.getBegin(), sf.getEnd());
1095
1096     if (mappedRange != null)
1097     {
1098       String group = sf.getFeatureGroup();
1099       int newBegin = Math.min(mappedRange[0], mappedRange[1]);
1100       int newEnd = Math.max(mappedRange[0], mappedRange[1]);
1101       SequenceFeature copy = new SequenceFeature(sf, newBegin, newEnd,
1102               group, sf.getScore());
1103       targetSequence.addSequenceFeature(copy);
1104     }
1105   }
1106
1107   /**
1108    * Formats a ranges map lookup key
1109    *
1110    * @param chromosome
1111    * @param species
1112    * @param fromRef
1113    * @param toRef
1114    * @return
1115    */
1116   protected static String makeRangesKey(String chromosome, String species,
1117           String fromRef, String toRef)
1118   {
1119     return species + EXCL + chromosome + EXCL + fromRef + EXCL
1120             + toRef;
1121   }
1122 }