PHi-C2

PHi-C2 allows for a physical interpretation of a Hi-C contact matrix. The phic package includes a suite of command line tools.

Installation (with Conda environment)

You can install phic in a clean environment as follows:

conda create -n phic python=3.12
conda activate phic
pip install phic

Without preparing a Python environment, PHi-C2 (=<2.0.13) runs on Google Colab.

Requirements

• PHi-C2 is based on python3.
• Python packages numpy, matplotlib, scipy, click, pandas, MDAnalysis, tqdm, psutil, hictkpy.

To visualize the simulated polymer dynamics and conformations, VMD can be used. Alternatively, the output PSF and DCD files can be viewed directly in a web browser using Mol* without any local installation.

Citation

If you use PHi-C2, please cite:

Soya Shinkai, Hiroya Itoga, Koji Kyoda, and Shuichi Onami. (2022). PHi-C2: interpreting Hi-C data as the dynamic 3D genome state. Bioinformatics 38(21) 4984–4986.

Quick Start

After the installation of phic and downloading of the directory demo, move to the directory demo:

demo/
  run.sh

Then, execute the following script:

./run.sh

This process may take a few minutes.

The demo uses Hi-C data of fission yeast Schizosaccharomyces pombe (chrI, 25-kb resolution, genome-wide KR normalization). The original experiment was reported by Kakui et al., and the included .hic file was generated from those data by Shinkai et al. in their Research Square.

Usage

phic needs a subcommand on the command line interface:

phic SUBCOMMAND [OPTIONS]

Subcommands:
fetch-fileinfo
      ↓
preprocessing
      ↓
optimization
      ├──> plot-optimization
      ├──> dynamics
      ├──> sampling
      ├──> msd
      │    └──> plot-msd
      └──> losstangent
           └──> plot-losstangent

As of version 2.2.2, most subcommands accept experimental --json/--no-json (on by default), --json-path, and --run-uuid options that append a structured analysis log to phic-metadata.json in the analysis directory (NAME/phic-metadata.json by default). These options are provided for internal pipeline integration; the schema and detailed usage will be documented in a future release.

0. fetch-fileinfo

phic fetch-fileinfo [OPTIONS]

Options:
  --input               TEXT     Input Hi-C file (.hic or .mcool format)  [required]

The fetch-fileinfo subcommand is used to inspect the basic metadata of a Hi-C data file. As of version 2.1.1, phic supports both .hic and .mcool formats as input.

Use this command to check available chromosomes, resolution levels, and indexing details in the input file before proceeding with further analysis. This ensures that downstream subcommands reference the correct chromosome names and binning resolutions.

This is a recommended first step when working with new input files.

Example:

phic fetch-fileinfo --input FILENAME.hic

1. preprocessing

phic preprocessing [OPTIONS]

Options:
  --input               TEXT     Input Hi-C file (.hic or .mcool format)  [required]
  --res                 INTEGER  Resolution of the bin size  [required]
  --plt-max-c           FLOAT    Maximum value of contact map  [required]
  --for-high-resolution FLAG     Normalization of contact map for high-resolution case (ex. 1-kb, 500-bp, 200-bp)  [default=False]
  --chr                 TEXT     Target chromosome  [required]
  --grs                 INTEGER  Start position of the target genomic region
  --gre                 INTEGER  End position of the target genomic region
  --norm                TEXT     Type of normalization to apply
  --tolerance           FLOAT    Threshold used to remove segments containing NaN values  [required]
  --name                TEXT     Output directory name (canonical path); if omitted, auto-named by script
  --help                         Show this message and exit.

In version 2.1.1 and later, the input data format has been changed to .hic or .mcool. Additionally, it is now possible to exclude rows and columns containing NaN values from the analysis by specifying their allowed proportion (ranging from 0 to 1) using the tolerance parameter.

As of version 2.2.4, bins with no read data—which hictkpy returns as 0.0 when the contact matrix is densified—are restored to NaN in the input contact matrix. They are therefore counted by --tolerance and excluded from the downstream normalization, instead of being treated as observed zero contacts. Note that sparse (e.g. high-resolution) data now yield much larger NaN proportions per row, so --tolerance may need to be set higher than in previous versions to retain the same segments.

