"""Module.
Examples
--------
pass
"""
__author__ = "Yeremia Gunawan Adhisantoso"
__license__ = "Clear BSD"
__email__ = "adhisant@tnt.uni-hannover.de"
import typing as t
from pydantic import validate_call
import os
import io
import re
import glob
import warnings
import numpy as np
import pandas as pd
from scipy import stats
from ... import loaders as cm_loaders
from ...consts import DataStructure
from ...structs.conflict_policy import ConflictPolicy
from ..mds.mds_numpy import comp_edm_from_p
from ...converters import convert_to_cm_coo
COO_DELIM = "\t"
PDB_COL_NAMES = [
"ATOM1",
"ATOM2",
"ATOM3",
"ATOM4",
"ATOM_SERIAL_NUMBER",
"X", "Y", "Z",
"A", "B"
]
PDB_USE_COLS = [
# "ATOM_SERIAL_NUMBER",
"X", "Y", "Z",
]
coord_pattern = re.compile(r'[-+]?\d*\.\d+|[-+]?\d+')
POINTS_FFORMAT = ".pdb"
MAPPING_PATTERN = "_coordinate_mapping.txt"
[docs]@validate_call
def gen_h3dg_coo(
chr_region: str,
bin_size_bp: int,
balancing: str,
input_fpath: str,
output_fpath: str,
on_conflict: ConflictPolicy = "error",
) -> None:
"""Generate COO format files for H3DG from a contact matrix file.
Notes
-----
- This function generates both raw and normalized COO files.
- The raw counts file is saved with the suffix `.raw`.
- The normalized counts file is saved with the suffix specified by the `balancing` parameter.
Parameters
----------
chr_region : str
The chromosome region to process.
bin_size_bp : int
The bin size (in base pairs).
balancing : str
The balancing method to use for normalization.
input_fpath : str
The file path to the input contact matrix data.
output_fpath : str
The base file path for the output COO files.
on_conflict : ConflictPolicy, optional
What to do when the output file already exists. One of
``"error"`` (raise), ``"overwrite"`` (replace),
``"skip"`` (return without writing). Defaults to ``"error"``.
Returns
-------
None
Examples
--------
pass
"""
#? Generate raw counts
convert_to_cm_coo(
input_fpath,
f"{output_fpath}.raw",
chr_region,
bin_size_bp,
"NONE",
output_delimiter=COO_DELIM,
on_conflict=on_conflict,
)
#? Generate normalized counts
convert_to_cm_coo(
input_fpath,
f"{output_fpath}.{balancing}",
chr_region,
bin_size_bp,
balancing,
output_delimiter=COO_DELIM,
on_conflict=on_conflict,
)
@validate_call
def find_h3dg_points_fpath(
path: str
) -> str:
"""Find the most recent points file in the specified directory.
Notes
-----
This function searches for files with the `.pdb` extension and returns the most recent one.
Parameters
----------
path : str
The directory path to search for points files.
Returns
-------
str
The file path to the most recent points file.
Examples
--------
pass
"""
points_fpath = sorted(glob.glob(f"{path}/*{POINTS_FFORMAT}"))[-1]
return points_fpath
@validate_call
def parse_h3dg_points(
points_fpath: str,
parser: str = 'regex',
) -> np.ndarray:
"""Parse 3D points from a PDB file using a specified parser.
Notes
-----
- The function supports two parsers: 'regex' and 'pandas'.
- The 'regex' parser uses regular expressions to extract coordinates.
- The 'pandas' parser uses pandas to read the file and extract coordinates.
Parameters
----------
points_fpath : str
The file path to the PDB file containing the points.
parser : str, optional
The parser to use ('regex' or 'pandas'), default is 'regex'.
Returns
-------
np.ndarray
A 2D numpy array of shape (N, 3) containing the parsed points.
Examples
--------
pass
"""
if parser == 'regex':
atoms = []
with open(points_fpath) as f:
for line in f:
if line.startswith('ATOM'):
coords = coord_pattern.findall(line)
atoms.append([float(x) for x in coords[2:2+3]])
num_points = len(atoms)
points = np.empty((num_points, 3))
for i, atom in enumerate(atoms):
points[i, :] = atom
elif parser == 'pandas':
points_lines = []
with open(points_fpath) as f:
for line in f:
if line.startswith('ATOM'):
points_lines.append(line)
f = io.StringIO("".join(points_lines))
points = (
pd.read_csv(
f,
# sep=r"\s+", #? Delimiter is all possible
delim_whitespace=True, #? Delimiter is all possible whitespace
header=None,
names=PDB_COL_NAMES,
usecols=PDB_USE_COLS,
index_col=None,
)
.to_numpy()
)
return points
@validate_call
def find_h3dg_mapping_fpath(
path: str
) -> str:
"""Find the coordinate mapping file in the specified directory.
Notes
-----
This function searches for files with the `_coordinate_mapping.txt` extension and returns the first one found.
Parameters
----------
path : str
The directory path to search for mapping files.
Returns
-------
str
The file path to the coordinate mapping file.
Examples
--------
pass
"""
mapping_fpath = glob.glob(f"{path}/*{MAPPING_PATTERN}")[-1]
return mapping_fpath
@validate_call
def parse_h3dg_mapping(
mapping_fpath: str,
bin_size_bp: int | None = None,
) -> np.ndarray:
"""Parse the coordinate mapping file and returns a mapping matrix.
Notes
-----
- If `bin_size_bp` is not provided, it is inferred from the differences in the loci.
- The function normalizes the loci by the bin size.
Parameters
----------
mapping_fpath : str
The file path to the coordinate mapping file.
bin_size_bp : Optional[int], optional
The bin size (in base pairs) to use for normalization, default is None.
Returns
-------
np.ndarray
A 2D numpy array of shape (N, 2) containing the parsed mapping.
Examples
--------
pass
"""
#? This is the mapping of points to loci (result of simulation)
mapping_mat = pd.read_csv(
mapping_fpath,
sep="\t",
header=None,
names=['loci', 'id'],
index_col=None,
).to_numpy()
if bin_size_bp is None:
bin_size_bp = np.diff(mapping_mat[:, 0]).min()
mapping_mat[:, 0] //= bin_size_bp
return mapping_mat
@validate_call
def load_h3dg_points(
res_path: str,
points_fpath: str | None = None,
parser: str = 'regex',
mapping_fpath: str | None = None,
bin_size_bp: int | None = None,
num_bins: int | None = None,
rc_ids: list | None = None,
def_coor: float = np.nan,
) -> np.ndarray:
"""Load and processes 3D points from a PDB file and a coordinate mapping file.
Notes
-----
- If `points_fpath` is not provided, the function finds the most recent points file in `res_path`.
- If `mapping_fpath` is not provided, the function finds the coordinate mapping file in `res_path`.
- The function handles missing data and invalid loci by setting them to `def_coor`.
Parameters
----------
res_path : str
The directory path containing the points and mapping files.
points_fpath : Optional[str], optional
The file path to the PDB file containing the points, default is None.
parser : str, optional
The parser to use ('regex' or 'pandas'), default is 'regex'.
mapping_fpath : Optional[str], optional
The file path to the coordinate mapping file, default is None.
bin_size_bp : Optional[int], optional
The bin size (in base pairs) to use for normalization, default is None.
num_bins : Optional[int], optional
The number of bins, default is None.
rc_ids : Optional[List]], optional
Row and column IDs to filter valid loci, default is None.
def_coor : float, optional
The default coordinate value for missing or invalid loci, default is np.nan.
Returns
-------
np.ndarray
A 2D numpy array of shape (N, 3) containing the processed points.
Examples
--------
pass
"""
if points_fpath is None:
points_fpath = find_h3dg_points_fpath(res_path)
points = parse_h3dg_points(
points_fpath,
parser=parser,
)
if mapping_fpath is None:
mapping_fpath = find_h3dg_mapping_fpath(res_path)
mapping_mat = parse_h3dg_mapping(
mapping_fpath,
bin_size_bp=bin_size_bp
)
if num_bins is None:
num_bins = mapping_mat[:, 0].max() + 1
new_points = np.full((num_bins, 3), def_coor)
new_points[mapping_mat[:, 0], :] = points[mapping_mat[:, 1], :]
points = new_points
if rc_ids is not None:
#? Create mapping from row/col ids to points
#? Reason: not all loci are valid, yet the points contains only valid points
row_ids, col_ids = rc_ids
unique_data_ids = np.unique([row_ids, col_ids])
mask = np.zeros(num_bins, dtype=bool)
mask[unique_data_ids] = True
points[~mask, :] = def_coor
return points
@validate_call
def comp_h3dg_obj_perf(
region1: str,
bin_size_bp: int,
balancing: str,
cm_fpath: str,
points_fpath: str,
mappings_fpath: str,
region2: str | None = None,
) -> dict:
"""Compute the performance metrics (Spearman and Pearson correlation) for Euclidean distances predicted by Hierarchical 3D Genome.
Notes
-----
This function computes the Spearman correlation between contact counts and Euclidean distances derived from 3D points.
It handles edge cases such as invalid loci and missing data.
Parameters
----------
region1 : str
The chromosome region for the first chromosome.
bin_size_bp : int
The bin size (in base pairs).
balancing : str
The balancing method to use.
cm_fpath : str
The file path to the input contact matrix data.
points_fpath : str
The file path to the points data (pdb or xyz).
mappings_fpath : str
The file path to the mappings data.
region2 : Optional[str], optional
The chromosome region for the second chromosome (default is None).
Returns
-------
dict
A dictionary containing the region, correlation, and data ratio.
Examples
--------
pass
"""
if not os.path.exists(points_fpath):
raise FileNotFoundError(f"Points file {points_fpath} not found.")
count_df = cm_loaders.load_cm_data(
cm_fpath,
bin_size_bp,
region1,
balancing=balancing,
region2=region2,
output_format=DataStructure.DF
)
row_ids = count_df[cm_loaders.ROW_IDS_COLNAME].to_numpy()
col_ids = count_df[cm_loaders.COL_IDS_COLNAME].to_numpy()
counts = count_df[cm_loaders.COUNTS_COLNAME].to_numpy()
#? Create mapping from row/col ids to points
#? Reason: not all loci are valid, yet the points contains only valid points
unique_data_ids = np.unique([row_ids, col_ids])
points_lines = []
with open(points_fpath) as f:
for line in f:
if line.startswith('ATOM'):
points_lines.append(line)
f = io.StringIO("".join(points_lines))
points = (
pd.read_csv(
f,
# sep=r"\s+", #? Delimiter is all possible
delim_whitespace=True, #? Delimiter is all possible whitespace
header=None,
names=PDB_COL_NAMES,
usecols=PDB_USE_COLS,
index_col=None,
)
.to_numpy()
)
#? This is the mapping of points to loci (result of simulation)
mapping_df = pd.read_csv(
mappings_fpath,
sep="\t",
header=None,
names=['loci', 'id'],
index_col=None,
).to_numpy()
mapping_df[:, 0] //= bin_size_bp
#? Filter out data that got removed during simulation
unique_data_ids_mask = np.isin(unique_data_ids, mapping_df[:, 0])
removed_data_ids = unique_data_ids[~unique_data_ids_mask]
valid_data_mask = ~np.isin(row_ids, removed_data_ids)
valid_data_mask &= ~np.isin(col_ids, removed_data_ids)
row_ids = row_ids[valid_data_mask]
col_ids = col_ids[valid_data_mask]
counts = counts[valid_data_mask]
unique_data_ids = unique_data_ids[unique_data_ids_mask]
mapping = np.searchsorted(unique_data_ids, np.arange(unique_data_ids.max()+1))
new_row_ids = mapping[row_ids]
new_col_ids = mapping[col_ids]
edm = comp_edm_from_p(
points,
new_row_ids,
new_col_ids
)
res = stats.spearmanr(counts, edm)
spearman_r = res.correlation
res = stats.pearsonr(counts, edm)
pearson_r = res.correlation
data_ratio = valid_data_mask.sum() / valid_data_mask.size
output_dict = {
'region': region1,
'spearman_r': spearman_r,
'pearson_r': pearson_r,
'data_ratio': data_ratio,
}
return output_dict