soft-masked blocks

Author: dpryan79

py2bit.c

@@ -356,6 +356,7 @@ static PyObject *py2bitHardMaskedBlocks(pyTwoBit_t *self, PyObject *args, PyObje
     // Form the output
     ret = PyList_New(totalBlocks);
     if(!ret) goto error;
+    if(totalBlocks == 0) return ret;
     for(i=0, j=0; i<tb->idx->nBlockCount[tid]; i++) {
         blockStart = tb->idx->nBlockStart[tid][i];
         blockEnd = blockStart + tb->idx->nBlockSizes[tid][i];
@@ -375,6 +376,80 @@ static PyObject *py2bitHardMaskedBlocks(pyTwoBit_t *self, PyObject *args, PyObje
     return NULL;
 }
 
+static PyObject *py2bitSoftMaskedBlocks(pyTwoBit_t *self, PyObject *args, PyObject *kwds) {
+    PyObject *ret = NULL, *tup = NULL;
+    TwoBit *tb = self->tb;
+    char *chrom;
+    unsigned long startl = 0, endl = 0, totalBlocks = 0, tid;
+    uint32_t start, end, len, blockStart, blockEnd, i, j;
+    static char *kwd_list[] = {"chrom", "start", "end", NULL};
+
+    if(!tb) {
+        PyErr_SetString(PyExc_RuntimeError, "The 2bit file handle is not open!");
+        return NULL;
+    }
+
+    if(!PyArg_ParseTupleAndKeywords(args, kwds, "s|kk", kwd_list, &chrom, &startl, &endl)) {
+        PyErr_SetString(PyExc_RuntimeError, "You must supply at least a chromosome!");
+        return NULL;
+    }
+
+    //Get the chromosome ID
+    for(i=0; i<tb->hdr->nChroms; i++) {
+        if(strcmp(tb->cl->chrom[i], chrom) == 0) {
+            tid = i;
+            break;
+        }
+    }
+
+    len = twobitChromLen(tb, chrom);
+    if(len == 0) {
+        PyErr_SetString(PyExc_RuntimeError, "The specified chromosome doesn't exist in the 2bit file!");
+        return NULL;
+    }
+    if(endl == 0) endl = len;
+    if(endl > len) endl = len;
+    end = (uint32_t) endl;
+    if(startl >= endl && startl > 0) {
+        PyErr_SetString(PyExc_RuntimeError, "The start value must be less then the end value (and the end of the chromosome");
+        return NULL;
+    }
+    start = (uint32_t) startl;
+
+    if(!tb->idx->maskBlockStart) {
+        PyErr_SetString(PyExc_RuntimeError, "The file was not opened with storeMasked=True! Consequently, there are no stored soft-masked regions.");
+        return NULL;
+    }
+    
+    // Count the total number of overlapping soft-masked blocks
+    for(i=0; i<tb->idx->maskBlockCount[tid]; i++) {
+        blockStart = tb->idx->maskBlockStart[tid][i];
+        blockEnd = blockStart + tb->idx->maskBlockSizes[tid][i];
+        if(blockStart < end && blockEnd > start) totalBlocks++;
+    }
+
+    // Form the output
+    ret = PyList_New(totalBlocks);
+    if(!ret) goto error;
+    if(totalBlocks == 0) return ret;
+    for(i=0, j=0; i<tb->idx->maskBlockCount[tid]; i++) {
+        blockStart = tb->idx->maskBlockStart[tid][i];
+        blockEnd = blockStart + tb->idx->maskBlockSizes[tid][i];
+        if(blockStart < end && blockEnd > start) {
+            tup = Py_BuildValue("(kk)", (unsigned long) blockStart, (unsigned long) blockEnd);
+            if(!tup) goto error;
+            if(PyList_SetItem(ret, j++, tup)) goto error;
+        }
+    }
+
+    return ret;
+error:
+    if(ret) Py_XDECREF(ret);
+    if(tup) Py_XDECREF(tup);
+    PyErr_SetString(PyExc_RuntimeError, "Received an error while constructing the output list and tuples!");
+    return NULL;
+}
+
 #if PY_MAJOR_VERSION >= 3
 PyMODINIT_FUNC PyInit_py2bit(void) {
     PyObject *res;

py2bit.h

@@ -17,6 +17,7 @@ static PyObject* py2bitChroms(pyTwoBit_t *pybw, PyObject *args);
 static PyObject *py2bitSequence(pyTwoBit_t *pybw, PyObject *args, PyObject *kwds);
 static PyObject *py2bitBases(pyTwoBit_t *pybw, PyObject *args, PyObject *kwds);
 static PyObject *py2bitHardMaskedBlocks(pyTwoBit_t *pybw, PyObject *args, PyObject *kwds);
+static PyObject *py2bitSoftMaskedBlocks(pyTwoBit_t *pybw, PyObject *args, PyObject *kwds);
 static void py2bitDealloc(pyTwoBit_t *pybw);
 
 static PyMethodDef tbMethods[] = {
@@ -163,7 +164,28 @@ Returns:\n\
 >>> import py2bit\n\
 >>> tb = py2bit.open(\"test/test.2bit\")\n\
 >>> print(tb.hardMaskedBlocks(\"chr1\")\n\
->>> \n\
+[(0, 50), (100, 150)]\n\
+>>> print(tb.hardMaskedBlocks(\"chr1\", 75, 100)\n\
+[]\n\
+>>> print(tb.hardMaskedBlocks(\"chr1\", 75, 101)\n\
+[(100, 150)]\n\
+>>> tb.close()"},
+    {"softMaskedBlocks", (PyCFunction)py2bitSoftMaskedBlocks, METH_VARARGS|METH_KEYWORDS,
+"Retrieve a list of soft-masked blocks on a single-chromosome (or range on it).\n\
+\n\
+Positional arguments:\n\
+    chr:   Chromosome name\n\
+\n\
+Optional keyword arguments:\n\
+    start: Starting position (0-based)\n\
+    end:   Ending position (1-based)\n\
+\n\
+Returns:\n\
+    A list of tuples, with items start and end.\n\
+\n\
+>>> import py2bit\n\
+>>> tb = py2bit.open(\"test/test.2bit\", storeMasked=True)\n\
+>>> print(tb.softMaskedBlocks(\"chr1\")\n\
 >>> tb.close()"},
     {"__enter__", (PyCFunction) py2bitEnter, METH_NOARGS, NULL},
     {"__exit__", (PyCFunction) py2bitClose, METH_VARARGS, NULL},

py2bitTest/test.py

@@ -56,3 +56,9 @@ def testHardMaskedBlocks(self):
         assert(tb.hardMaskedBlocks("chr1", 75, 101) == [(100, 150)])
         assert(tb.hardMaskedBlocks("chr2") == [(50, 100)])
         tb.close()
+
+    def testSoftMaskedBlocks(self):
+        tb = py2bit.open(self.fname, storeMasked=True)
+        assert(tb.softMaskedBlocks("chr1") == [(62, 70)])
+        assert(tb.softMaskedBlocks("chr1", 0, 50) == [])
+        tb.close()