When using the preprocessing subcommand, a directory will be automatically created based on the input Hi-C file name, chromosome number, genomic region of interest (optional), resolution, and normalization method. All subsequent analysis results will be stored in this directory. In the following explanations, we refer to this directory as NAME.

The outputs are as follows:

NAME/
├── C_normalized.npz
├── C_normalized.svg
├── P_normalized.npz
├── P_normalized.svg
├── phic-metadata.json
└── _meta_data/

phic-metadata.json is created here and is subsequently appended to by optimization, plot-optimization, msd, and losstangent as those subcommands are run on NAME (see the note on --json/--no-json above).

The two .npz files can be loaded with numpy.load using the keys below. Here, N is the phic-analysis matrix size (after removing NaN-band rows/columns per --tolerance), and N_input is the input matrix size (the target genomic region binned at the --res resolution, before removing NaN-band rows/columns):

Example:

phic preprocessing --input FILENAME.hic --res 25000 --plt-max-c 0.05 --chr 2 --grs 40000000 --gre 65000000 --norm KR --tolerance 0.4
phic preprocessing --input FILENAME.hic --res 100000 --plt-max-c 0.05 --chr 2 --norm KR --tolerance 0.8

2. optimization

phic optimization [OPTIONS]

Options:
  --name                      TEXT   Target directory name  [required]
  --init-k-backbone           FLOAT  Initial parameter of K_i,i+1  [default=0.5]
  --stop-condition-parameter  FLOAT  Parameter for the stop condition  [default=1e-7]
  --backtracking-factor       FLOAT  Backtracking factor  [default=0.7]
  --gradient-degree           INT    Gradient used for optimizing of K  [default=2]
  --help                             Show this message and exit.

The outputs are the followings:

NAME/
└── data_optimization/
    ├── K_optimized.npz
    ├── C_optimized.npz
    ├── P_optimized.npz
    └── optimization.log

As of version 2.2.1, optimization also produces C_optimized.npz and P_optimized.npz (previously generated by plot-optimization), so the downstream analyses no longer require running plot-optimization first. As of version 2.2.2, K_optimized.npz and C_optimized.npz also include a _with_nan_bands variant of each key, expanded to the full input genomic coordinate system (matching P_normalized.npz); plot-optimization now reads these directly instead of re-inserting NaN rows/columns itself. The three .npz files can be loaded with numpy.load using the keys below:

optimization.log is a tab-separated text file that records the per-iteration trajectory of the optimization. It has a header row followed by one row per accepted step:

It is consumed by plot-optimization to render Cost.svg and Eta.svg, and can also be loaded directly with numpy.loadtxt(path, delimiter="\t", skiprows=1) for custom analysis.

Example:

phic optimization --name NAME

3-1. plot-optimization

phic plot-optimization [OPTIONS]

Options:
  --name        TEXT      Target directory name  [required]
  --plt-max-c   FLOAT     Maximum value of contact map  [required]
  --plt-max-k   FLOAT     Maximum and minimum values of optimized K map  [required]
  --no-figures            Skip all matplotlib figure generation (faster; use when figures are not needed)
  --help                  Show this message and exit.

As of version 2.2.1, plot-optimization reads the arrays pre-computed by optimization and only renders figures; it no longer outputs any .npz files. The --res option has been removed (the resolution is taken from P_normalized.npz).

The outputs are the followings:

NAME/
└── data_optimization/
    ├── C.svg
    ├── Correlation.png
    ├── Correlation_distance_corrected.png
    ├── Cost.svg
    ├── Eta.svg
    ├── K.svg
    └── P.svg

Example:

phic plot-optimization --name NAME --plt-max-c 0.05 --plt-max-k 0.01

3-2. dynamics

phic dynamics [OPTIONS]

Options:
  --name        TEXT      Target directory name  [required]
  --eps         FLOAT     Time step between output frames (exact OU dynamics)
                          [default=1e-3]
  --frame       INTEGER   The number of output frames  [required]
  --sample      INTEGER   The number of output dynamics  [default=1]
  --com-motion            Allow free diffusion of the center-of-mass mode
                          [default=off]
  --seed        INTEGER   Seed of the random numbers  [default=12345678]
  --help                  Show this message and exit.

