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