src/jalview/ext/htsjdk/VCFReader.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
   3  * Copyright (C) $$Year-Rel$$ The Jalview Authors
   4  *
   5  * This file is part of Jalview.
   6  *
   7  * Jalview is free software: you can redistribute it and/or
   8  * modify it under the terms of the GNU General Public License
   9  * as published by the Free Software Foundation, either version 3
  10  * of the License, or (at your option) any later version.
  11  *
  12  * Jalview is distributed in the hope that it will be useful, but
  13  * WITHOUT ANY WARRANTY; without even the implied warranty
  14  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  15  * PURPOSE.  See the GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
  19  * The Jalview Authors are detailed in the 'AUTHORS' file.
  20  */
  21 package jalview.ext.htsjdk;
  22
  23 import jalview.bin.Cache;
  24
  25 import java.io.Closeable;
  26 import java.io.File;
  27 import java.io.IOException;
  28
  29 import htsjdk.samtools.util.CloseableIterator;
  30 import htsjdk.variant.variantcontext.VariantContext;
  31 import htsjdk.variant.vcf.VCFFileReader;
  32 import htsjdk.variant.vcf.VCFHeader;
  33
  34 /**
  35  * A thin wrapper for htsjdk classes to read either plain, or compressed, or
  36  * compressed and indexed VCF files
  37  */
  38 public class VCFReader implements Closeable, Iterable<VariantContext>
  39 {
  40   private static final String GZ = "gz";
  41
  42   private static final String TBI_EXTENSION = ".tbi";
  43
  44   private static final String CSI_EXTENSION = ".csi";
  45
  46   private boolean indexed;
  47
  48   private VCFFileReader reader;
  49
  50   /**
  51    * Constructor given a raw or compressed VCF file or a (csi or tabix) index file
  52    * <p>
  53    * If the file path ends in ".tbi" or ".csi", <em>or</em> appending one of these
  54    * extensions gives a valid file path, open as indexed, else as unindexed.
  55    *
  56    * @param f
  57    * @throws IOException
  58    */
  59   public VCFReader(String filePath) throws IOException
  60   {
  61     indexed = false;
  62     if (filePath.endsWith(TBI_EXTENSION)
  63             || filePath.endsWith(CSI_EXTENSION))
  64     {
  65       indexed = true;
  66       filePath = filePath.substring(0, filePath.length() - 4);
  67     }
  68     else if (new File(filePath + TBI_EXTENSION).exists())
  69     {
  70       indexed = true;
  71     }
  72     else if (new File(filePath + CSI_EXTENSION).exists())
  73     {
  74       indexed = true;
  75     }
  76
  77     /*
  78      * we pass the name of the unindexed file to htsjdk,
  79      * with a flag to assert whether it is indexed
  80      */
  81     File file = new File(filePath);
  82     if (file.exists())
  83     {
  84       reader = new VCFFileReader(file, indexed);
  85     }
  86     else
  87     {
  88       Cache.error("File not found: " + filePath);
  89     }
  90   }
  91
  92   @Override
  93   public void close() throws IOException
  94   {
  95     if (reader != null)
  96     {
  97       reader.close();
  98     }
  99   }
 100
 101   /**
 102    * Returns an iterator over VCF variants in the file. The client should call
 103    * close() on the iterator when finished with it.
 104    */
 105   @Override
 106   public CloseableIterator<VariantContext> iterator()
 107   {
 108     return reader == null ? null : reader.iterator();
 109   }
 110
 111   /**
 112    * Queries for records overlapping the region specified. Note that this method
 113    * is performant if the VCF file is indexed, and may be very slow if it is
 114    * not.
 115    * <p>
 116    * Client code should call close() on the iterator when finished with it.
 117    *
 118    * @param chrom
 119    *          the chromosome to query
 120    * @param start
 121    *          query interval start
 122    * @param end
 123    *          query interval end
 124    * @return
 125    */
 126   public CloseableIterator<VariantContext> query(final String chrom,
 127           final int start, final int end)
 128   {
 129     if (reader == null)
 130     {
 131       return null;
 132     }
 133     if (indexed)
 134     {
 135       return reader.query(chrom, start, end);
 136     }
 137     else
 138     {
 139       return queryUnindexed(chrom, start, end);
 140     }
 141   }
 142
 143   /**
 144    * Returns an iterator over variant records read from a flat file which
 145    * overlap the specified chromosomal positions. Call close() on the iterator
 146    * when finished with it!
 147    *
 148    * @param chrom
 149    * @param start
 150    * @param end
 151    * @return
 152    */
 153   protected CloseableIterator<VariantContext> queryUnindexed(
 154           final String chrom, final int start, final int end)
 155   {
 156     final CloseableIterator<VariantContext> it = reader.iterator();
 157
 158     return new CloseableIterator<VariantContext>()
 159     {
 160       boolean atEnd = false;
 161
 162       // prime look-ahead buffer with next matching record
 163       private VariantContext next = findNext();
 164
 165       private VariantContext findNext()
 166       {
 167         if (atEnd)
 168         {
 169           return null;
 170         }
 171         VariantContext variant = null;
 172         while (it.hasNext())
 173         {
 174           variant = it.next();
 175           int vstart = variant.getStart();
 176
 177           if (vstart > end)
 178           {
 179             atEnd = true;
 180             close();
 181             return null;
 182           }
 183
 184           int vend = variant.getEnd();
 185           // todo what is the undeprecated way to get
 186           // the chromosome for the variant?
 187           if (chrom.equals(variant.getContig()) && (vstart <= end)
 188                   && (vend >= start))
 189           {
 190             return variant;
 191           }
 192         }
 193         return null;
 194       }
 195
 196       @Override
 197       public boolean hasNext()
 198       {
 199         boolean hasNext = !atEnd && (next != null);
 200         if (!hasNext)
 201         {
 202           close();
 203         }
 204         return hasNext;
 205       }
 206
 207       @Override
 208       public VariantContext next()
 209       {
 210         /*
 211          * return the next match, and then re-prime
 212          * it with the following one (if any)
 213          */
 214         VariantContext temp = next;
 215         next = findNext();
 216         return temp;
 217       }
 218
 219       @Override
 220       public void remove()
 221       {
 222         // not implemented
 223       }
 224
 225       @Override
 226       public void close()
 227       {
 228         it.close();
 229       }
 230     };
 231   }
 232
 233   /**
 234    * Returns an object that models the VCF file headers
 235    *
 236    * @return
 237    */
 238   public VCFHeader getFileHeader()
 239   {
 240     return reader == null ? null : reader.getFileHeader();
 241   }
 242
 243   /**
 244    * Answers true if we are processing a tab-indexed VCF file, false if it is a
 245    * plain text (uncompressed) file.
 246    *
 247    * @return
 248    */
 249   public boolean isIndex()
 250   {
 251     return indexed;
 252   }
 253 }