Add options to apply SNP filter on allc files using filter-allc; add multiprocess option in filter-allc

Author: yupenghe

methylpy/parser.py

@@ -246,11 +246,14 @@ def parse_args():
                                 no_reindex=args.no_reindex)
      elif args.command == "filter-allc":
           from methylpy.utilities import filter_allc_file
-          filter_allc_file(allc_file=args.allc_file,
-                           output_file=args.output_file,
+          filter_allc_file(allc_files=args.allc_files,
+                           output_files=args.output_files,
+                           num_procs=args.num_procs,
                            mc_type=args.mc_type,
                            chroms=args.chroms,
                            compress_output=args.compress_output,
+                           max_mismatch=args.max_mismatch,
+                           max_mismatch_frac=args.max_mismatch_frac,
                            min_cov=args.min_cov)
 
      elif args.command == "allc-to-bigwig":
@@ -1480,22 +1483,29 @@ def add_filter_allc_subparser(subparsers):
                                      help="Filter allc file")
 
      filter_allc_req = filter_allc.add_argument_group("required inputs")
-     filter_allc_req.add_argument("--allc-file",
-                             type=str,
-                             required=True,
-                             help="Allc file to filter.")
+     filter_allc_req.add_argument("--allc-files",
+                                  type=str,
+                                  required=True,
+                                  nargs="+",
+                                  help="allc files to filter.")
+
+     filter_allc_req.add_argument("--output-files",
+                                  type=str,
+                                  required=True,
+                                  nargs="+",
+                                  help="Name of output files. Each output file matches each allc file.")
+
+     filter_allc_req.add_argument("--num-procs",
+                                  type=int,
+                                  default=1,
+                                  help="Number of processors you wish to use to parallelize this function")
 
-     filter_allc_req.add_argument("--output-file",
-                             type=str,
-                             required=True,
-                             help="Name of output file.")
-     
      filter_allc_opt = filter_allc.add_argument_group("optional inputs")
 
      filter_allc_opt.add_argument("--mc-type",
                                   type=str,
                                   nargs="+",
-                                  default=["CGN"],
+                                  default=None,
                                   help="List of space separated cytosine nucleotide contexts for "
                                   + "sites to be included in output file. These classifications "
                                   + "may use the wildcards H (indicating anything but a G) and "
@@ -1507,6 +1517,28 @@ def add_filter_allc_subparser(subparsers):
                                   help="Minimum number of reads that must cover a site for it to be "
                                   + "included in output file.")
 
