Repository: hasindu2008/slow5lib (language: python)

Issue:
cleaning



Discussion / hints:

Fix patch:
diff --git a/.github/workflows/c-cpp.yml b/.github/workflows/c-cpp.yml
index 2bb9e042..256e3d7c 100644
--- a/.github/workflows/c-cpp.yml
+++ b/.github/workflows/c-cpp.yml
@@ -9,7 +9,7 @@ on:
 jobs:
   ubuntu_14:
     name: ubuntu_14
-    runs-on: ubuntu-16.04
+    runs-on: ubuntu-18.04
     container: ubuntu:14.04
     steps:
     - uses: actions/checkout@v2
@@ -23,11 +23,12 @@ jobs:
       run: make test
   ubuntu_16:
     name: Ubuntu 16
-    runs-on: ubuntu-16.04
+    runs-on: ubuntu-18.04
+    container: ubuntu:16.04
     steps:
     - uses: actions/checkout@v2
     - name: install packages
-      run: sudo apt-get update && sudo apt-get install -y zlib1g-dev
+      run: apt-get update && apt-get install -y zlib1g-dev gcc g++ autoconf automake make
     - name: build
       run: make -j8
     - name: test
diff --git a/.github/workflows/python.yml b/.github/workflows/python.yml
index 1185bd1f..a5a1f92c 100644
--- a/.github/workflows/python.yml
+++ b/.github/workflows/python.yml
@@ -9,7 +9,7 @@ on:
 jobs:
   ubuntu_14:
     name: pyslow5 Ubuntu_14
-    runs-on: ubuntu-16.04
+    runs-on: ubuntu-18.04
     container: ubuntu:14.04
     steps:
     - uses: actions/checkout@v2
@@ -23,11 +23,12 @@ jobs:
       run: python3 -m unittest -v python/test.py
   ubuntu_16:
     name: pyslow5 Ubuntu 16
-    runs-on: ubuntu-16.04
+    runs-on: ubuntu-18.04
+    container: ubuntu:16.04
     steps:
     - uses: actions/checkout@v2
     - name: install packages
-      run: sudo apt-get update && sudo apt-get install -y zlib1g-dev python3-pip python3-setuptools && pip3 install cython "numpy>=1.14.5,<1.19.0"
+      run: apt-get update && apt-get install -y lsb-release zlib1g-dev python3-pip python3-setuptools && pip3 install cython "numpy>=1.14.5,<1.19.0"
     - name: env
       run: lsb_release -a &&  ldd --version && gcc --version && g++ --version  && python3 --version
     - name: build
diff --git a/CMakeLists.txt b/CMakeLists.txt
index b0661179..2872395e 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -38,6 +38,7 @@ add_subdirectory(thirdparty/streamvbyte)
 #
 # Build a shared lib
 add_library(slow5 SHARED ${slow5_} ${slow5_idx} ${slow5_misc} ${slow5_press})
+target_link_libraries(slow5 streamvbyte_slow5)
 #
 #just to get rid of Clion warnings
 file(GLOB_RECURSE C_SOURCES RELATIVE ${CMAKE_SOURCE_DIR} "test/*.c")
@@ -46,9 +47,9 @@ add_executable(slow5test test/unit_test_two_rg.c) #run make test to test all tes
 #
 # link with dependencies
 IF (${DEPLOY_PLATFORM} STREQUAL "x86")
-    target_link_libraries(slow5test slow5 streamvbyte -lz -ldl -lm -lpthread -lrt)
-    target_link_libraries(slow5exmp slow5 streamvbyte -lz -ldl -lm -lpthread -lrt -fopenmp)
+    target_link_libraries(slow5test slow5 streamvbyte_slow5 -lz -ldl -lm -lpthread -lrt)
+    target_link_libraries(slow5exmp slow5 streamvbyte_slow5 -lz -ldl -lm -lpthread -lrt -fopenmp)
 ELSE()
-    target_link_libraries(slow5test slow5 -lz -ldl -lm)
-    target_link_libraries(slow5exmp slow5 -lz -ldl -lm)
+    target_link_libraries(slow5test slow5 streamvbyte_slow5 -lz -ldl -lm)
+    target_link_libraries(slow5exmp slow5 streamvbyte_slow5 -lz -ldl -lm)
 ENDIF()
diff --git a/Makefile b/Makefile
index d92f35e8..fa10a117 100644
--- a/Makefile
+++ b/Makefile
@@ -17,7 +17,7 @@ endif
 ifeq ($(zstd_local),)
 else
 CFLAGS		+= -DSLOW5_USE_ZSTD
