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