Merge master

CHANGELOG.md

@@ -2,6 +2,7 @@ Version 0.2.0:
 
    * Changed the package name from PileOMeth to MethylDackel. It is unfortunate that the temporary "PileOMeth" name stuck around for so long.
    * Fixed the plotting where the sometimes the plotted lines escape the bounds of the graph. The cause for this was that read #2 was being ignored when the graph bounds were being computed.
+   * Added the `--requireFlags`/`-r` option, which is equivalent to the -f option in samtools. The default is 0, which requires nothing.
 
 Version 0.1.13:
 

MBias.c

@@ -202,6 +202,10 @@ void mbias_usage() {
 "                  0xF00 or 3840 in decimal. If you would like to change that,\n"
 "                  you can specify a new value here.\n"
 "                  ignored. Specifying this causes them to be included.\n"
+" -R, --requireFlags   Require each alignment to have all bits in this value\n"
+"                  present, or else the alignment is ignored. This is equivalent\n"
+"                  to the -f option in samtools. The default is 0, which\n"
+"                  includes all alignments.\n"
 " --txt            Output tab separated metrics to the screen. These can be\n"
 "                  imported into R or another program for manual plotting and\n"
 "                  analysis.\n"
@@ -230,6 +234,7 @@ int mbias_main(int argc, char *argv[]) {
     config.bedName = NULL;
     config.bed = NULL;
     config.ignoreFlags = 0xF00;
+    config.requireFlags = 0;
     for(i=0; i<16; i++) config.bounds[i] = 0;
     for(i=0; i<16; i++) config.absoluteBounds[i] = 0;
 
@@ -243,6 +248,7 @@ int mbias_main(int argc, char *argv[]) {
         {"txt",          0, NULL,   7},
         {"noSVG",        0, NULL,   8},
         {"ignoreFlags",  1, NULL, 'F'},
+        {"requireFlags", 1, NULL, 'R'},
         {"help",         0, NULL, 'h'},
         {"version",      0, NULL, 'v'},
         {0,              0, NULL,   0}
@@ -292,6 +298,9 @@ int mbias_main(int argc, char *argv[]) {
         case 'F' :
             config.ignoreFlags = atoi(optarg);
             break;
+        case 'R' :
+            config.requireFlags = atoi(optarg);
+            break;
         case 'q' :
             config.minMapq = atoi(optarg);
             break;

Makefile

@@ -11,7 +11,7 @@ OPTS ?= -Wall -g -O3
 all: lib MethylDackel
 
 OBJS = common.o bed.o svg.o pileup.o extract.o MBias.o mergeContext.o
-VERSION = 0.1.13
+VERSION = 0.2.0
 
 #If we're building from a git repo, then append the most recent tag
 ifneq "$(wildcard .git)" ""

MethylDackel.h

@@ -47,6 +47,7 @@ typedef struct {
  @field keepDiscordant	0: Do not include discordantly aligned reads when calculating metrics
  @field	keepSingleton	0: Do not include singletons when calculating metrics
  @field ignoreFlags     Mask that's logically &ed with and ignored if > 0. Defaults to 0xF00.
+ @field requireFlags    Mask that's logically &ed with and ignored if < mask. Defaults to 0, which means ignore.
  @field merge   1: Merge Cs in either a CpG or CHG context into single entries
  @field methylKit       Output in a format compatible with methylKit
  @field output_fp	Output file pointers (to CpG, CHG, and CHH, respectively)
@@ -59,7 +60,7 @@ typedef struct {
 typedef struct {
     int keepCpG, keepCHG, keepCHH;
     int minMapq, minPhred, keepDupes, maxDepth, minDepth;
-    int keepDiscordant, keepSingleton, ignoreFlags;
+    int keepDiscordant, keepSingleton, ignoreFlags, requireFlags;
     int merge, methylKit;
     int fraction, counts, logit;
     FILE **output_fp;

common.c

@@ -182,6 +182,7 @@ int filter_func(void *data, bam1_t *b) {
         if(b->core.tid == -1 || b->core.flag & BAM_FUNMAP) continue; //Unmapped
         if(b->core.qual < ldata->config->minMapq) continue; //-q
         if(b->core.flag & ldata->config->ignoreFlags) continue; //By default: secondary alignments, QC failed, PCR duplicates, and supplemental alignments
+        if(ldata->config->requireFlags && (b->core.flag & ldata->config->requireFlags) != ldata->config->requireFlags) continue;
         if(!ldata->config->keepDupes && b->core.flag & BAM_FDUP) continue;
         p = bam_aux_get(b, "NH");
         if(p != NULL) {

extract.c

@@ -310,6 +310,10 @@ void extract_usage() {
 "                  0xF00 or 3840 in decimal. If you would like to change that,\n"
 "                  you can specify a new value here.\n"
 "                  ignored. Specifying this causes them to be included.\n"
+" -R, --requireFlags   Require each alignment to have all bits in this value\n"
+"                  present, or else the alignment is ignored. This is equivalent\n"
+"                  to the -f option in samtools. The default is 0, which\n"
+"                  includes all alignments.\n"
 " --noCpG          Do not output CpG context methylation metrics\n"
 " --CHG            Output CHG context methylation metrics\n"
 " --CHH            Output CHH context methylation metrics\n"
@@ -374,6 +378,7 @@ int extract_main(int argc, char *argv[]) {
     config.counts = 0;
     config.logit = 0;
     config.ignoreFlags = 0xF00;
+    config.requireFlags = 0;
     for(i=0; i<16; i++) config.bounds[i] = 0;
     for(i=0; i<16; i++) config.absoluteBounds[i] = 0;
 
@@ -400,11 +405,12 @@ int extract_main(int argc, char *argv[]) {
         {"nCTOT",        1, NULL,  15},
         {"nCTOB",        1, NULL,  16},
         {"ignoreFlags",  1, NULL, 'F'},
+        {"requireFlags", 1, NULL, 'R'},
         {"help",         0, NULL, 'h'},
         {"version",      0, NULL, 'v'},
         {0,              0, NULL,   0}
     };
-    while((c = getopt_long(argc, argv, "hvq:p:r:l:o:D:f:c:m:d:F:", lopts,NULL)) >=0){
+    while((c = getopt_long(argc, argv, "hvq:p:r:l:o:D:f:c:m:d:F:R:", lopts,NULL)) >=0){
         switch(c) {
         case 'h' :
             extract_usage();
@@ -482,6 +488,9 @@ int extract_main(int argc, char *argv[]) {
         case 'F' :
             config.ignoreFlags = atoi(optarg);
             break;
+        case 'R':
+            config.requireFlags = atoi(optarg);
+            break;
         case 'q' :
             config.minMapq = atoi(optarg);
             break;

tests/test.py

@@ -61,6 +61,13 @@ def rm(f):
 assert lines == 49
 rm('cg_aln_CpG.bedGraph')
 
+# Check that --requireFlags is working
+check_call('../MethylDackel extract --requireFlags 0xD00 cg100.fa cg_aln.bam -q 2', shell=True)
+assert op.exists('cg_aln_CpG.bedGraph')
+lines = sum(1 for _ in open('cg_aln_CpG.bedGraph'))
+assert lines == 49
+rm('cg_aln_CpG.bedGraph')
+
 # Check absolute trimming bounds
 check_call('../MethylDackel extract --nOT 50,50,40,40 cg100.fa cg_aln.bam -q 2', shell=True)
 assert op.exists('cg_aln_CpG.bedGraph')