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