feat: LDSC genetic correlation logic

[xyg123]

✨ Context

LDSC genetic correlation will provide context for assessing genetic relatedness between diseases and traits, the results will be informative of vertical/horizontal pleiotropy.

🛠 What does this PR implement

This PR implements the logic for computing pairwise genetic correlations from raw sumstats, it also includes an overhaul or the heritability estimates to hide the orchestration logic from the gentropy business layer. As well as a general refactoring of the codebase for a modular structure, similar to the original python implementation of LDSC.

🚦 Before submitting

• Do these changes cover one single feature (one change at a time)?
• Did you read the contributor guideline?
• Did you make sure to update the documentation with your changes?
• Did you make sure there is no commented out code in this PR?
• Did you follow conventional commits standards in PR title and commit messages?
• Did you make sure the branch is up-to-date with the dev branch?
• Did you write any new necessary tests?
• Did you make sure the changes pass local tests (make test)?
• Did you make sure the changes pass pre-commit rules (e.g uv run pre-commit run --all-files)?

[project-defiant]

I had to catch up with the Step before I could look over the Genetic Correlation. Overall I think this PR will not be merged, as it has a lot of obsolete, or unnecessary changes. I did not review the individual modules in the ldsc, just the rg.

[project-defiant]

why changed to 20M?

[project-defiant]

import warnings
warnings.warn("import from gentropy.method.ldsc instead. This file will be removed in a future release.", DeprecationWarning,
              stacklevel=2)

# normal module code

See https://stackoverflow.com/questions/30093490/how-to-declare-a-module-deprecated-in-python

[project-defiant]

How do you know that t1 is correctly flipped ?

[project-defiant]

You can reduce the amount of ops

t2_direct = prepared_2.select(
    "chromosome", 
    "position, 
    "ref", 
    "alt", 
    f.col("beta").alias("beta2"), 
    f.col("standardErrro").alias("se")
)
t2_flipped = t2_direct.select(
    "chromosome",
    "position",
    f.col("alt").alias("ref"),
    f.col("ref").alias("alt"),
    (f.col("beta2") * -1).alias("beta2"),
    "se2",
    "n2"
)
t2_combined = t2.unionByName(t2_flipped)

Couple of notes:

• This way you do not need to compute the (beta2 * -1) over n_row = N * 2.
• I would argue that this does not need dedup, as ref != alt and you already deduplicated on the sumstat level in previous functions.
• The flipped column is redundant
• You should definately persist and force the compute before flipping. This will speed up.

[project-defiant]

I think this is not a part of genetic correlation ?

[project-defiant]

Instead of passing 4 x float array (N,1) pass 2 with Zscores. Will reduce the space complexity.

[project-defiant]

Given this is public method, that shall be called by external users, this is the only spot where the shape validation should happen. Every other method downstream, called in this ones call stack do not need to re-verify these assumptions.

Suggested change

	for arr in (beta1, se1, N1, beta2, se2, N2, ld, w_ld):
	arr = np.asarray(arr, dtype=float)
	if arr.ndim != 1:
	raise ValueError("All input arrays must be 1D.")
	beta1 = np.asarray(beta1, dtype=float)
	se1 = np.asarray(se1, dtype=float)
	N1 = np.asarray(N1, dtype=float)
	beta2 = np.asarray(beta2, dtype=float)
	se2 = np.asarray(se2, dtype=float)
	N2 = np.asarray(N2, dtype=float)
	ld = np.asarray(ld, dtype=float)
	w_ld = np.asarray(w_ld, dtype=float)
	z1 = beta1 / se1
	z2 = beta2 / se2
	n = z1.shape[0]
	M_mat = np.array([[float(M_ldsc_scalar)]])
	x = ld.reshape((n, 1))
	w = w_ld.reshape((n, 1))
	# Step 1: Estimate h2 for each trait
	def _run_hsq(z: np.ndarray, N: np.ndarray) -> tuple[float, float]:
	"""Run Hsq regression on a single trait and return (h2, h2_se).
	Args:
	z (np.ndarray): Z-scores for the trait, shape (n_snp,).
	N (np.ndarray): Sample sizes, shape (n_snp,).
	Returns:
	tuple[float, float]: Heritability estimate and its standard error.
	"""
	chisq = z ** 2
	y_h = chisq.reshape((n, 1))
	N_mat = N.reshape((n, 1))
	step_cutoff = twostep if intercept is None else None
	hsqhat = Hsq(
	y_h, x, w, N_mat, M_mat,
	n_blocks=n_blocks,
	intercept=intercept,
	twostep=step_cutoff,
	)
	h2 = float(max(hsqhat.tot, 0.0))
	h2_se = float(hsqhat.tot_se)
	return h2, h2_se
	h2_1, h2_1_se = _run_hsq(z1, N1)
	h2_2, h2_2_se = _run_hsq(z2, N2)
	arrays = (beta1, se1, N1, beta2, se2, N2, ld, w_ld)
	expected_shape = beta1.shape
	if len(expected_shape) != 1 or any(arr.shape != expected_shape for arr in arrays):
	raise ValueError("All input arrays must have the same 1D shape.")
	n = expected_shape[0]
	z1 = (beta1 / se1).reshape((n, 1))
	z2 = (beta2 / se2).reshape((n, 1))
	N1 = N1.reshape((n, 1))
	N2 = N2.reshape((n, 1))
	x = ld.reshape((n, 1))
	w = w_ld.reshape((n, 1))
	M_mat = np.array([[float(M_ldsc_scalar)]])
	def _run_hsq(z: np.ndarray, N: np.ndarray) -> tuple[float, float]:
	chisq = z**2
	step_cutoff = twostep if intercept is None else None
	hsqhat = Hsq(
	chisq,
	x,
	w,
	N,
	M_mat,
	n_blocks=n_blocks,
	intercept=intercept,
	twostep=step_cutoff,
	)
	h2 = float(max(hsqhat.tot, 0.0))
	h2_se = float(hsqhat.tot_se)
	return h2, h2_se
	h2_1, h2_1_se = _run_hsq(z1, N1)
	h2_2, h2_2_se = _run_hsq(z2, N2)

[project-defiant]

I can not find any reference to the delta method, where did you took it from?