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