dynamics integrates the normal modes using Gillespie's (1996) exact one-step solution for the Ornstein--Uhlenbeck process, replacing the previous Heun (predictor--corrector) discretization. The exact scheme has no stability upper bound on the time step and converts normal modes to real-space coordinates only when each output frame is written. The time between output frames is therefore eps; --interval is no longer used. To reproduce the frame spacing from an earlier release, set --eps to the previous eps * interval.

By default, the center-of-mass mode is fixed at the origin. Use --com-motion to allow it to diffuse freely.

The outputs are the followings:

NAME/
└── data_dynamics/
    ├── polymer_N{NUMBER-OF-BEADS}.psf
    ├── sample{SAMPLE-NUMBER}.dcd
    └── sample{SAMPLE-NUMBER}.xyz

Example:

phic dynamics --name NAME --eps 1e-2 --frame 100

3-3. sampling

phic sampling [OPTIONS]

Options:
  --name    TEXT      Target directory name  [required]
  --sample  INTEGER   The number of output conformations  [required]
  --seed    INTEGER   Seed of the random numbers  [default=12345678]
  --help              Show this message and exit.

The outputs are the followings:

NAME/
└── data_sampling/
    ├── polymer_N{NUMBER-OF-BEADS}.psf
    ├── conformations.dcd
    └── conformations.xyz

Example:

phic sampling --name NAME --sample 100

3-4-1. msd

phic msd [OPTIONS]

Options:
  --name  TEXT     Target directory name  [required]
  --help           Show this message and exit.

As of version 2.2.1, the exponent range of the normalized time is automatically determined from the eigenvalues of the Laplacian matrix induced from the optimized polymer network interaction matrix, so the --upper and --lower options have been removed.

The output is the following:

NAME/
└── data_MSD/
    └── MSD_matrix.npz

MSD_matrix.npz contains three arrays and can be loaded with numpy.load using the keys below:

Example:

phic msd --name NAME

3-4-2. plot-msd

phic plot-msd [OPTIONS]

Options:
  --name        TEXT     Target directory name  [required]
  --plt-upper   INTEGER  Upper value of the exponent of the normalized time in the spectrum  [required]
  --plt-lower   INTEGER  Lower value of the exponent of the normalized time in the spectrum  [required]
  --plt-max-log FLOAT    Maximum value of log10 MSD  [required]
  --plt-min-log FLOAT    Minimum value of log10 MSD  [required]
  --aspect      FLOAT    Aspect ratio of the spectrum  [default=0.8]
  --help                 Show this message and exit.

The outputs are the followings:

NAME/
└── data_MSD/
    ├── fig_MSD_curves.png
    └── fig_MSD_spectrum.svg

Example:

phic plot-msd --name NAME --plt-upper 3 --plt-lower 0 --plt-max-log 2.0 --plt-min-log 0.5 --aspect 1.5

3-5-1. losstangent

phic losstangent [OPTIONS]

Options:
  --name    TEXT      Target directory name  [required]
  --help              Show this message and exit.

As of version 2.2.1, the exponent range of the angular frequency is automatically determined from the eigenvalues of the Laplacian matrix induced from the optimized polymer network interaction matrix, so the --upper and --lower options have been removed.

The output is the following:

NAME/
└── data_losstangent/
    └── losstangent_matrix.npz

losstangent_matrix.npz contains three arrays and can be loaded with numpy.load using the keys below:

Example:

phic losstangent --name NAME

3-5-2. plot-losstangent

phic plot-losstangent [OPTIONS]

Options:
  --name          TEXT      Target directory name  [required]
  --plt-upper     INTEGER   Upper value of the exponent of the angular frequency in the spectrum  [required]
  --plt-lower     INTEGER   Lower value of the exponent of the angular frequency in the spectrum  [required]
  --plt-max-log   FLOAT     Maximum value of log10 tanδ  [required]
  --aspect        FLOAT     Aspect ratio of the spectrum  [default=0.8]
  --help                    Show this message and exit.

The output is the following:

NAME/
└── data_losstangent/
    └── fig_losstangent_spectrum.svg

Example:

phic plot-losstangent --name NAME --plt-upper 0 --plt-lower -3 --plt-max-log 0.3 --aspect 1.5