Ability to call directly from signal file

.gitignore

@@ -0,0 +1,7 @@
+.venv/
+scratch/
+__pycache__/
+bin/
+*h5
+inv_cdf*
+fseq2/idr_2_0_3/inv_cdf.c

README.md

@@ -91,8 +91,8 @@ Output directory path. Default is current directory.
 Prefix for all output files. This overrides exisiting files. Default is `fseq2_result`.
 
 ##### `-sig_format`
-Signal format for reconstructed signal. Available format `wig`, `bigwig`, `np_array`. Note if choose `np_array`, arrays 
-for each chrom are stored in [`NAME_sig.h5`](#name_sigh5) with `chrom` as key, and no gaussian smooth applied. Default is False, without output signal.
+Signal format for reconstructed signal. Available format `wig`, `bigwig`, `np_array` `memmap_np`. Note if choose `np_array`, arrays 
+for each chrom are stored in [`NAME_sig.h5`](#name_sigh5) with `chrom` as key, and no gaussian smooth applied, and peak calling parameters as attributes. If `memmap_np`, signal is stored in a directory with one file per chromosome. Peaks can be recalled from either `memmap_np` or `np_array` using `callpeak_sig`. Default is False, without output signal.
 
 ##### `-sort_by`
 Sort peaks and summits by `pValue` or `chromAndStart`. Default is `chromAndStart`.
@@ -162,6 +162,26 @@ Useful for IDR analysis: in `callpeak_idr`, we set it to the minimum distance be
 Maximum number of peaks called. Default is not set. If set, overrides `p_thr` and `q_thr`.
 
 
+## `callpeak_sig`
+Call peaks from the signal file. Using an equivalent p/q threshold, these results should exactly match those from `callpeak`.
+#### Command line input:
+
+##### `signal_file`
+The signal file generated by `callpeak` when using `-sig_format np_array` or `-sig_format memmap_np`. Note if using multiple cpus, the `np_array` must be read wholey into memory first, while the `memmap_np` can be read in one chromosome at a time.
+
+##### `-num_peaks`
+Maximum number of peaks called. Default is not set. If set, overrides `p_thr` and `q_thr`.
+
+##### `-p_thr`
+P value threshold. Default is 0.01. Consider to relax it to 0.05 when without control data or calling broad peaks. 
+To resemble F-Seq1 results, specify `-p_thr False`, then filter out peaks whose signalValue 
+(7th column in `.narrowPeak`) below est. threshold.
+
+##### `-q_thr`
+Q value (FDR) threshold. Default is not set and use `p_thr`. If set, only use `q_thr`.
+
+##### `-cpus`
+Number of cores to use. Default is 1.
 
 ## `callpeak_idr`
 #### Command line input:

fseq2/callpeak_main.py

@@ -6,6 +6,7 @@
 import multiprocessing as mp
 import sys
 
+import h5py
 from numpy import ptp
 from pandas import concat, DataFrame
 
@@ -25,6 +26,7 @@ def main(args):
         chrom_ls_to_process = list(chr_size_dic)
     else:
         chrom_ls_to_process = False
+        chr_size_dic = {}
         if (args.sig_format == 'bigwig'):
             sys.exit('Error: please specify chrom_size_file to output bigwig signal.')
 
@@ -109,6 +111,10 @@ def main(args):
         gaussian_smooth_sigma = int(1000 / 6)
         if args.sig_format == 'np_array':
             treatment_np_tmp_name = f'{args.o}/{args.name}_sig.h5'
+        elif args.sig_format == 'memmap_np':
+            # Use temp HDF5 during processing, then convert to individual .npy files
+            treatment_np_tmp_name = f'{args.temp_dir_name}/{args.name}_sig.h5'
+            memmap_output_dir = f'{args.o}/{args.name}_signal_tracks'
         elif (args.sig_format == 'wig') or (args.sig_format == 'bigwig'):
             treatment_np_tmp_name = f'{args.temp_dir_name}/{args.name}_sig.h5'
 
@@ -208,9 +214,40 @@ def main(args):
                             standard_narrowpeak=args.standard_narrowpeak)
 
     if args.sig_format:
+        # Save all parameters to HDF5 file for exact reproducibility with callpeak_sig
+        if args.sig_format in ('np_array', 'memmap_np'):
+            with h5py.File(treatment_np_tmp_name, mode='a', libver='latest') as sig_file:
+                # Core peak calling parameters
+                sig_file.attrs['threshold'] = threshold
+                sig_file.attrs['peak_region_threshold'] = peak_region_threshold
+                sig_file.attrs['min_distance'] = min_distance
+                sig_file.attrs['min_prominence'] = min_prominence
+                sig_file.attrs['window_size'] = window_size
+                sig_file.attrs['lambda_bg_lower_bound'] = lambda_bg_lower_bound
+                # KDE parameters
+                sig_file.attrs['bandwidth'] = bandwidth
+                sig_file.attrs['fragment_offset'] = fragment_offset
+                sig_file.attrs['scaling_factor'] = scaling_factor
+                sig_file.attrs['fragment_size'] = fragment_size
+                sig_file.attrs['feature_length'] = feature_length
+                # Data parameters
+                sig_file.attrs['ncuts'] = ncuts
+                sig_file.attrs['sparse_data'] = args.sparse_data
+                # Control-related parameters (if applicable)
+                if args.control_file:
+                    sig_file.attrs['has_control'] = True
+                    sig_file.attrs['bandwidth_control'] = bandwidth_control
+                    sig_file.attrs['fragment_offset_control'] = fragment_offset_control
+                    sig_file.attrs['ncuts_control'] = ncuts_control
+                else:
+                    sig_file.attrs['has_control'] = False
+
+        # Pass memmap output directory if using memmap_np format
+        memmap_out_dir = memmap_output_dir if args.sig_format == 'memmap_np' else None
         fseq2.output_sig(sig_format=args.sig_format, treatment_np_tmp_name=treatment_np_tmp_name,
                    out_dir=args.o, out_name=args.name, chr_size_dic=chr_size_dic,
-                   sig_float_precision=sig_float_precision, gaussian_smooth_sigma=gaussian_smooth_sigma) #note if output np_array, then no gaussian smooth
+                   sig_float_precision=sig_float_precision, gaussian_smooth_sigma=gaussian_smooth_sigma,
+                   memmap_out_dir=memmap_out_dir)  # note if output np_array, then no gaussian smooth
 
     if args.v:
         print('#3: Done', flush=True)