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
   9 import jalview.analysis.AlignmentUtils;
  10 import jalview.analysis.Dna;
  11 import jalview.api.AlignViewControllerGuiI;
  12 import jalview.datamodel.AlignmentI;
  13 import jalview.datamodel.DBRefEntry;
  14 import jalview.datamodel.GeneLoci;
  15 import jalview.datamodel.Mapping;
  16 import jalview.datamodel.SequenceFeature;
  17 import jalview.datamodel.SequenceI;
  18 import jalview.ext.ensembl.EnsemblMap;
  19 import jalview.ext.htsjdk.VCFReader;
  20 import jalview.io.gff.Gff3Helper;
  21 import jalview.io.gff.SequenceOntologyI;
  22 import jalview.util.MapList;
  23 import jalview.util.MappingUtils;
  24
  25 import java.io.IOException;
  26 import java.util.HashMap;
  27 import java.util.List;
  28 import java.util.Map;
  29 import java.util.Map.Entry;
  30
  31 /**
  32  * A class to read VCF data (using the htsjdk) and add variants as sequence
  33  * features on dna and any related protein product sequences
  34  *
  35  * @author gmcarstairs
  36  */
  37 public class VCFLoader
  38 {
  39   private static final String EXCL = "!";
  40
  41   /*
  42    * the alignment we are associated VCF data with
  43    */
  44   private AlignmentI al;
  45
  46   /*
  47    * mappings between VCF and sequence reference assembly regions, as
  48    * key = "species!chromosome!fromAssembly!toAssembly
  49    * value = Map{fromRange, toRange}
  50    */
  51   private Map<String, Map<int[], int[]>> assemblyMappings;
  52
  53   /**
  54    * Constructor given an alignment context
  55    *
  56    * @param alignment
  57    */
  58   public VCFLoader(AlignmentI alignment)
  59   {
  60     al = alignment;
  61
  62     // map of species!chromosome!fromAssembly!toAssembly to {fromRange, toRange}
  63     assemblyMappings = new HashMap<String, Map<int[], int[]>>();
  64   }
  65
  66   /**
  67    * Loads VCF on to an alignment - provided it can be related to one or more
  68    * sequence's chromosomal coordinates.
  69    * <p>
  70    * This method is not thread safe - concurrent threads should use separate
  71    * instances of this class.
  72    *
  73    * @param filePath
  74    * @param status
  75    */
  76   public void loadVCF(String filePath, AlignViewControllerGuiI status)
  77   {
  78     VCFReader reader = null;
  79
  80     try
  81     {
  82       // long start = System.currentTimeMillis();
  83       reader = new VCFReader(filePath);
  84
  85       VCFHeader header = reader.getFileHeader();
  86       VCFHeaderLine ref = header
  87               .getOtherHeaderLine(VCFHeader.REFERENCE_KEY);
  88       // check if reference is wrt assembly19 (GRCh37)
  89       // todo may need to allow user to specify reference assembly?
  90       boolean isRefGrch37 = (ref != null && ref.getValue().contains(
  91               "assembly19"));
  92
  93       int varCount = 0;
  94       int seqCount = 0;
  95
  96       /*
  97        * query for VCF overlapping each sequence in turn
  98        */
  99       for (SequenceI seq : al.getSequences())
 100       {
 101         int added = loadVCF(seq, reader, isRefGrch37);
 102         if (added > 0)
 103         {
 104           seqCount++;
 105           varCount += added;
 106           computePeptideVariants(seq);
 107         }
 108       }
 109       // long elapsed = System.currentTimeMillis() - start;
 110       String msg = String.format("Added %d VCF variants to %d sequence(s)",
 111               varCount, seqCount);
 112       if (status != null)
 113       {
 114         status.setStatus(msg);
 115       }
 116     } catch (Throwable e)
 117     {
 118       System.err.println("Error processing VCF: " + e.getMessage());
 119       e.printStackTrace();
 120     } finally
 121     {
 122       if (reader != null)
 123       {
 124         try
 125         {
 126           reader.close();
 127         } catch (IOException e)
 128         {
 129           // ignore
 130         }
 131       }
 132     }
 133   }
 134
 135   /**
 136    * (Re-)computes peptide variants from dna variants, for any protein sequence
 137    * to which the dna sequence has a mapping. Note that although duplicate
 138    * features may get computed, they will not be added, since duplicate sequence
 139    * features are ignored in Sequence.addSequenceFeature.
 140    *
 141    * @param dnaSeq
 142    */
 143   protected void computePeptideVariants(SequenceI dnaSeq)
 144   {
 145     DBRefEntry[] dbrefs = dnaSeq.getDBRefs();
 146     if (dbrefs == null)
 147     {
 148       return;
 149     }
 150     for (DBRefEntry dbref : dbrefs)
 151     {
 152       Mapping mapping = dbref.getMap();
 153       if (mapping == null || mapping.getTo() == null
 154               || mapping.getMap().getFromRatio() != 3)
 155       {
 156         continue;
 157       }
 158       AlignmentUtils.computeProteinFeatures(dnaSeq, mapping.getTo(),
 159               mapping.getMap());
 160     }
 161   }
 162
 163   /**
 164    * Tries to add overlapping variants read from a VCF file to the given
 165    * sequence, and returns the number of overlapping variants found. Note that
 166    * this requires the sequence to hold information as to its chromosomal
 167    * positions and reference, in order to be able to map the VCF variants to the
 168    * sequence.
 169    *
 170    * @param seq
 171    * @param reader
 172    * @param isVcfRefGrch37
 173    * @return
 174    */
 175   protected int loadVCF(SequenceI seq, VCFReader reader,
 176           boolean isVcfRefGrch37)
 177   {
 178     int count = 0;
 179     GeneLoci seqCoords = seq.getGeneLoci();
 180     if (seqCoords == null)
 181     {
 182       return 0;
 183     }
 184
 185     List<int[]> seqChromosomalContigs = seqCoords.mapping.getToRanges();
 186     for (int[] range : seqChromosomalContigs)
 187     {
 188       count += addVcfVariants(seq, reader, range, isVcfRefGrch37);
 189     }
 190
 191     return count;
 192   }
 193
 194   /**
 195    * Queries the VCF reader for any variants that overlap the given chromosome
 196    * region of the sequence, and adds as variant features. Returns the number of
 197    * overlapping variants found.
 198    *
 199    * @param seq
 200    * @param reader
 201    * @param range
 202    *          start-end range of a sequence region in its chromosomal
 203    *          coordinates
 204    * @param isVcfRefGrch37
 205    *          true if the VCF is with reference to GRCh37
 206    * @return
 207    */
 208   protected int addVcfVariants(SequenceI seq, VCFReader reader,
 209           int[] range, boolean isVcfRefGrch37)
 210   {
 211     GeneLoci seqCoords = seq.getGeneLoci();
 212
 213     String chromosome = seqCoords.chromosome;
 214     String seqRef = seqCoords.assembly;
 215     String species = seqCoords.species;
 216
 217     // TODO handle species properly
 218     if ("".equals(species))
 219     {
 220       species = "human";
 221     }
 222
 223     /*
 224      * map chromosomal coordinates from GRCh38 (sequence) to
 225      * GRCh37 (VCF) if necessary
 226      */
 227     // TODO generalise for other assemblies and species
 228     int offset = 0;
 229     String fromRef = "GRCh38";
 230     if (fromRef.equalsIgnoreCase(seqRef) && isVcfRefGrch37)
 231     {
 232       String toRef = "GRCh37";
 233       int[] newRange = mapReferenceRange(range, chromosome, species,
 234               fromRef, toRef);
 235       if (newRange == null)
 236       {
 237         System.err.println(String.format(
 238                 "Failed to map %s:%s:%s:%d:%d to %s", species, chromosome,
 239                 fromRef, range[0], range[1], toRef));
 240         return 0;
 241       }
 242       offset = newRange[0] - range[0];
 243       range = newRange;
 244     }
 245
 246     boolean forwardStrand = range[0] <= range[1];
 247
 248     /*
 249      * query the VCF for overlaps
 250      * (convert a reverse strand range to forwards)
 251      */
 252     int count = 0;
 253     MapList mapping = seqCoords.mapping;
 254
 255     int fromLocus = Math.min(range[0], range[1]);
 256     int toLocus = Math.max(range[0], range[1]);
 257     CloseableIterator<VariantContext> variants = reader.query(chromosome,
 258             fromLocus, toLocus);
 259     while (variants.hasNext())
 260     {
 261       /*
 262        * get variant location in sequence chromosomal coordinates
 263        */
 264       VariantContext variant = variants.next();
 265
 266       /*
 267        * we can only process SNP variants (which can be reported
 268        * as part of a MIXED variant record
 269        */
 270       if (!variant.isSNP() && !variant.isMixed())
 271       {
 272         continue;
 273       }
 274
 275       count++;
 276       int start = variant.getStart() - offset;
 277       int end = variant.getEnd() - offset;
 278
 279       /*
 280        * convert chromosomal location to sequence coordinates
 281        * - null if a partially overlapping feature
 282        */
 283       int[] seqLocation = mapping.locateInFrom(start, end);
 284       if (seqLocation != null)
 285       {
 286         addVariantFeatures(seq, variant, seqLocation[0], seqLocation[1],
 287                 forwardStrand);
 288       }
 289     }
 290
 291     variants.close();
 292
 293     return count;
 294   }
 295
 296   /**
 297    * Inspects the VCF variant record, and adds variant features to the sequence.
 298    * Only SNP variants are added, not INDELs.
 299    * <p>
 300    * If the sequence maps to the reverse strand of the chromosome, reference and
 301    * variant bases are recorded as their complements (C/G, A/T).
 302    *
 303    * @param seq
 304    * @param variant
 305    * @param featureStart
 306    * @param featureEnd
 307    * @param forwardStrand
 308    */
 309   protected void addVariantFeatures(SequenceI seq, VariantContext variant,
 310           int featureStart, int featureEnd, boolean forwardStrand)
 311   {
 312     byte[] reference = variant.getReference().getBases();
 313     if (reference.length != 1)
 314     {
 315       /*
 316        * sorry, we don't handle INDEL variants
 317        */
 318       return;
 319     }
 320
 321     /*
 322      * for now we extract allele frequency as feature score; note
 323      * this attribute is String for a simple SNP, but List<String> if
 324      * multiple alleles at the locus; we extract for the simple case only
 325      */
 326     Object af = variant.getAttribute("AF");
 327     float score = 0f;
 328     if (af instanceof String)
 329     {
 330       try
 331       {
 332         score = Float.parseFloat((String) af);
 333       } catch (NumberFormatException e)
 334       {
 335         // leave as 0
 336       }
 337     }
 338
 339     StringBuilder sb = new StringBuilder();
 340     sb.append(forwardStrand ? (char) reference[0] : complement(reference));
 341
 342     /*
 343      * inspect alleles and record SNP variants (as the variant
 344      * record could be MIXED and include INDEL and SNP alleles)
 345      * warning: getAlleles gives no guarantee as to the order
 346      * in which they are returned
 347      */
 348     for (Allele allele : variant.getAlleles())
 349     {
 350       if (!allele.isReference())
 351       {
 352         byte[] alleleBase = allele.getBases();
 353         if (alleleBase.length == 1)
 354         {
 355           sb.append(",").append(
 356                   forwardStrand ? (char) alleleBase[0]
 357                           : complement(alleleBase));
 358         }
 359       }
 360     }
 361     String alleles = sb.toString(); // e.g. G,A,C
 362
 363     String type = SequenceOntologyI.SEQUENCE_VARIANT;
 364
 365     SequenceFeature sf = new SequenceFeature(type, alleles, featureStart,
 366             featureEnd, score, "VCF");
 367
 368     sf.setValue(Gff3Helper.ALLELES, alleles);
 369
 370     Map<String, Object> atts = variant.getAttributes();
 371     for (Entry<String, Object> att : atts.entrySet())
 372     {
 373       sf.setValue(att.getKey(), att.getValue());
 374     }
 375     seq.addSequenceFeature(sf);
 376   }
 377
 378   /**
 379    * A convenience method to complement a dna base and return the string value
 380    * of its complement
 381    *
 382    * @param reference
 383    * @return
 384    */
 385   protected String complement(byte[] reference)
 386   {
 387     return String.valueOf(Dna.getComplement((char) reference[0]));
 388   }
 389
 390   /**
 391    * Determines the location of the query range (chromosome positions) in a
 392    * different reference assembly.
 393    * <p>
 394    * If the range is just a subregion of one for which we already have a mapping
 395    * (for example, an exon sub-region of a gene), then the mapping is just
 396    * computed arithmetically.
 397    * <p>
 398    * Otherwise, calls the Ensembl REST service that maps from one assembly
 399    * reference's coordinates to another's
 400    *
 401    * @param queryRange
 402    *          start-end chromosomal range in 'fromRef' coordinates
 403    * @param chromosome
 404    * @param species
 405    * @param fromRef
 406    *          assembly reference for the query coordinates
 407    * @param toRef
 408    *          assembly reference we wish to translate to
 409    * @return the start-end range in 'toRef' coordinates
 410    */
 411   protected int[] mapReferenceRange(int[] queryRange, String chromosome,
 412           String species, String fromRef, String toRef)
 413   {
 414     /*
 415      * first try shorcut of computing the mapping as a subregion of one
 416      * we already have (e.g. for an exon, if we have the gene mapping)
 417      */
 418     int[] mappedRange = findSubsumedRangeMapping(queryRange, chromosome,
 419             species, fromRef, toRef);
 420     if (mappedRange != null)
 421     {
 422       return mappedRange;
 423     }
 424
 425     /*
 426      * call (e.g.) http://rest.ensembl.org/map/human/GRCh38/17:45051610..45109016:1/GRCh37
 427      */
 428     EnsemblMap mapper = new EnsemblMap();
 429     int[] mapping = mapper.getMapping(species, chromosome, fromRef, toRef,
 430             queryRange);
 431
 432     if (mapping == null)
 433     {
 434       // mapping service failure
 435       return null;
 436     }
 437
 438     /*
 439      * save mapping for possible future re-use
 440      */
 441     String key = makeRangesKey(chromosome, species, fromRef, toRef);
 442     if (!assemblyMappings.containsKey(key))
 443     {
 444       assemblyMappings.put(key, new HashMap<int[], int[]>());
 445     }
 446
 447     assemblyMappings.get(key).put(queryRange, mapping);
 448
 449     return mapping;
 450   }
 451
 452   /**
 453    * If we already have a 1:1 contiguous mapping which subsumes the given query
 454    * range, this method just calculates and returns the subset of that mapping,
 455    * else it returns null. In practical terms, if a gene has a contiguous
 456    * mapping between (for example) GRCh37 and GRCh38, then we assume that its
 457    * subsidiary exons occupy unchanged relative positions, and just compute
 458    * these as offsets, rather than do another lookup of the mapping.
 459    * <p>
 460    * If in future these assumptions prove invalid (e.g. for bacterial dna?!),
 461    * simply remove this method or let it always return null.
 462    * <p>
 463    * Warning: many rapid calls to the /map service map result in a 429 overload
 464    * error response
 465    *
 466    * @param queryRange
 467    * @param chromosome
 468    * @param species
 469    * @param fromRef
 470    * @param toRef
 471    * @return
 472    */
 473   protected int[] findSubsumedRangeMapping(int[] queryRange, String chromosome,
 474           String species, String fromRef, String toRef)
 475   {
 476     String key = makeRangesKey(chromosome, species, fromRef, toRef);
 477     if (assemblyMappings.containsKey(key))
 478     {
 479       Map<int[], int[]> mappedRanges = assemblyMappings.get(key);
 480       for (Entry<int[], int[]> mappedRange : mappedRanges.entrySet())
 481       {
 482         int[] fromRange = mappedRange.getKey();
 483         int[] toRange = mappedRange.getValue();
 484         if (fromRange[1] - fromRange[0] == toRange[1] - toRange[0])
 485         {
 486           /*
 487            * mapping is 1:1 in length, so we trust it to have no discontinuities
 488            */
 489           if (MappingUtils.rangeContains(fromRange, queryRange))
 490           {
 491             /*
 492              * fromRange subsumes our query range
 493              */
 494             int offset = queryRange[0] - fromRange[0];
 495             int mappedRangeFrom = toRange[0] + offset;
 496             int mappedRangeTo = mappedRangeFrom + (queryRange[1] - queryRange[0]);
 497             return new int[] { mappedRangeFrom, mappedRangeTo };
 498           }
 499         }
 500       }
 501     }
 502     return null;
 503   }
 504
 505   /**
 506    * Formats a ranges map lookup key
 507    *
 508    * @param chromosome
 509    * @param species
 510    * @param fromRef
 511    * @param toRef
 512    * @return
 513    */
 514   protected static String makeRangesKey(String chromosome, String species,
 515           String fromRef, String toRef)
 516   {
 517     return species + EXCL + chromosome + EXCL + fromRef + EXCL
 518             + toRef;
 519   }
 520 }