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