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