-CPPFLAGS += -I $(zstd_local)
+CPPFLAGS 	+= -I $(zstd_local)
 endif
 BUILD_DIR	= lib
 
diff --git a/README.md b/README.md
index 18851c8c..7da956cc 100644
--- a/README.md
+++ b/README.md
@@ -8,6 +8,7 @@ SLOW5 is a new file format for storing signal data from Oxford Nanopore Technolo
 
 Full documentation: https://hasindu2008.github.io/slow5lib  
 Pre-print: https://www.biorxiv.org/content/10.1101/2021.06.29.450255v1  
+Publication: https://www.nature.com/articles/s41587-021-01147-4  
 SLOW5 specification: https://hasindu2008.github.io/slow5specs
 
 [![BioConda Install](https://img.shields.io/conda/dn/bioconda/pyslow5.svg?style=flag&label=BioConda%20install)](https://anaconda.org/bioconda/pyslow5)
diff --git a/docs/getting_started.md b/docs/getting_started.md
index f3949031..5481d2ee 100644
--- a/docs/getting_started.md
+++ b/docs/getting_started.md
@@ -6,6 +6,7 @@
 SLOW5 is a new file format for storing signal data from Oxford Nanopore Technologies (ONT) devices. SLOW5 was developed to overcome inherent limitations in the standard FAST5 signal data format that prevent efficient, scalable analysis and cause many headaches for developers. SLOW5 can be encoded in human-readable ASCII format, or a more compact and efficient binary format (BLOW5) - this is analogous to the seminal SAM/BAM format for storing DNA sequence alignments. The BLOW5 binary format supports  *zlib* (DEFLATE) compression, or other compression methods, thereby minimising the data storage footprint while still permitting efficient parallel access. Detailed benchmarking experiments have shown that SLOW5 format is an order of magnitude faster and significantly smaller than FAST5.
 
 Pre-print: https://www.biorxiv.org/content/10.1101/2021.06.29.450255v1  
+Publication: https://www.nature.com/articles/s41587-021-01147-4  
 SLOW5 specification: https://hasindu2008.github.io/slow5specs
 
 
diff --git a/include/slow5/slow5.h b/include/slow5/slow5.h
index 2a36ecf9..5055538c 100644
--- a/include/slow5/slow5.h
+++ b/include/slow5/slow5.h
@@ -296,6 +296,7 @@ struct slow5_file_meta {
     const char *pathname;       ///< file path
     int fd;                     ///< file descriptor
     uint64_t start_rec_offset;  ///< offset (in bytes) of the first SLOW5 record (skipping the SLOW5 header; used for indexing)
+    char *fread_buffer;         ///< buffer for fread
 };
 typedef struct slow5_file_meta slow5_file_meta_t;
 
diff --git a/python/pyslow5.pyx b/python/pyslow5.pyx
index 1e2cb4ac..a5e5fe74 100644
--- a/python/pyslow5.pyx
+++ b/python/pyslow5.pyx
@@ -22,6 +22,7 @@ cdef class Open:
     cdef pyslow5.slow5_file_t *s5
     cdef pyslow5.slow5_rec_t *rec
     cdef pyslow5.slow5_rec_t *read
+    cdef bint index_state
     cdef pyslow5.uint64_t head_len
     cdef pyslow5.uint64_t aux_len
     cdef pyslow5.slow5_aux_type *s5_aux_type
@@ -65,6 +66,9 @@ cdef class Open:
         self.s5 = NULL
         self.rec = NULL
         self.read = NULL
+        self.p = ""
+        self.m = ""
+        self.index_state = False
         self.s5_aux_type = NULL
         self.aux_get_err = 1
         self.aux_get_len = 0
@@ -116,21 +120,22 @@ cdef class Open:
                             -7: ["SLOW5_ERR_NOTFOUND", "read id not found"],
                             }
         p = str.encode(pathname)
-        self.p = p
+        self.p = strdup(p)
         m = str.encode(mode)
-        self.m = m
+        self.m = strdup(m)
         # print binary strings of filepath and mode
         self.logger.debug("FILE: {}, mode: {}".format(self.p.decode(), self.m.decode()))
+        self.logger.debug("FILE: {}, mode: {}".format(self.p, self.m))
         # opens file and creates slow5 object
         self.s5 = pyslow5.slow5_open(self.p, self.m)
         if self.s5 is NULL:
             raise MemoryError()
         # load or create and load index
         self.logger.debug("Number of read_groups: {}".format(self.s5.header.num_read_groups))
-        self.logger.debug("Creating/loading index...")
-        ret = slow5_idx_load(self.s5)
-        if ret != 0:
-            self.logger.warning("slow5_idx_load return not 0: {}: {}".format(ret, self.error_codes[ret]))
+        # self.logger.debug("Creating/loading index...")
+        # ret = slow5_idx_load(self.s5)
+        # if ret != 0:
+        #     self.logger.warning("slow5_idx_load return not 0: {}: {}".format(ret, self.error_codes[ret]))
 
 
     def __init__(self, pathname, mode, DEBUG=0):
@@ -159,10 +164,11 @@ cdef class Open:
         range = d['range']
         offset = d['offset']
         raw_unit = range / digitisation
-        new_raw = []
-        for i in d['signal']:
-            j = (i + offset) * raw_unit
-            new_raw.append(float("{0:.2f}".format(round(j,2))))
+        # new_raw = []
+        # for i in d['signal']:
+        #     j = (i + offset) * raw_unit
+        #     new_raw.append(float("{0:.2f}".format(round(j,2))))
+        new_raw = np.array(raw_unit * (d['signal'] + offset), dtype=np.float32)
         return new_raw
 
 
@@ -173,6 +179,15 @@ cdef class Open:
         aux = str 'name'/'all' or list of names of auxiliary fields added to return dictionary
         returns dic = dictionary of main fields for read_id, with any aux fields added
         '''
+        if not self.index_state:
+            self.logger.debug("FILE: {}, mode: {}".format(self.p.decode(), self.m.decode()))
+            self.logger.debug("FILE: {}, mode: {}".format(self.p, self.m))
+            self.logger.debug("Creating/loading index...")
+            ret = slow5_idx_load(self.s5)
+            if ret != 0:
+                self.logger.warning("slow5_idx_load return not 0: {}: {}".format(ret, self.error_codes[ret]))
+            else:
+                self.index_state = True
         dic = {}
         aux_dic = {}
         ID = str.encode(read_id)
diff --git a/setup.py b/setup.py
index ed0b7214..f5c3a5b7 100644
--- a/setup.py
+++ b/setup.py
@@ -87,7 +87,7 @@ def readme():
     ext_modules=extensions,
     cmdclass= {'build_ext': build_ext},
     classifiers = [
-        'Development Status :: 4 - Beta',
+        'Development Status :: 5 - Production/Stable',
         'License :: OSI Approved :: MIT License',
         'Programming Language :: C',
         'Programming Language :: Cython',
diff --git a/src/slow5.c b/src/slow5.c
index c7a6df03..af0aae69 100644
--- a/src/slow5.c
+++ b/src/slow5.c
@@ -65,6 +65,8 @@ KSORT_INIT(str_slow5, ksstr_t, ks_lt_str)
 #define SLOW5_AUX_ARRAY_CAP_INIT (256) /* Initial capacity for parsing auxiliary array: 2^8 */
 #define SLOW5_AUX_ARRAY_STR_CAP_INIT (1024) /* Initial capacity for storing auxiliary array string: 2^10 */
 
+#define SLOW5_FSTREAM_BUFF_SIZE (131072)  /* buffer size for freads and fwrites */
+
 static inline void slow5_free(struct slow5_file *s5p);
 static int slow5_rec_aux_parse(char *tok, char *read_mem, uint64_t offset, size_t read_size, struct slow5_rec *read, enum slow5_fmt format, struct slow5_aux_meta *aux_meta);
 static inline khash_t(slow5_s2a) *slow5_rec_aux_init(void);
@@ -117,11 +119,29 @@ struct slow5_file *slow5_init(FILE *fp, const char *pathname, enum slow5_fmt for
         slow5_errno = SLOW5_ERR_UNK;
         return NULL;
     }
+
+    char *fread_buff = (char *)calloc(SLOW5_FSTREAM_BUFF_SIZE, sizeof(char));
+    if (!fread_buff) {
+        SLOW5_MALLOC_ERROR();
+        slow5_errno = SLOW5_ERR_MEM;
+        return NULL;
+    }
+
+    if (setvbuf(fp, fread_buff, _IOFBF, SLOW5_FSTREAM_BUFF_SIZE) != 0){
+        SLOW5_WARNING("Could not set a large buffer for file stream of '%s': %s.", pathname, strerror(errno));;
+        free(fread_buff);
+        fread_buff = NULL;
+    }
+    else {
+        SLOW5_LOG_DEBUG("Buffer for file stream of '%s' was set to %d.", pathname, SLOW5_FSTREAM_BUFF_SIZE);
+    }
+
     // TODO Attempt to determine from magic number
 
     slow5_press_method_t method;
     struct slow5_hdr *header = slow5_hdr_init(fp, format, &method);
     if (!header) {
+        free(fread_buff);
         SLOW5_ERROR("Parsing slow5 header of file '%s' failed.", pathname);
         return NULL;
     }
@@ -130,6 +150,7 @@ struct slow5_file *slow5_init(FILE *fp, const char *pathname, enum slow5_fmt for
     if (!s5p) {
         SLOW5_MALLOC_ERROR();
         slow5_hdr_free(header);
+        free(fread_buff);
         slow5_errno = SLOW5_ERR_MEM;
         return NULL;
     }
@@ -137,9 +158,11 @@ struct slow5_file *slow5_init(FILE *fp, const char *pathname, enum slow5_fmt for
     s5p->fp = fp;
     s5p->format = format;
     s5p->header = header;
+    s5p->meta.fread_buffer = fread_buff;
 
     s5p->compress = slow5_press_init(method);
     if (!s5p->compress) {
+        free(fread_buff);
         slow5_hdr_free(header);
         free(s5p);
         return NULL;
@@ -147,6 +170,7 @@ struct slow5_file *slow5_init(FILE *fp, const char *pathname, enum slow5_fmt for
 
     if ((s5p->meta.fd = fileno(fp)) == -1) {
         SLOW5_ERROR("Obtaining file descriptor with fileno() failed: %s.", strerror(errno));
+        free(fread_buff);
         slow5_press_free(s5p->compress);
         slow5_hdr_free(header);
         free(s5p);
@@ -157,6 +181,7 @@ struct slow5_file *slow5_init(FILE *fp, const char *pathname, enum slow5_fmt for
     s5p->meta.pathname = pathname;
     if ((s5p->meta.start_rec_offset = ftello(fp)) == -1) {
         SLOW5_ERROR("Obtaining file offset with ftello() failed: %s.", strerror(errno));
+        free(fread_buff);
         slow5_press_free(s5p->compress);
         slow5_hdr_free(header);
         free(s5p);
@@ -338,6 +363,7 @@ static inline void slow5_free(struct slow5_file *s5p) {
         slow5_press_free(s5p->compress);
         slow5_hdr_free(s5p->header);
         slow5_idx_free(s5p->index);
+        free(s5p->meta.fread_buffer);
         free(s5p);
     }
 }
@@ -787,9 +813,10 @@ void *slow5_hdr_to_mem(struct slow5_hdr *header, enum slow5_fmt format, slow5_pr
                     if (value != NULL) {
                         len_to_cp = strlen(value);
 
-                        //special case for "."
-                        if (strlen(value)==0) {
-                            len_to_cp++;
+                        // special case for SLOW5_ASCII_MISSING_CHAR
+                        bool is_empty = (len_to_cp == 0);
+                        if (is_empty) {
+                            len_to_cp ++;
                         }
 
                         // Realloc if necessary
@@ -799,35 +826,32 @@ void *slow5_hdr_to_mem(struct slow5_hdr *header, enum slow5_fmt format, slow5_pr
                             SLOW5_MALLOC_CHK(mem);
                         }
 
-                        if (strlen(value)==0) { //special case for "."
-                            memcpy(mem + len, ".", len_to_cp);
+                        if (is_empty) { // special case for SLOW5_ASCII_MISSING_CHAR
+                            mem[len] = SLOW5_ASCII_MISSING_CHAR;
                         } else {
                             memcpy(mem + len, value, len_to_cp);
                         }
                         len += len_to_cp;
-                    }
-                    else{
-                        // I added this here - hasindu
+                    } else {
                         // Realloc if necessary
-                        if (len + 1 >= cap) { // +1 for "."
+                        if (len + 1 >= cap) { // +1 for SLOW5_ASCII_MISSING_CHAR
                             cap *= 2;
                             mem = (char *) realloc(mem, cap * sizeof *mem);
                             SLOW5_MALLOC_CHK(mem);
                         }
-                        mem[len] = '.';
+                        mem[len] = SLOW5_ASCII_MISSING_CHAR;
                         ++ len;
                     }
                 } else {
                     // Realloc if necessary
-                    if (len + 2 >= cap) { // +2 for . and SLOW5_SEP_COL_CHAR
+                    if (len + 2 >= cap) { // +2 for SLOW5_ASCII_MISSING_CHAR and SLOW5_SEP_COL_CHAR
                         cap *= 2;
                         mem = (char *) realloc(mem, cap * sizeof *mem);
                         SLOW5_MALLOC_CHK(mem);
                     }
                     mem[len] = SLOW5_SEP_COL_CHAR;
-                    ++ len;
-                    mem[len] = '.';
-                    ++ len;
+                    mem[len + 1] = SLOW5_ASCII_MISSING_CHAR;
+                    len += 2;
                 }
             }
 
@@ -3598,7 +3622,7 @@ void *slow5_rec_to_mem(struct slow5_rec *read, struct slow5_aux_meta *aux_meta,
             for (uint64_t i = 0; i < aux_meta->num; ++ i) {
                 struct slow5_rec_aux_data aux_data = { 0 };
 
-                bool hacky_malloc_flag = 0;
+                bool malloc_flag = false;
 
                 khint_t pos = kh_get(slow5_s2a, read->aux_map, aux_meta->attrs[i]);
                 if (pos != kh_end(read->aux_map)) {
@@ -3607,7 +3631,7 @@ void *slow5_rec_to_mem(struct slow5_rec *read, struct slow5_aux_meta *aux_meta,
                     aux_data.len = 1;
                     aux_data.bytes = SLOW5_AUX_TYPE_META[aux_meta->types[i]].size;
                     aux_data.data = (uint8_t *) malloc(aux_data.bytes);
-                    hacky_malloc_flag = 1;
+                    malloc_flag = true;
                     slow5_memcpy_null_type(aux_data.data, aux_meta->types[i]);
                 }
 
@@ -3633,7 +3657,7 @@ void *slow5_rec_to_mem(struct slow5_rec *read, struct slow5_aux_meta *aux_meta,
                 }
                 curr_len += aux_data.bytes;
 
-                if(hacky_malloc_flag){
+                if (malloc_flag) {
                     free(aux_data.data);
                 }
             }
diff --git a/thirdparty/streamvbyte/CMakeLists.txt b/thirdparty/streamvbyte/CMakeLists.txt
index 7a363deb..407f71e8 100644
--- a/thirdparty/streamvbyte/CMakeLists.txt
+++ b/thirdparty/streamvbyte/CMakeLists.txt
@@ -1,6 +1,6 @@
 cmake_minimum_required(VERSION 3.14)
 
-project(streamvbyte)
+project(streamvbyte_slow5)
 #todo add processor specific flags
 set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -g -Wall -fPIC -march=native -std=c99 -O3 -Wall -Wextra -pedantic -Wshadow")
 
@@ -9,4 +9,4 @@ set(streamvbyte_encode src/streamvbyte_encode.c)
 set(streamvbyte_zigzag src/streamvbyte_zigzag.c)
 
 # Build a shared lib
-add_library(streamvbyte STATIC ${streamvbyte_decode} ${streamvbyte_encode} ${streamvbyte_zigzag})
\ No newline at end of file
+add_library(streamvbyte_slow5 STATIC ${streamvbyte_decode} ${streamvbyte_encode} ${streamvbyte_zigzag})


Test patch:
diff --git a/python/test.py b/python/test.py
index 3508ff77..87e008e0 100644
--- a/python/test.py
+++ b/python/test.py
@@ -69,7 +69,7 @@ def test_len_raw_signal(self):
     def test_pylen_signal(self):
         self.assertEqual(len(self.read['signal']), 59670)
     def test_signal(self):
-        self.assertEqual(sum(self.read['signal'][:10]), sum([190.26, 108.92, 109.46, 109.1, 107.67, 108.39, 108.75, 109.1, 111.07, 108.39]))
+        self.assertEqual(int(sum([round(i, 2) for i in self.read['signal'][:10]])), int(sum([190.26, 108.92, 109.46, 109.1, 107.67, 108.39, 108.75, 109.1, 111.07, 108.39])))
 
 
 class TestAUX(unittest.TestCase):
@@ -215,7 +215,7 @@ def test_seq_reads_pA_aux_all(self):
         for i, read in enumerate(reads):
             with self.subTest(i=i, attr=read['read_id']):
                 self.assertEqual(read['read_id'], results[i][0])
-                self.assertEqual(sum(read['signal'][:10]), results[i][1])
+                self.assertEqual(int(sum([round(i, 2) for i in read['signal'][:10]])), int(results[i][1]))
                 self.assertEqual(read['start_time'] ,results[i][2])