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.htsjdk.VCFReader;
  18 import jalview.io.gff.SequenceOntologyI;
  19 import jalview.util.MapList;
  20
  21 import java.util.List;
  22 import java.util.Map;
  23 import java.util.Map.Entry;
  24
  25 public class VCFLoader
  26 {
  27   AlignmentI al;
  28
  29   /**
  30    * Constructor given an alignment context
  31    *
  32    * @param alignment
  33    */
  34   public VCFLoader(AlignmentI alignment)
  35   {
  36     al = alignment;
  37   }
  38
  39   /**
  40    * Loads VCF on to an alignment - provided it can be related to one or more
  41    * sequence's chromosomal coordinates
  42    *
  43    * @param filePath
  44    */
  45   public void loadVCF(String filePath)
  46   {
  47     VCFReader reader = null;
  48
  49     try
  50     {
  51       long start = System.currentTimeMillis();
  52       reader = new VCFReader(filePath);
  53
  54       VCFHeader header = reader.getFileHeader();
  55       VCFHeaderLine ref = header
  56               .getOtherHeaderLine(VCFHeader.REFERENCE_KEY);
  57       // check if reference is wrt assembly19 (GRCh37)
  58       // todo may need to allow user to specify reference assembly?
  59       boolean isRefGrch37 = (ref != null && ref.getValue().contains(
  60               "assembly19"));
  61
  62       int varCount = 0;
  63       int seqCount = 0;
  64
  65       /*
  66        * query for VCF overlapping each sequence in turn
  67        */
  68       for (SequenceI seq : al.getSequences())
  69       {
  70         int added = loadVCF(seq, reader, isRefGrch37);
  71         if (added > 0)
  72         {
  73           seqCount++;
  74           varCount += added;
  75           computePeptideVariants(seq);
  76         }
  77       }
  78       long elapsed = System.currentTimeMillis() - start;
  79       System.out.println(String.format(
  80               "Added %d VCF variants to %d sequence(s) (%dms)", varCount,
  81               seqCount, elapsed));
  82
  83       reader.close();
  84     } catch (Throwable e)
  85     {
  86       System.err.println("Error processing VCF: " + e.getMessage());
  87       e.printStackTrace();
  88     }
  89   }
  90
  91   /**
  92    * (Re-)computes peptide variants from dna variants, for any protein sequence
  93    * to which the dna sequence has a mapping. Note that although duplicate
  94    * features may get computed, they will not be added, since duplicate sequence
  95    * features are ignored in Sequence.addSequenceFeature.
  96    *
  97    * @param dnaSeq
  98    */
  99   protected void computePeptideVariants(SequenceI dnaSeq)
 100   {
 101     DBRefEntry[] dbrefs = dnaSeq.getDBRefs();
 102     if (dbrefs == null)
 103     {
 104       return;
 105     }
 106     for (DBRefEntry dbref : dbrefs)
 107     {
 108       Mapping mapping = dbref.getMap();
 109       if (mapping == null || mapping.getTo() == null
 110               || mapping.getMap().getFromRatio() != 3)
 111       {
 112         continue;
 113       }
 114       AlignmentUtils.computeProteinFeatures(dnaSeq, mapping.getTo(),
 115               mapping.getMap());
 116     }
 117   }
 118
 119   /**
 120    * Tries to add overlapping variants read from a VCF file to the given
 121    * sequence, and returns the number of overlapping variants found. Note that
 122    * this requires the sequence to hold information as to its chromosomal
 123    * positions and reference, in order to be able to map the VCF variants to the
 124    * sequence.
 125    *
 126    * @param seq
 127    * @param reader
 128    * @param isVcfRefGrch37
 129    * @return
 130    */
 131   protected int loadVCF(SequenceI seq, VCFReader reader, boolean isVcfRefGrch37)
 132   {
 133     int count = 0;
 134     GeneLoci seqCoords = ((Sequence) seq).getGeneLoci();
 135     if (seqCoords == null)
 136     {
 137       return 0;
 138     }
 139
 140     MapList mapping = seqCoords.getMapping();
 141     List<int[]> seqChromosomalContigs = mapping.getToRanges();
 142     for (int[] range : seqChromosomalContigs)
 143     {
 144       count += addVcfVariants(seq, reader, range, isVcfRefGrch37);
 145     }
 146
 147     return count;
 148   }
 149
 150   /**
 151    * Queries the VCF reader for any variants that overlap the given chromosome
 152    * region of the sequence, and adds as variant features. Returns the number of
 153    * overlapping variants found.
 154    *
 155    * @param seq
 156    * @param reader
 157    * @param range
 158    *          start-end range of a sequence region in its chromosomal
 159    *          coordinates
 160    * @param isVcfRefGrch37
 161    *          true if the VCF is with reference to GRCh37
 162    * @return
 163    */
 164   protected int addVcfVariants(SequenceI seq, VCFReader reader,
 165           int[] range, boolean isVcfRefGrch37)
 166   {
 167     GeneLoci seqCoords = ((Sequence) seq).getGeneLoci();
 168
 169     String chromosome = seqCoords.getChromosome();
 170     String seqRef = seqCoords.getReference();
 171     String species = seqCoords.getSpecies();
 172
 173     /*
 174      * map chromosomal coordinates from GRCh38 (sequence) to
 175      * GRCh37 (VCF) if necessary
 176      */
 177     int offset = 0;
 178     if ("GRCh38".equalsIgnoreCase(seqRef) && isVcfRefGrch37)
 179     {
 180       int[] newRange = mapReferenceRange(range, chromosome, species,
 181               "GRCh38",
 182               "GRCh37");
 183       offset = newRange[0] - range[0];
 184       range = newRange;
 185     }
 186
 187     /*
 188      * query the VCF for overlaps
 189      */
 190     int count = 0;
 191     MapList mapping = seqCoords.getMapping();
 192
 193     CloseableIterator<VariantContext> variants = reader.query(chromosome,
 194             range[0], range[1]);
 195     while (variants.hasNext())
 196     {
 197       /*
 198        * get variant location in sequence chromosomal coordinates
 199        */
 200       VariantContext variant = variants.next();
 201       count++;
 202       int start = variant.getStart() - offset;
 203       int end = variant.getEnd() - offset;
 204
 205       /*
 206        * get location in sequence coordinates
 207        */
 208       int[] seqLocation = mapping.locateInFrom(start, end);
 209       int featureStart = seqLocation[0];
 210       int featureEnd = seqLocation[1];
 211       addVariantFeatures(seq, variant, featureStart, featureEnd);
 212     }
 213
 214     variants.close();
 215
 216     return count;
 217   }
 218
 219   /**
 220    * Inspects the VCF variant record, and adds variant features to the sequence
 221    *
 222    * @param seq
 223    * @param variant
 224    * @param featureStart
 225    * @param featureEnd
 226    */
 227   protected void addVariantFeatures(SequenceI seq, VariantContext variant,
 228           int featureStart, int featureEnd)
 229   {
 230     StringBuilder sb = new StringBuilder();
 231     sb.append(variant.getReference().getBaseString());
 232
 233     int alleleCount = 0;
 234     for (Allele allele : variant.getAlleles())
 235     {
 236       if (!allele.isReference())
 237       {
 238         sb.append(",").append(allele.getBaseString());
 239         alleleCount++;
 240       }
 241     }
 242     String alleles = sb.toString(); // e.g. G,A,C
 243
 244     String type = SequenceOntologyI.SEQUENCE_VARIANT;
 245     float score = 0f;
 246     if (alleleCount == 1)
 247     {
 248       try
 249       {
 250         score = (float) variant.getAttributeAsDouble("AF", 0d);
 251       } catch (NumberFormatException e)
 252       {
 253         // leave score as 0
 254       }
 255     }
 256     SequenceFeature sf = new SequenceFeature(type, alleles, featureStart,
 257             featureEnd, score, "VCF");
 258
 259     /*
 260      * only add 'alleles' property if a SNP, as we can
 261      * only handle SNPs when computing peptide variants
 262      */
 263     if (variant.isSNP())
 264     {
 265       sf.setValue("alleles", alleles);
 266     }
 267
 268     Map<String, Object> atts = variant.getAttributes();
 269     for (Entry<String, Object> att : atts.entrySet())
 270     {
 271       sf.setValue(att.getKey(), att.getValue());
 272     }
 273     seq.addSequenceFeature(sf);
 274   }
 275
 276   /**
 277    * Call the Ensembl REST service that maps from one assembly reference's
 278    * coordinates to another's
 279    *
 280    * @param range
 281    * @param chromosome
 282    * @param species
 283    * @param fromRef
 284    * @param toRef
 285    * @return
 286    */
 287   protected int[] mapReferenceRange(int[] range, String chromosome,
 288           String species, String fromRef, String toRef)
 289   {
 290     // TODO call
 291     // http://rest.ensembl.org/map/species/fromRef/chromosome:range[0]..range[1]:1/toRef?content-type=application/json
 292     // and parse the JSON response
 293
 294     // 1922632 is the difference between 37 and 38 for chromosome 17
 295     return new int[] { range[0] - 1922632, range[1] - 1922632 };
 296   }
 297 }