diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..1584a56 --- /dev/null +++ b/.gitignore @@ -0,0 +1,7 @@ +.venv/ +scratch/ +__pycache__/ +bin/ +*h5 +inv_cdf* +fseq2/idr_2_0_3/inv_cdf.c diff --git a/README.md b/README.md index f94f9d1..97a4db7 100644 --- a/README.md +++ b/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: diff --git a/fseq2/callpeak_main.py b/fseq2/callpeak_main.py index cabb576..7316b58 100644 --- a/fseq2/callpeak_main.py +++ b/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) diff --git a/fseq2/callpeak_sig_main.py b/fseq2/callpeak_sig_main.py new file mode 100644 index 0000000..039bf5d --- /dev/null +++ b/fseq2/callpeak_sig_main.py @@ -0,0 +1,361 @@ +#!/usr/bin/env python + +"""F-Seq Version2 call peaks from signal file main script. + +This module provides functionality to call peaks directly from pre-computed +signal files produced by the callpeak command with -sig_format np_array +or -sig_format memmap_np options. +""" + +import json +import multiprocessing as mp +import os +import time + +import h5py +import numpy as np +from pandas import concat, DataFrame + +from fseq2 import fseq2 + + +def main(args): + """Main entry point for callpeak_sig subcommand. + + Args: + args: argparse.Namespace with signal file path and parameters + """ + ### 1. Setup and validation ### + if args.v: + print('=====================================', flush=True) + print(f'F-Seq Version {fseq2.__version__}', flush=True) + print('=====================================', flush=True) + print('#1: Loading signal file and calculating parameters', flush=True) + + signal_file_path = args.signal_file + + # Detect input format: HDF5 file or memmap_np directory + is_memmap_dir = os.path.isdir(signal_file_path) + + if is_memmap_dir: + # memmap_np directory format + metadata_path = os.path.join(signal_file_path, 'metadata.json') + if not os.path.exists(metadata_path): + raise ValueError(f"No metadata.json found in {signal_file_path}. " + "Is this a valid memmap_np signal directory?") + + with open(metadata_path, 'r') as f: + metadata = json.load(f) + + chrom_ls = list(metadata['chromosomes'].keys()) + + if not chrom_ls: + raise ValueError(f"No chromosome data found in {signal_file_path}") + + # Always read parameters from file + params = metadata.get('parameters', {}) + if 'threshold' not in params: + raise ValueError( + "Signal directory does not contain callpeak parameters. " + "Regenerate with a recent version of fseq2 callpeak -sig_format memmap_np." + ) + threshold = float(params['threshold']) + peak_region_threshold = float(params['peak_region_threshold']) + min_prominence = float(params['min_prominence']) + lambda_bg_lower_bound = float(params['lambda_bg_lower_bound']) + sparse_data = bool(params['sparse_data']) + # Allow command-line overrides for min_distance and window_size + min_distance = args.min_distance if args.min_distance is not None else int(params['min_distance']) + window_size = args.window_size if args.window_size is not None else int(params['window_size']) + else: + # HDF5 file format + with h5py.File(signal_file_path, 'r') as sig_file: + chrom_ls = list(sig_file.keys()) + + if not chrom_ls: + raise ValueError(f"No chromosome data found in {signal_file_path}") + + # Always read parameters from file + if 'threshold' not in sig_file.attrs: + raise ValueError( + "Signal file does not contain callpeak parameters. " + "Regenerate with a recent version of fseq2 callpeak -sig_format np_array." + ) + threshold = float(sig_file.attrs['threshold']) + peak_region_threshold = float(sig_file.attrs['peak_region_threshold']) + min_prominence = float(sig_file.attrs['min_prominence']) + lambda_bg_lower_bound = float(sig_file.attrs['lambda_bg_lower_bound']) + sparse_data = bool(sig_file.attrs['sparse_data']) + # Allow command-line overrides for min_distance and window_size + min_distance = args.min_distance if args.min_distance is not None else int(sig_file.attrs['min_distance']) + window_size = args.window_size if args.window_size is not None else int(sig_file.attrs['window_size']) + + if args.v: + format_type = "memmap_np directory" if is_memmap_dir else "HDF5 file" + print(f'\tSignal source: {signal_file_path} ({format_type})', flush=True) + print(f'\tChromosomes found: {len(chrom_ls)}', flush=True) + print(f'\tThreshold: {threshold:.3f}', flush=True) + print(f'\tPeak region threshold: {peak_region_threshold:.3f}', flush=True) + print(f'\tMin distance: {min_distance}', flush=True) + print(f'\tWindow size: {window_size}', flush=True) + print(f'\tMin prominence: {min_prominence:.3f}', flush=True) + print(f'\tLambda bg lower bound: {lambda_bg_lower_bound:.3f}', flush=True) + print(f'\tSparse data: {sparse_data}', flush=True) + print('#1: Done', flush=True) + print('-------------------------------------', flush=True) + + ### 2. Process each chromosome ### + if args.v: + print('#2: Calling peaks from signal\n', flush=True) + + cpus = args.cpus + ctx = mp.get_context('spawn') + + if is_memmap_dir: + # For memmap_np: pass file paths, workers will memory-map the files + chrom_data = [] + for chrom in chrom_ls: + chrom_info = metadata['chromosomes'][chrom] + npy_path = os.path.join(signal_file_path, chrom_info['file']) + first_cut = chrom_info['first_cut'] + chrom_data.append(( + chrom, + npy_path, # Pass path instead of array + first_cut, + threshold, + peak_region_threshold, + min_distance, + min_prominence, + window_size, + lambda_bg_lower_bound, + sparse_data, + args.v + )) + + with ctx.Pool(processes=cpus) as pool: + results = pool.starmap(process_chrom_from_memmap, chrom_data) + else: + # For HDF5: pre-load all data into memory (h5py files can't be pickled) + chrom_data = [] + with h5py.File(signal_file_path, 'r') as sig_file: + for chrom in chrom_ls: + signal_array = sig_file[chrom][:].astype(np.float32) + first_cut = int(sig_file.attrs[chrom]) + chrom_data.append(( + chrom, + signal_array, + first_cut, + threshold, + peak_region_threshold, + min_distance, + min_prominence, + window_size, + lambda_bg_lower_bound, + sparse_data, + args.v + )) + + with ctx.Pool(processes=cpus) as pool: + results = pool.starmap(process_chrom_from_signal, chrom_data) + + # Filter out empty results + results = [r for r in results if r is not None and not r.empty] + + if not results: + print("Warning: No peaks found in any chromosome.", flush=True) + # Write empty output files + DataFrame().to_csv(f'{args.o}/{args.name}_summits.narrowPeak', + sep='\t', header=None, index=None) + DataFrame().to_csv(f'{args.o}/{args.name}_peaks.narrowPeak', + sep='\t', header=None, index=None) + return + + result_df = concat(results) + + if args.v: + print(f'\n\tTotal peaks before filtering: {result_df.shape[0]}', flush=True) + + ### 3. Calculate p-value and q-value ### + if not args.skip_stats: + if args.v: + print('#2: Computing p-value and q-value', flush=True) + + result_df = fseq2.interpolate_poisson_p_value(result_df) + result_df = fseq2.calculate_q_value( + result_df=result_df, + p_thr=args.p_thr, + q_thr=args.q_thr, + num_peaks=args.num_peaks + ) + else: + # Provide dummy statistical values when skipping stats + result_df['-log10_p_value_interpolated'] = result_df['score'] + result_df['q_value'] = result_df['score'] + + if args.v: + print(f'\tPeaks after filtering: {result_df.shape[0]}', flush=True) + print('#2: Done', flush=True) + print('-------------------------------------', flush=True) + + ### 4. Write output ### + if args.v: + print('#3: Writing output', flush=True) + + fseq2.narrowPeak_writer( + result_df=result_df, + peak_type='summit', + name=args.name, + out_dir=args.o, + prior_pad_summit=args.prior_pad_summit, + sort_by=args.sort_by + ) + fseq2.narrowPeak_writer( + result_df=result_df, + peak_type='peak', + name=args.name, + out_dir=args.o, + sort_by=args.sort_by, + standard_narrowpeak=args.standard_narrowpeak + ) + + if args.v: + print(f'\tOutput: {args.o}/{args.name}_summits.narrowPeak', flush=True) + print(f'\tOutput: {args.o}/{args.name}_peaks.narrowPeak', flush=True) + print(f'#3: Done - {result_df.shape[0]} peaks written', flush=True) + print('-------------------------------------', flush=True) + print(f'Thanks for using F-seq{fseq2.__version__}!\n', flush=True) + + return + + +def process_chrom_from_signal(chrom, signal_array, first_cut, threshold, + peak_region_threshold, min_distance, min_prominence, + window_size, lambda_bg_lower_bound, sparse_data, + verbose=False): + """Process a single chromosome from signal array. + + Args: + chrom: Chromosome name + signal_array: numpy array of signal values (float32) + first_cut: Starting genomic position + threshold: Minimum peak height threshold + peak_region_threshold: Threshold for contiguous peak regions + min_distance: Minimum distance between peaks + min_prominence: Minimum prominence for peaks (from original callpeak) + window_size: Window size for lambda calculation + lambda_bg_lower_bound: Lower bound for background lambda + sparse_data: Whether to use sparse data mode for local lambda calculation + verbose: Whether to print progress + + Returns: + pd.DataFrame with peak data or empty DataFrame if no peaks found + """ + start_time = time.time() + + if signal_array.size == 0: + if verbose: + print(f'\t{chrom}: No signal data - skipping', flush=True) + return DataFrame() + + last_cut = first_cut + signal_array.size + + # Call peaks using existing function with all parameters from original callpeak + result_df = fseq2.call_peaks( + chrom=chrom, + first_cut=first_cut, + kdepy_result=signal_array, + min_height=threshold, + peak_region_threshold=peak_region_threshold, + min_distance=min_distance, + min_prominence=min_prominence + ) + + if result_df.empty: + if verbose: + print(f'\t{chrom}: No peaks found', flush=True) + return result_df + + # Calculate query_value and lambda_local for statistical testing + summit_abs_pos_array = result_df['summit'].values - first_cut + + # Calculate query value (average signal around summit) + query_value = fseq2.calculate_query_value( + result_df=result_df, + kdepy_result=signal_array, + summit_abs_pos_array=summit_abs_pos_array, + window_size=window_size, + use_max=False + ) + + # Calculate lambda_bg (background estimate) using the same lower bound as original + lambda_bg = fseq2.calculate_lambda_bg( + kdepy_result_control=signal_array, + window_size=window_size, + lambda_bg_lower_bound=lambda_bg_lower_bound + ) + + # Calculate lambda_local (local background for each peak) with same sparse_data setting + lambda_local = fseq2.find_local_lambda( + control_np_tmp=signal_array, + control_np_tmp_name=False, # Signal is in memory, not file + summit_abs_pos_array=summit_abs_pos_array, + lambda_bg=lambda_bg, + window_size=window_size, + sparse_data=sparse_data, + use_max=False + ) + + result_df['query_value'] = query_value + result_df['lambda_local'] = lambda_local + + end_time = time.time() + if verbose: + print(f'\t{chrom}: first={first_cut}, last={last_cut}, ' + f'peaks={result_df.shape[0]}, completed in {end_time - start_time:.3f} seconds.', + flush=True) + + return result_df + + +def process_chrom_from_memmap(chrom, npy_path, first_cut, threshold, + peak_region_threshold, min_distance, min_prominence, + window_size, lambda_bg_lower_bound, sparse_data, + verbose=False): + """Process a single chromosome from memory-mapped .npy file. + + This function loads the signal array using memory-mapping, which allows + multiple processes to share the same memory without copying. + + Args: + chrom: Chromosome name + npy_path: Path to the .npy file containing signal data + first_cut: Starting genomic position + threshold: Minimum peak height threshold + peak_region_threshold: Threshold for contiguous peak regions + min_distance: Minimum distance between peaks + min_prominence: Minimum prominence for peaks (from original callpeak) + window_size: Window size for lambda calculation + lambda_bg_lower_bound: Lower bound for background lambda + sparse_data: Whether to use sparse data mode for local lambda calculation + verbose: Whether to print progress + + Returns: + pd.DataFrame with peak data or empty DataFrame if no peaks found + """ + # Load signal array with memory mapping (read-only, shared across processes) + signal_array = np.load(npy_path, mmap_mode='r') + + # Delegate to the main processing function + return process_chrom_from_signal( + chrom=chrom, + signal_array=signal_array, + first_cut=first_cut, + threshold=threshold, + peak_region_threshold=peak_region_threshold, + min_distance=min_distance, + min_prominence=min_prominence, + window_size=window_size, + lambda_bg_lower_bound=lambda_bg_lower_bound, + sparse_data=sparse_data, + verbose=verbose + ) diff --git a/fseq2/fseq2.py b/fseq2/fseq2.py index c357e05..2d296ab 100644 --- a/fseq2/fseq2.py +++ b/fseq2/fseq2.py @@ -9,6 +9,7 @@ __version__ = '2.0.4' import functools +import json import math import os import sys @@ -250,20 +251,18 @@ def run_kde(cuts_array, start_array, end_array, strand_array, with h5py.File(params.treatment_np_tmp_name, mode='a', libver='latest') as sig_file: try: sig_file.create_dataset(chrom, - data=np.round(np.divide(kdepy_result, + data=np.divide(kdepy_result, kdepy_result_control, out=np.zeros_like(kdepy_result), - where=kdepy_result_control>lambda_bg_lower_bound_ind), - params.sig_float_precision).astype(np.float16), - dtype='f2')#, compression='gzip') + where=kdepy_result_control>lambda_bg_lower_bound_ind).astype(np.float32), + dtype='f4')#, compression='gzip') sig_file.attrs[chrom] = first_cut except RuntimeError: del sig_file[chrom] sig_file.create_dataset(chrom, - data=np.round(np.divide(kdepy_result, + data=np.divide(kdepy_result, kdepy_result_control, out=np.zeros_like(kdepy_result), - where=kdepy_result_control>lambda_bg_lower_bound_ind), - params.sig_float_precision).astype(np.float16), - dtype='f2')#, compression='gzip') + where=kdepy_result_control>lambda_bg_lower_bound_ind).astype(np.float32), + dtype='f4')#, compression='gzip') sig_file.attrs[chrom] = first_cut # 2. calculate lambda_bg which needs all control signal in memory @@ -377,14 +376,14 @@ def run_kde_wo_control(cuts_array, start_array, end_array, strand_array, chrom, with h5py.File(params.treatment_np_tmp_name, mode='a', libver='latest') as sig_file: try: sig_file.create_dataset(chrom, - data=np.round(kdepy_result, params.sig_float_precision).astype(np.float16), - dtype='f2')#, compression='gzip') + data=kdepy_result, + dtype='f4')#, compression='gzip') sig_file.attrs[chrom] = first_cut except RuntimeError: del sig_file[chrom] sig_file.create_dataset(chrom, - data=np.round(kdepy_result, params.sig_float_precision).astype(np.float16), - dtype='f2')#, compression='gzip') + data=kdepy_result, + dtype='f4')#, compression='gzip') sig_file.attrs[chrom] = first_cut @@ -987,7 +986,8 @@ def gaussian_kernel_fft(sig, sigma): return fftconvolve(sig, gaussian_kernel1d(sigma, 0, lw), mode='same') -def output_sig(sig_format, treatment_np_tmp_name, out_dir, out_name, chr_size_dic, sig_float_precision, gaussian_smooth_sigma=False): +def output_sig(sig_format, treatment_np_tmp_name, out_dir, out_name, chr_size_dic, sig_float_precision, + gaussian_smooth_sigma=False, memmap_out_dir=None): #start_time = time.time() @@ -1021,11 +1021,56 @@ def output_sig(sig_format, treatment_np_tmp_name, out_dir, out_name, chr_size_di step=1) else: for chrom_line in chrom_line_ls: - #kdepy_result = sig_file[chrom_line][:].astype(np.float64) - kdepy_result = np.round(sig_file[chrom_line][:], sig_float_precision).astype(np.float16) - #kdepy_result[kdepy_result == 0] = np.NaN + kdepy_result = sig_file[chrom_line][:].astype(np.float64) + #kdepy_result = np.round(sig_file[chrom_line][:], sig_float_precision).astype(np.float16) + kdepy_result[kdepy_result == 0] = np.NaN output_bw.addEntries(chrom_line, int(sig_file.attrs[chrom_line]), values=kdepy_result, span=1, step=1) + elif sig_format == 'memmap_np': + # Create output directory for memory-mapped numpy files + os.makedirs(memmap_out_dir, exist_ok=True) + + # Build metadata dictionary with chromosome info and parameters + metadata = { + 'chromosomes': {}, + 'parameters': {} + } + + # Copy parameters from HDF5 attributes + param_keys = ['threshold', 'peak_region_threshold', 'min_distance', 'min_prominence', + 'window_size', 'lambda_bg_lower_bound', 'bandwidth', 'fragment_offset', + 'scaling_factor', 'fragment_size', 'feature_length', 'ncuts', 'sparse_data', + 'has_control', 'bandwidth_control', 'fragment_offset_control', 'ncuts_control'] + for key in param_keys: + if key in sig_file.attrs: + val = sig_file.attrs[key] + # Convert numpy types to Python types for JSON serialization + if isinstance(val, (np.integer, np.floating)): + val = val.item() + elif isinstance(val, np.bool_): + val = bool(val) + metadata['parameters'][key] = val + + # Write individual .npy files for each chromosome + for chrom_line in chrom_line_ls: + signal_array = sig_file[chrom_line][:].astype(np.float32) + first_cut = int(sig_file.attrs[chrom_line]) + + # Save as .npy file + npy_path = f'{memmap_out_dir}/{chrom_line}.npy' + np.save(npy_path, signal_array) + + # Store chromosome metadata + metadata['chromosomes'][chrom_line] = { + 'first_cut': first_cut, + 'length': len(signal_array), + 'file': f'{chrom_line}.npy' + } + + # Write metadata JSON + metadata_path = f'{memmap_out_dir}/metadata.json' + with open(metadata_path, 'w') as f: + json.dump(metadata, f, indent=2) #end_time = time.time() #print(f'time = {end_time-start_time:.3f} seconds.') diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..689dda2 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,3 @@ +[build-system] +requires = ["setuptools", "wheel", "cython", "numpy < 2.0"] +build-backend = "setuptools.build_meta" diff --git a/setup.py b/setup.py index 915beb7..8256c61 100644 --- a/setup.py +++ b/setup.py @@ -5,7 +5,6 @@ from setuptools import setup, find_packages, Extension import sys from pathlib import Path -import numpy # with open('README.rst') as readme_file: # readme = readme_file.read() @@ -30,12 +29,14 @@ def readme(): try: from Cython.Build import cythonize + import numpy extensions = cythonize([ Extension("fseq2.idr_2_0_3.inv_cdf", ["fseq2/idr_2_0_3/inv_cdf.pyx", ], include_dirs=[numpy.get_include()]), ]) except ImportError: + import numpy extensions = [ Extension("fseq2.idr_2_0_3.inv_cdf", ["fseq2/idr_2_0_3/inv_cdf.c", ],