+     filter_allc_opt.add_argument("--max-mismatch",
+                                  type=int,
+                                  nargs="+",
+                                  default=None,
+                                  help="Maximum numbers of mismatch basecalls allowed in each nucleotide in "
+                                  +"the sequence context of a site for it to be included in output file. If "
+                                  +"the sequence context has three nucleotides, an example of this option is \"0 1 2\". "
+                                  +"It requires no mismatch basecall at the first nucleotide, at most one mismatch "
+                                  +"basecall at the second nucleotide, and at most two at the third nucleotide for a site "
+                                  +"to be reported.")
+
+     filter_allc_opt.add_argument("--max-mismatch-frac",
+                                  type=float,
+                                  nargs="+",
+                                  default=None,
+                                  help="Maximum fraction of mismatch basecalls out of unambiguous basecalls allowed "
+                                  +"in each nucleotide in the sequence context of a site for it to be included "
+                                  +" in output file. If the sequence context has three nucleotides, an example "
+                                  +"of this option is \"0 0 0.1\". It requires no mismatch basecall at the first "
+                                  +"and second nucleotide, and at most 10% mismatches out of unambiguous basecalls "
+                                  +"at the third nucleotide for a site to be reported.")
+
      filter_allc_opt.add_argument("--compress-output",
                                   type=str2bool,
                                   default=True,

methylpy/utilities.py

@@ -191,15 +191,68 @@ def convert_allc_to_bigwig(input_allc_file,
     subprocess.check_call(shlex.split("rm "+output_file+".wig "+output_file+".chrom_size"))
 
 
-def filter_allc_file(allc_file,
-                     output_file,
-                     mc_type="CGN",
+def filter_allc_file(allc_files,
+                     output_files,
+                     num_procs=1,
+                     mc_type=None,
                      chroms = None,
                      compress_output=False,
                      min_cov=0,
+                     max_mismatch=None,
+                     max_mismatch_frac=None,
                      buffer_line_number=100000):
 
-    mc_class = expand_nucleotide_code(mc_type)
+    # User input checks
+    if not isinstance(allc_files, list):
+        exit("allc files must be a list of string(s)")
+    if not isinstance(output_files, list):
+        exit("output files must be a list of string(s)")
+    if len(allc_files) != len(output_files):
+        exit("Number of allc files does not match number of output files")
+
+    if num_procs > 1:
+        pool = multiprocessing.Pool(min(num_procs,len(allc_files)))
+        print_checkpoint("Filtering allc files using "
+                         +str(min(num_procs,len(allc_files)))
+                         +" node(s).")
+        for allc_file,output_file in zip(allc_files,output_files):
+            pool.apply_async(filter_allc_file_worker,
+                             (),
+                             {"allc_file":allc_file,
+                              "output_file":output_file,
+                              "mc_type":mc_type,
+                              "chroms":chroms,
+                              "compress_output":compress_output,
+                              "min_cov":min_cov,
+                              "max_mismatch":max_mismatch,
+                              "max_mismatch_frac":max_mismatch_frac
+                             })
+        pool.close()
+        pool.join()
+    else:
+        print_checkpoint("Filtering allc files using single node.")
+        for allc_file,output_file in zip(allc_files,output_files):
+            filter_allc_file_worker(allc_file=allc_file,
+                                    output_file=output_file,
+                                    mc_type=mc_type,
+                                    chroms=chroms,
+                                    compress_output=compress_output,
+                                    min_cov=min_cov,
+                                    max_mismatch=max_mismatch,
+                                    max_mismatch_frac=max_mismatch_frac)
+    
+def filter_allc_file_worker(allc_file,
+                            output_file,
+                            mc_type=None,
+                            chroms = None,
+                            compress_output=False,
+                            min_cov=0,
+                            max_mismatch=None,
+                            max_mismatch_frac=None,
+                            buffer_line_number=100000):
+
+    if mc_type is not None:
+        mc_class = expand_nucleotide_code(mc_type)
     # input & output
     f = open_allc_file(allc_file)
     if compress_output:
@@ -214,14 +267,48 @@ def filter_allc_file(allc_file,
     for line in f:
         fields = line.split("\t")
         if (chroms is None or fields[0] in chroms) \
-           and fields[3] in mc_class \
+           and (mc_type is None or fields[3] in mc_class) \
            and int(fields[5]) >= min_cov:
-            line_counts += 1
-            out += line
+            pass
+        else:
+            continue
+        
+        if max_mismatch is not None or max_mismatch_frac is not None:
+            try:
+                matches = map(int,fields[7].split(","))
+                mismatches = map(int,fields[8].split(","))
+            except:
+                print_error("allc file "+allc_file+"does not contain SNP information used "
+                            +"for applying mismatch-based filtering!\n")
+            pass_snp_filter = True
+            if max_mismatch is not None:
+                try:
+                    for index,cutoff in enumerate(max_mismatch):
+                        if mismatches[index] > cutoff:
+                            pass_snp_filter = False
+                            break
+                except:
+                    pass
+            if max_mismatch_frac is not None:
+                try:
+                    for index,cutoff in enumerate(max_mismatch_frac):
+                        cov = mismatches[index] + matches[index]
+                        if cov > 0 and float(mismatches[index])/float(cov) > cutoff:
+                            pass_snp_filter = False
+                            break
+                except:
+                    pass
+            if not pass_snp_filter:
+                continue
+            
+        line_counts += 1
+        out += line
+        # print out once reach buffer limit
         if line_counts > buffer_line_number:
             output_fhandler.write(out)
             line_counts = 0
             out = ""
+    # clear buffer
     if line_counts > 0:
         output_fhandler.write(out)
         out = ""
@@ -382,16 +469,17 @@ def merge_allc_files_worker(allc_files,
     fhandles = []
     chrom_pointer = []
     chroms = set([])
-    try:
+    #try:
+    if True:
         for index,allc_file in enumerate(allc_files):
             fhandles.append(open_allc_file(allc_file))
             chrom_pointer.append(read_allc_index(allc_file))
             for chrom in chrom_pointer[index].keys():
                 chroms.add(chrom)
-    except:
-        for f in fhandles:
-            f.close()
-        exit() # exit due to failure of openning all allc files at once
+   # except:
+   #     for f in fhandles:
+   #         f.close()
+   #     exit() # exit due to failure of openning all allc files at once
     if query_chroms is not None:
         if isinstance(query_chroms,list):
             chroms = query_chroms