adds a feature to limit by insert size

MethylDackel.h

@@ -89,7 +89,7 @@ typedef struct {
 */
 typedef struct {
  int keepCpG, keepCHG, keepCHH;
- int minMapq, minPhred, keepDupes, minDepth;
+ int minMapq, minPhred, keepDupes, minDepth, minInsertSize, maxInsertSize;
  int keepDiscordant, keepSingleton, ignoreFlags, requireFlags;
  int merge, methylKit, minOppositeDepth;
  int ignoreNH;

extract.c

@@ -580,6 +580,8 @@ void extract_usage() {
 " -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"
+" --minInsertSize INT Minimum insert size to include a read (default 0)\n"
+" --maxInsertSize INT Maximum insert size to include a read (default 2000)\n"
 " -r STR Region string in which to extract methylation\n"
 " -l FILE A BED file listing regions for inclusion.\n"
 " --keepStrand If a BED file is specified, then this option will cause the\n"
@@ -749,6 +751,8 @@ int extract_main(int argc, char *argv[]) {
  config.cytosine_report = 0;
  config.noBAM = 0;
  config.minConversionEfficiency = 0.0;
+ config.minInsertSize = 0;
+ config.maxInsertSize = 2000;
  for(i=0; i<16; i++) config.bounds[i] = 0;
  for(i=0; i<16; i++) config.absoluteBounds[i] = 0;
 
@@ -781,6 +785,8 @@ int extract_main(int argc, char *argv[]) {
  {"cytosine_report", 0, NULL, 21},
  {"minConversionEfficiency", 1, NULL, 22},
  {"ignoreNH", 0, NULL, 23},
+ {"minInsertSize",1, NULL, 24},
+ {"maxInsertSize",1, NULL, 25},
  {"ignoreFlags", 1, NULL, 'F'},
  {"requireFlags", 1, NULL, 'R'},
  {"help", 0, NULL, 'h'},
@@ -893,6 +899,12 @@ int extract_main(int argc, char *argv[]) {
  case 23:
  config.ignoreNH = 1;
  break;
+ case 24:
+ config.minInsertSize = atoi(optarg);
+ break;
+ case 25: 
+ config.maxInsertSize = atoi(optarg);
+ break;
  case 'M':
  config.BWName = optarg;
  break;
@@ -1338,6 +1350,8 @@ int extract_main(int argc, char *argv[]) {
 
  }
 
+
+
 
 
  //Output files

perRead.c

@@ -191,6 +191,7 @@ void *perReadMetrics(void *foo) {
  if(config->requireFlags && (config->requireFlags & b->core.flag) != config->requireFlags) continue;
  if(config->ignoreFlags && (config->ignoreFlags & b->core.flag) != 0) continue;
  if(b->core.qual < config->minMapq) continue;
+ if(llabs(b->core.isize) < config->minInsertSize || llabs(b->core.isize) > config->maxInsertSize) continue;
  processRead(config, b, seq, localPos2, seqlen, &nmethyl, &nunmethyl);
  addRead(os, b, hdr, nmethyl, nunmethyl);
  }

tests/test.py

@@ -143,5 +143,13 @@ def rm(f):
 lines = sum(1 for _ in open('test15_CpG.bedGraph'))
 assert lines == 49
 rm('test15_CpG.bedGraph')
-print("Finished correctly")
 
+# Test min/max InsertSize
+rm('test16_CpG.bedGraph')
+check_call([MPath, 'extract', '-o', 'test16', '--minInsertSize', '0', '--maxInsertSize', '200', '-q', '1', 'cg100.fa', 'cg_aln.bam'])
+assert op.exists('test16_CpG.bedGraph')
+lines = sum(1 for _ in open('test16_CpG.bedGraph'))
+assert lines == 49
+rm('test16_CpG.bedGraph')
+
+print("Finished correctly")