"""Module.
Examples
--------
pass
"""
__author__ = "Yeremia Gunawan Adhisantoso"
__license__ = "Clear BSD"
__email__ = "adhisant@tnt.uni-hannover.de"
import typing as t
import os
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 ..mds.mds_numpy import comp_edm_from_p
from ...structs.conflict_policy import ConflictPolicy
from ...converters import convert_to_cm_coo
SHNEIGH_COOR_COL_NAMES = ["IDS", "X", "Y", "Z"]
COO_DELIM = "\t"
[docs]def gen_shneigh_coo(
chr_region: str,
bin_size_bp: int,
balancing: str,
input_fpath: str,
output_fpath: str,
on_conflict: ConflictPolicy = "error",
):
"""Generate a COO format file for SHNeigh.
Notes
-----
This function converts the input data to a COO format file suitable for SHNeigh.
Parameters
----------
chr_region : str
The chromosome region.
bin_size_bp : int
The bin size (in base pairs).
balancing : str
The balancing method.
input_fpath : str
The path to the input contact matrix.
output_fpath : str
The path to the output COO file.
on_conflict : ConflictPolicy, optional
What to do when ``output_fpath`` already exists. One of
``"error"`` (raise), ``"overwrite"`` (replace),
``"skip"`` (return without writing). Defaults to ``"error"``.
Returns
-------
None
"""
convert_to_cm_coo(
chr_region,
bin_size_bp,
balancing,
input_fpath,
output_fpath,
on_conflict=on_conflict,
gen_pseudo_weights=False,
output_delimiter=COO_DELIM,
)
def comp_shneigh_obj_perf(
region1: str,
bin_size_bp: int,
balancing: str,
input_fpath: str,
points_fpath: str,
region2: str | None = None,
) -> dict:
"""Compute the performance metrics (Spearman and Pearson correlation) for Euclidean distances predicted by SHNeigh.
Notes
-----
This function computes the Spearman rank correlation between contact counts and Euclidean distances
derived from 3D coordinates. It handles cases where some loci are invalid and ensures that only valid
data points are used in the computation. If the points file does not exist, a `FileNotFoundError` is raised.
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.
input_fpath : str
The path to the input file containing contact counts.
points_fpath : str
The path to the file containing 3D coordinates.
region2 : Optional[str], optional
The chromosome region for the second chromosome, by default None.
Returns
-------
dict
A dictionary containing the region, Spearman rank correlation, Pearson correlation, and data ratio.
"""
if not os.path.exists(points_fpath):
raise FileNotFoundError(f"Points file {points_fpath} not found.")
count_df = cm_loaders.load_cm_data(
input_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_df = pd.read_csv(
points_fpath,
delim_whitespace=True, #? Delimiter is all possible whitespace
header=None,
names=SHNEIGH_COOR_COL_NAMES,
index_col=None,
)
points = points_df.iloc[:, 1:4].to_numpy()
unique_data_ids_mask = np.isin(unique_data_ids, points_df["IDS"])
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