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 ? reference : Dna
 341             .getComplement((char) reference[0]));
 342
 343     /*
 344      * inspect alleles and record SNP variants (as the variant
 345      * record could be MIXED and include INDEL and SNP alleles)
 346      * warning: getAlleles gives no guarantee as to the order
 347      * in which they are returned
 348      */
 349     for (Allele allele : variant.getAlleles())
 350     {
 351       if (!allele.isReference())
 352       {
 353         byte[] alleleBase = allele.getBases();
 354         if (alleleBase.length == 1)
 355         {
 356           sb.append(",").append(
 357                   forwardStrand ? alleleBase : Dna
 358                           .getComplement((char) alleleBase[0]));
 359         }
 360       }
 361     }
 362     String alleles = sb.toString(); // e.g. G,A,C
 363
 364     String type = SequenceOntologyI.SEQUENCE_VARIANT;
 365
 366     SequenceFeature sf = new SequenceFeature(type, alleles, featureStart,
 367             featureEnd, score, "VCF");
 368
 369     sf.setValue(Gff3Helper.ALLELES, alleles);
 370
 371     Map<String, Object> atts = variant.getAttributes();
 372     for (Entry<String, Object> att : atts.entrySet())
 373     {
 374       sf.setValue(att.getKey(), att.getValue());
 375     }
 376     seq.addSequenceFeature(sf);
 377   }
 378
 379   /**
 380    * Determines the location of the query range (chromosome positions) in a
 381    * different reference assembly.
 382    * <p>
 383    * If the range is just a subregion of one for which we already have a mapping
 384    * (for example, an exon sub-region of a gene), then the mapping is just
 385    * computed arithmetically.
 386    * <p>
 387    * Otherwise, calls the Ensembl REST service that maps from one assembly
 388    * reference's coordinates to another's
 389    *
 390    * @param queryRange
 391    *          start-end chromosomal range in 'fromRef' coordinates
 392    * @param chromosome
 393    * @param species
 394    * @param fromRef
 395    *          assembly reference for the query coordinates
 396    * @param toRef
 397    *          assembly reference we wish to translate to
 398    * @return the start-end range in 'toRef' coordinates
 399    */
 400   protected int[] mapReferenceRange(int[] queryRange, String chromosome,
 401           String species, String fromRef, String toRef)
 402   {
 403     /*
 404      * first try shorcut of computing the mapping as a subregion of one
 405      * we already have (e.g. for an exon, if we have the gene mapping)
 406      */
 407     int[] mappedRange = findSubsumedRangeMapping(queryRange, chromosome,
 408             species, fromRef, toRef);
 409     if (mappedRange != null)
 410     {
 411       return mappedRange;
 412     }
 413
 414     /*
 415      * call (e.g.) http://rest.ensembl.org/map/human/GRCh38/17:45051610..45109016:1/GRCh37
 416      */
 417     EnsemblMap mapper = new EnsemblMap();
 418     int[] mapping = mapper.getMapping(species, chromosome, fromRef, toRef,
 419             queryRange);
 420
 421     if (mapping == null)
 422     {
 423       // mapping service failure
 424       return null;
 425     }
 426
 427     /*
 428      * save mapping for possible future re-use
 429      */
 430     String key = makeRangesKey(chromosome, species, fromRef, toRef);
 431     if (!assemblyMappings.containsKey(key))
 432     {
 433       assemblyMappings.put(key, new HashMap<int[], int[]>());
 434     }
 435
 436     assemblyMappings.get(key).put(queryRange, mapping);
 437
 438     return mapping;
 439   }
 440
 441   /**
 442    * If we already have a 1:1 contiguous mapping which subsumes the given query
 443    * range, this method just calculates and returns the subset of that mapping,
 444    * else it returns null. In practical terms, if a gene has a contiguous
 445    * mapping between (for example) GRCh37 and GRCh38, then we assume that its
 446    * subsidiary exons occupy unchanged relative positions, and just compute
 447    * these as offsets, rather than do another lookup of the mapping.
 448    * <p>
 449    * If in future these assumptions prove invalid (e.g. for bacterial dna?!),
 450    * simply remove this method or let it always return null.
 451    * <p>
 452    * Warning: many rapid calls to the /map service map result in a 429 overload
 453    * error response
 454    *
 455    * @param queryRange
 456    * @param chromosome
 457    * @param species
 458    * @param fromRef
 459    * @param toRef
 460    * @return
 461    */
 462   protected int[] findSubsumedRangeMapping(int[] queryRange, String chromosome,
 463           String species, String fromRef, String toRef)
 464   {
 465     String key = makeRangesKey(chromosome, species, fromRef, toRef);
 466     if (assemblyMappings.containsKey(key))
 467     {
 468       Map<int[], int[]> mappedRanges = assemblyMappings.get(key);
 469       for (Entry<int[], int[]> mappedRange : mappedRanges.entrySet())
 470       {
 471         int[] fromRange = mappedRange.getKey();
 472         int[] toRange = mappedRange.getValue();
 473         if (fromRange[1] - fromRange[0] == toRange[1] - toRange[0])
 474         {
 475           /*
 476            * mapping is 1:1 in length, so we trust it to have no discontinuities
 477            */
 478           if (MappingUtils.rangeContains(fromRange, queryRange))
 479           {
 480             /*
 481              * fromRange subsumes our query range
 482              */
 483             int offset = queryRange[0] - fromRange[0];
 484             int mappedRangeFrom = toRange[0] + offset;
 485             int mappedRangeTo = mappedRangeFrom + (queryRange[1] - queryRange[0]);
 486             return new int[] { mappedRangeFrom, mappedRangeTo };
 487           }
 488         }
 489       }
 490     }
 491     return null;
 492   }
 493
 494   /**
 495    * Formats a ranges map lookup key
 496    *
 497    * @param chromosome
 498    * @param species
 499    * @param fromRef
 500    * @param toRef
 501    * @return
 502    */
 503   protected static String makeRangesKey(String chromosome, String species,
 504           String fromRef, String toRef)
 505   {
 506     return species + EXCL + chromosome + EXCL + fromRef + EXCL
 507             + toRef;
 508   }
 509 }