Added experimental support for reading bcf files.

vcf/parser.py

@@ -6,6 +6,7 @@
 import os
 import re
 import sys
+import subprocess
 
 try:
     from collections import OrderedDict
@@ -261,8 +262,15 @@ def __init__(self, fsock=None, filename=None, compressed=False, prepend_chr=Fals
                 filename = fsock.name
                 compressed = compressed or filename.endswith('.gz')
         elif filename:
-            compressed = compressed or filename.endswith('.gz')
-            self._reader = open(filename, 'rb' if compressed else 'rt')
+            if filename.endswith(".bcf"):
+                if which("bcftools") is None:
+                    raise Exception('bcftools not available. Please install bcftools')
+                if not os.path.isfile(filename + ".csi"):
+                    subprocess.Popen(['bcftools','index', filename],stdout=subprocess.PIPE)
+                self._reader = subprocess.Popen(['bcftools','view', filename],stdout=subprocess.PIPE).stdout
+            else:
+                compressed = compressed or filename.endswith('.gz')
+                self._reader = open(filename, 'rb' if compressed else 'rt')
         self.filename = filename
         if compressed:
             self._reader = gzip.GzipFile(fileobj=self._reader)
@@ -617,13 +625,18 @@ def fetch(self, chrom, start=None, end=None):
         if not self.filename:
             raise Exception('Please provide a filename (or a "normal" fsock)')
 
-        if not self._tabix:
+        if not self._tabix and self.filename.endswith(".bcf") == False:
             self._tabix = pysam.Tabixfile(self.filename)
 
         if self._prepend_chr and chrom[:3] == 'chr':
             chrom = chrom[3:]
 
-        self.reader = self._tabix.fetch(chrom, start, end)
+        if not self.filename.endswith(".bcf"):
+            self.reader = self._tabix.fetch(chrom, start, end)
+        else:
+            print locals()
+            region = "{chrom}:{start}-{end}".format(**locals())
+            self.reader = subprocess.Popen(['bcftools','view','-H', self.filename, region],stdout=subprocess.PIPE).stdout
         return self
 
 
@@ -758,6 +771,23 @@ def __update_readme():
     file('README.rst', 'w').write(vcf.__doc__)
 
 
+def which(program):
+    """ returns the path to an executable or None if it can't be found
+     http://stackoverflow.com/questions/377017/test-if-executable-exists-in-python
+     """
+    def is_exe(fpath):
+        return os.path.isfile(fpath) and os.access(fpath, os.X_OK)
+    fpath, fname = os.path.split(program)
+    if fpath:
+        if is_exe(program):
+            return program
+    else:
+        for path in os.environ["PATH"].split(os.pathsep):
+            exe_file = os.path.join(path, program)
+            if is_exe(exe_file):
+                return exe_file
+    return None
+
 # backwards compatibility
 VCFReader = Reader
 VCFWriter = Writer