Merge pull request #113 from dpryan79/fix99

.azure-pipelines/test-template.yml

@@ -1,10 +1,11 @@
 steps:
 - bash: |
-    source activate foo
+    source activate MethylDackel
+    export CONDA_BUILD_SYSROOT=''
     make install CC=$CC CFLAGS="-O3 -Wall -I$CONDA_PREFIX/include" LIBS="-L$CONDA_PREFIX/lib" LIBBIGWIG="$CONDA_PREFIX/lib/libBigWig.a" prefix=$CONDA_PREFIX/bin
   displayName: Installing methylDackel
 - bash: |
-    source activate foo
+    source activate MethylDackel
     export LD_LIBRARY_PATH=$CONDA_PREFIX/lib
-    make test CC=$CC CFLAGS="-O3 -Wall -I$CONDA_PREFIX/include" LIBS="-L$CONDA_PREFIX/lib" LIBBIGWIG="$CONDA_PREFIX/lib/libBigWig.a"
+    python tests/test.py
   displayName: Test methylDackel

CHANGELOG.md

@@ -1,3 +1,7 @@
+Version 0.5.2:
+
+   * Rewrote how read-pair overlap handling is performed. It now uses the constructor/destructor mechanism from htslib instead of using internal htslib structures and functions. This allows supporting newer htslib versions. Currently 1.11 is the only tested and working version, due to changes in the pileup constructor interface in it. (issue #99)
+
 Version 0.5.1:
 
    * Fixed an issue in `MethylDackel mbias` due to an uninitialized value (issue #93).

Makefile

@@ -10,8 +10,8 @@ CFLAGS ?= -Wall -g -O3 -pthread
 
 all: MethylDackel
 
-OBJS = common.o bed.o svg.o pileup.o extract.o MBias.o mergeContext.o perRead.o
-VERSION = 0.5.1
+OBJS = common.o bed.o svg.o overlaps.o extract.o MBias.o mergeContext.o perRead.o
+VERSION = 0.5.2
 
 version.h:
 	echo '#define VERSION "$(VERSION)"' > $@
@@ -28,7 +28,7 @@ lib: libMethylDackel.a
 MethylDackel: libMethylDackel.a version.h $(OBJS)
 	$(CC) $(CFLAGS) $(LIBS) -o MethylDackel $(OBJS) main.c libMethylDackel.a $(LIBBIGWIG) -lm -lz -lpthread -lhts -lcurl
 
-test: MethylDackel 
+test: MethylDackel
 	python tests/test.py
 
 clean:

MethylDackel.h

@@ -5,16 +5,21 @@
 #include "htslib/faidx.h"
 #include "bigWig.h"
 
+#if HTS_VERSION < 101100
+#error "The minimum supported version of htslib is 1.11!"
+#endif
+
+
 //These are needed to handle multiple threads
-pthread_mutex_t positionMutex;
-pthread_mutex_t bwMutex;
-pthread_mutex_t outputMutex;
-uint32_t globalTid;
-uint32_t globalPos;
-uint32_t globalEnd;
-uint32_t bin;
-uint32_t outputBin;
-uint64_t globalnVariantPositions;
+extern pthread_mutex_t positionMutex;
+extern pthread_mutex_t bwMutex;
+extern pthread_mutex_t outputMutex;
+extern uint32_t globalTid;
+extern uint32_t globalPos;
+extern uint32_t globalEnd;
+extern uint32_t bin;
+extern uint32_t outputBin;
+extern uint64_t globalnVariantPositions;
 
 /*! @typedef
  @abstract Structure to hold one region defined in a BED file
@@ -121,13 +126,15 @@ typedef struct {
  @field	config:	The Config* structure containing the settings
  @field hdr:	The input header
  @field iter:	The alignment iterator that should be traversed.
+ @field ohash:  A pointer to the hash table needed for overlap detection
  @field bedIdx: The last index into the BED file
 */
 typedef struct {
     Config *config;
     htsFile *fp;
     bam_hdr_t *hdr;
     hts_itr_t *iter;
+    void *ohash;
     int32_t bedIdx;
 } mplp_data;
 
@@ -224,3 +231,9 @@ void parseBounds(char *s2, int *vals, int mult);
 
 //Used internally to not split CpGs/CHGs between threads
 void adjustBounds(Config *config, bam_hdr_t *hdr, faidx_t *fai, uint32_t *localTid, uint32_t *localPos, uint32_t *localEnd);
+
+// Read-pair overlap handling functions
+int custom_overlap_constructor(void *data, const bam1_t *b, bam_pileup_cd *cd);
+int custom_overlap_destructor(void *data, const bam1_t *b, bam_pileup_cd *cd);
+void *initOlapHash();
+void destroyOlapHash(void *ohash);

README.md

@@ -7,7 +7,7 @@
 Prerequisites
 =============
 
-A C compiler, such as gcc, htslib (versions 1.4 through 1.9 are known to be compatible) and libBigWig are required. For libBigWig, the static library is used.
+A C compiler, such as gcc, htslib (at least versions 1.11, earlier versions are not compatible) and libBigWig are required. For libBigWig, the static library is used.
 
 Compilation
 ===========

azure-pipeline.yml

@@ -8,37 +8,39 @@ jobs:
     vmImage: 'ubuntu-16.04'
   strategy:
     matrix:
-      Python37:
-        python.version: '3.7'
+      htslib111:
+        htslib_version: '1.11'
     maxParallel: 1
 
   steps:
   - bash: echo "##vso[task.prependpath]/usr/share/miniconda/bin"
     displayName: Add conda to PATH
-  - bash: conda create -n foo -q --yes -c conda-forge -c bioconda htslib python gcc_impl_linux-64 binutils_impl_linux-64 binutils_linux-64 gcc_linux-64 libgcc-ng libbigwig
+  - bash: conda create -n MethylDackel -q --yes -c conda-forge -c bioconda htslib=$(HTSLIB_VERSION) python gcc_impl_linux-64 binutils_impl_linux-64 binutils_linux-64 gcc_linux-64 libgcc-ng libbigwig
     displayName: Installing dependencies
   - template: .azure-pipelines/test-template.yml
 
 #This tries to link against 32-bit libraries for some reason
-#- job: 'OSX'
-#  pool:
-#    vmImage: 'macOS-10.13'
-#  strategy:
-#    matrix:
-#      Python37:
-#        python.version: '3.7'
-#    maxParallel: 1
-#
-#  steps:
-#  - script: |
-#      echo "Removing homebrew from Azure to avoid conflicts."
-#      curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/uninstall > ~/uninstall_homebrew
-#      chmod +x ~/uninstall_homebrew
-#      ~/uninstall_homebrew -fq
-#      rm ~/uninstall_homebrew
-#    displayName: Remove homebrew
-#  - bash: echo "##vso[task.prependpath]$CONDA/bin"
-#    displayName: Add conda to PATH
-#  - bash: conda create -n foo -q --yes -c conda-forge -c bioconda htslib python clang clang_osx-64 clangxx cctools compiler-rt ld64
-#    displayName: Installing dependencies
-#  - template: .azure-pipelines/test-template.yml
+- job: 'OSX'
+  pool:
+    vmImage: 'macOS-latest'
+  strategy:
+    matrix:
+      htslib111:
+        htslib_version: '1.11'
+    maxParallel: 1
+
+  steps:
+  - script: |
+      echo "Removing homebrew from Azure to avoid conflicts."
+      curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/uninstall > ~/uninstall_homebrew
+      chmod +x ~/uninstall_homebrew
+      ~/uninstall_homebrew -fq
+      rm ~/uninstall_homebrew
+    displayName: Remove homebrew
+  - bash: echo "##vso[task.prependpath]$CONDA/bin"
+    displayName: Add conda to PATH
+  - bash: sudo chown -R 501:20 /usr/local/miniconda/pkgs
+    displayName: Fix permissions
+  - bash: conda create -n MethylDackel -q --yes -c conda-forge -c bioconda htslib=$(HTSLIB_VERSION) python clang clang_osx-64 clangxx cctools compiler-rt ld64 libbigwig
+    displayName: Installing dependencies
+  - template: .azure-pipelines/test-template.yml

extract.c

@@ -28,7 +28,6 @@ struct lastCall{
     uint32_t nmethyl, nunmethyl;
 };
 
-const char *quux[2] = {"foo", "bar"};
 const char *TriNucleotideContexts[25] = {"CAA", "CAC", "CAG", "CAT", "CAN", \
                                          "CCA", "CCC", "CCG", "CCT", "CCN", \
                                          "CGA", "CGC", "CGG", "CGT", "CGN", \
@@ -247,7 +246,8 @@ void *extractCalls(void *foo) {
     Config *config = (Config*) foo;
     bam_hdr_t *hdr;
     bam_mplp_t iter;
-    int ret, tid, pos, i, seqlen, type, rv, o = 0;
+    int ret, tid, i, seqlen, type, rv, o = 0;
+    hts_pos_t pos;
     int32_t bedIdx = 0;
     int n_plp; //This will need to be modified for multiple input files
     int strand, direction;
@@ -314,7 +314,7 @@ void *extractCalls(void *foo) {
     data->fp = fp;
     data->bedIdx = bedIdx;
 
-    plp = calloc(1, sizeof(bam_pileup1_t *)); //This will have to be modified for multiple input files
+    plp = calloc(1, sizeof(bam_pileup1_t *));
     if(plp == NULL) {
         fprintf(stderr, "Couldn't allocate space for the plp structure in extractCalls()!\n");
         return NULL;
@@ -385,10 +385,13 @@ void *extractCalls(void *foo) {
         }
 
         //Start the pileup
+        data->ohash = initOlapHash();
         iter = bam_mplp_init(1, filter_func, (void **) &data);
-        bam_mplp_init_overlaps(iter);
         bam_mplp_set_maxcnt(iter, INT_MAX);
-        while((ret = cust_mplp_auto(iter, &tid, &pos, &n_plp, plp)) > 0) {
+        bam_mplp_constructor(iter, custom_overlap_constructor);
+        bam_mplp_destructor(iter, custom_overlap_destructor);
+
+        while((ret = bam_mplp64_auto(iter, &tid, &pos, &n_plp, plp)) > 0) {
             if(pos < localPos || pos >= localEnd) continue; // out of the region requested
 
             if(config->bed) { //Handle -l
@@ -525,6 +528,7 @@ void *extractCalls(void *foo) {
             pthread_mutex_unlock(&outputMutex);
             break;
         }
+        destroyOlapHash(data->ohash);
     }
 
     free(os_CpG->s); free(os_CpG);
@@ -567,7 +571,7 @@ void extract_usage() {
 " -q INT           Minimum MAPQ threshold to include an alignment (default 10)\n"
 " -p INT           Minimum Phred threshold to include a base (default 5). This\n"
 "                  must be >0.\n"
-" -D INT           Maximum per-base depth (default 2000)\n"
+" -D INT           Ignored, kept only for backward compatibility.\n"
 " -d INT           Minimum per-base depth for reporting output. If you use\n"
 "                  --mergeContext, this then applies to the merged CpG/CHG.\n"
 "                  (default 1)\n"

main.c

@@ -1,8 +1,19 @@
 #include <string.h>
 #include <stdio.h>
+#include <pthread.h>
 #include "htslib/hts.h"
 #include "version.h" //This has the VERSION define
 
+pthread_mutex_t positionMutex;
+pthread_mutex_t bwMutex;
+pthread_mutex_t outputMutex;
+uint32_t globalTid = 0;
+uint32_t globalPos = 0;
+uint32_t globalEnd = 0;
+uint32_t bin = 0;
+uint32_t outputBin = 0;
+uint64_t globalnVariantPositions = 0;
+
 int mbias_main(int argc, char *argv[]);
 int extract_main(int argc, char *argv[]);
 int mergeContext_main(int argc, char *argv[]);