OSCR

Locus- and Gene-Level Insights into the Inverse Association Between Alzheimer's Disease and Cancer.

Code ↔ Paper

1 match between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 1 match
  1. [1] § 4. Methods › 4.2. SNP-Based Heritability Analysis (Observed Scale) ↔ ldscore/regressions.py, lines 538 to 677 · score 0.61 · Score Regression, LD scores, obs, heritability, partition, variance · Report a mismatch

Paper

Loaded from Europe PMC by your browser; OSCR does not store it: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 743 lines · 29 KB · GPL-3.0 · 1 match

  1. '''
  2. (c) 2014 Brendan Bulik-Sullivan and Hilary Finucane
  4. Estimators of heritability and genetic correlation.
  6. Shape convention is (n_snp, n_annot) for all classes.
  7. Last column = intercept.
  9. '''
  10. from __future__ import division
  11. import numpy as np
  12. import pandas as pd
  13. from scipy.stats import norm, chi2
  14. import jackknife as jk
  15. from irwls import IRWLS
  16. from scipy.stats import t as tdist
  17. from collections import namedtuple
  18. np.seterr(divide='raise', invalid='raise')
  20. s = lambda x: remove_brackets(str(np.matrix(x)))
  23. def update_separators(s, ii):
  24. '''s are separators with ii masked. Returns unmasked separators.'''
  25. maplist = np.arange(len(ii))[np.squeeze(ii)]
  26. mask_to_unmask = lambda i: maplist[i]
  27. t = np.apply_along_axis(mask_to_unmask, 0, s[1:-1])
  28. t = np.hstack(((0), t, (len(ii))))
  29. return t
  32. def p_z_norm(est, se):
  33. '''Convert estimate and se to Z-score and P-value.'''
  34. try:
  35. Z = est / se
  36. except (FloatingPointError, ZeroDivisionError):
  37. Z = float('inf')
  39. P = chi2.sf(Z ** 2, 1, loc=0, scale=1) # 0 if Z=inf
  40. return P, Z
  43. def remove_brackets(x):
  44. '''Get rid of brackets and trailing whitespace in numpy arrays.'''
  45. return x.replace('[', '').replace(']', '').strip()
  48. def append_intercept(x):
  49. '''
  50. Appends an intercept term to the design matrix for a linear regression.
  52. Parameters
  53. ----------
  54. x : np.matrix with shape (n_row, n_col)
  55. Design matrix. Columns are predictors; rows are observations.
  57. Returns
  58. -------
  59. x_new : np.matrix with shape (n_row, n_col+1)
  60. Design matrix with intercept term appended.
  62. '''
  63. n_row = x.shape[0]
  64. intercept = np.ones((n_row, 1))
  65. x_new = np.concatenate((x, intercept), axis=1)
  66. return x_new
  69. def remove_intercept(x):
  70. '''Removes the last column.'''
  71. n_col = x.shape[1]
  72. return x[:, 0:n_col - 1]
  75. def gencov_obs_to_liab(gencov_obs, P1, P2, K1, K2):
  76. '''
  77. Converts genetic covariance on the observed scale in an ascertained sample to genetic
  78. covariance on the liability scale in the population
  80. Parameters
  81. ----------
  82. gencov_obs : float
  83. Genetic covariance on the observed scale in an ascertained sample.
  84. P1, P2 : float in (0,1)
  85. Prevalences of phenotypes 1,2 in the sample.
  86. K1, K2 : float in (0,1)
  87. Prevalences of phenotypes 1,2 in the population.
  89. Returns
  90. -------
  91. gencov_liab : float
  92. Genetic covariance between liabilities in the population.
  94. Note: if a trait is a QT, set P = K = None.
  96. '''
  97. c1 = 1
  98. c2 = 1
  99. if P1 is not None and K1 is not None:
  100. c1 = np.sqrt(h2_obs_to_liab(1, P1, K1))
  101. if P2 is not None and K2 is not None:
  102. c2 = np.sqrt(h2_obs_to_liab(1, P2, K2))
  104. return gencov_obs * c1 * c2
  107. def h2_obs_to_liab(h2_obs, P, K):
  108. '''
  109. Converts heritability on the observed scale in an ascertained sample to heritability
  110. on the liability scale in the population.
  112. Parameters
  113. ----------
  114. h2_obs : float
  115. Heritability on the observed scale in an ascertained sample.
  116. P : float in (0,1)
  117. Prevalence of the phenotype in the sample.
  118. K : float in (0,1)
  119. Prevalence of the phenotype in the population.
  121. Returns
  122. -------
  123. h2_liab : float
  124. Heritability of liability in the population.
  126. '''
  127. if np.isnan(P) and np.isnan(K):
  128. return h2_obs
  129. if K <= 0 or K >= 1:
  130. raise ValueError('K must be in the range (0,1)')
  131. if P <= 0 or P >= 1:
  132. raise ValueError('P must be in the range (0,1)')
  134. thresh = norm.isf(K)
  135. conversion_factor = K ** 2 * \
  136. (1 - K) ** 2 / (P * (1 - P) * norm.pdf(thresh) ** 2)
  137. return h2_obs * conversion_factor
  140. class LD_Score_Regression(object):
  142. def __init__(self, y, x, w, N, M, n_blocks, intercept=None, slow=False, step1_ii=None, old_weights=False):
  143. for i in [y, x, w, M, N]:
  144. try:
  145. if len(i.shape) != 2:
  146. raise TypeError('Arguments must be 2D arrays.')
  147. except AttributeError:
  148. raise TypeError('Arguments must be arrays.')
  150. n_snp, self.n_annot = x.shape
  151. if any(i.shape != (n_snp, 1) for i in [y, w, N]):
  152. raise ValueError(
  153. 'N, weights and response (z1z2 or chisq) must have shape (n_snp, 1).')
  154. if M.shape != (1, self.n_annot):
  155. raise ValueError('M must have shape (1, n_annot).')
  157. M_tot = float(np.sum(M))
  158. x_tot = np.sum(x, axis=1).reshape((n_snp, 1))
  159. self.constrain_intercept = intercept is not None
  160. self.intercept = intercept
  161. self.n_blocks = n_blocks
  162. tot_agg = self.aggregate(y, x_tot, N, M_tot, intercept)
  163. initial_w = self._update_weights(
  164. x_tot, w, N, M_tot, tot_agg, intercept)
  165. Nbar = np.mean(N) # keep condition number low
  166. x = np.multiply(N, x) / Nbar
  167. if not self.constrain_intercept:
  168. x, x_tot = append_intercept(x), append_intercept(x_tot)
  169. yp = y
  170. else:
  171. yp = y - intercept
  172. self.intercept_se = 'NA'
  173. del y
  174. self.twostep_filtered = None
  175. if step1_ii is not None and self.constrain_intercept:
  176. raise ValueError(
  177. 'twostep is not compatible with constrain_intercept.')
  178. elif step1_ii is not None and self.n_annot > 1:
  179. raise ValueError(
  180. 'twostep not compatible with partitioned LD Score yet.')
  181. elif step1_ii is not None:
  182. n1 = np.sum(step1_ii)
  183. self.twostep_filtered = n_snp - n1
  184. x1 = x[np.squeeze(step1_ii), :]
  185. yp1, w1, N1, initial_w1 = map(
  186. lambda a: a[step1_ii].reshape((n1, 1)), (yp, w, N, initial_w))
  187. update_func1 = lambda a: self._update_func(
  188. a, x1, w1, N1, M_tot, Nbar, ii=step1_ii)
  189. step1_jknife = IRWLS(
  190. x1, yp1, update_func1, n_blocks, slow=slow, w=initial_w1)
  191. step1_int, _ = self._intercept(step1_jknife)
  192. yp = yp - step1_int
  193. x = remove_intercept(x)
  194. x_tot = remove_intercept(x_tot)
  195. update_func2 = lambda a: self._update_func(
  196. a, x_tot, w, N, M_tot, Nbar, step1_int)
  197. s = update_separators(step1_jknife.separators, step1_ii)
  198. step2_jknife = IRWLS(
  199. x, yp, update_func2, n_blocks, slow=slow, w=initial_w, separators=s)
  200. c = np.sum(np.multiply(initial_w, x)) / \
  201. np.sum(np.multiply(initial_w, np.square(x)))
  202. jknife = self._combine_twostep_jknives(
  203. step1_jknife, step2_jknife, M_tot, c, Nbar)
  204. elif old_weights:
  205. initial_w = np.sqrt(initial_w)
  206. x = IRWLS._weight(x, initial_w)
  207. y = IRWLS._weight(yp, initial_w)
  208. jknife = jk.LstsqJackknifeFast(x, y, n_blocks)
  209. else:
  210. update_func = lambda a: self._update_func(
  211. a, x_tot, w, N, M_tot, Nbar, intercept)
  212. jknife = IRWLS(
  213. x, yp, update_func, n_blocks, slow=slow, w=initial_w)
  215. self.coef, self.coef_cov, self.coef_se = self._coef(jknife, Nbar)
  216. self.cat, self.cat_cov, self.cat_se =\
  217. self._cat(jknife, M, Nbar, self.coef, self.coef_cov)
  219. self.tot, self.tot_cov, self.tot_se = self._tot(self.cat, self.cat_cov)
  220. self.prop, self.prop_cov, self.prop_se =\
  221. self._prop(jknife, M, Nbar, self.cat, self.tot)
  223. self.enrichment, self.M_prop = self._enrichment(
  224. M, M_tot, self.cat, self.tot)
  225. if not self.constrain_intercept:
  226. self.intercept, self.intercept_se = self._intercept(jknife)
  228. self.jknife = jknife
  229. self.tot_delete_values = self._delete_vals_tot(jknife, Nbar, M)
  230. self.part_delete_values = self._delete_vals_part(jknife, Nbar, M)
  231. if not self.constrain_intercept:
  232. self.intercept_delete_values = jknife.delete_values[
  233. :, self.n_annot]
  235. self.M = M
  237. @classmethod
  238. def aggregate(cls, y, x, N, M, intercept=None):
  239. if intercept is None:
  240. intercept = cls.__null_intercept__
  242. num = M * (np.mean(y) - intercept)
  243. denom = np.mean(np.multiply(x, N))
  244. return num / denom
  246. def _update_func(self, x, ref_ld_tot, w_ld, N, M, Nbar, intercept=None, ii=None):
  247. raise NotImplementedError
  249. def _delete_vals_tot(self, jknife, Nbar, M):
  250. '''Get delete values for total h2 or gencov.'''
  251. n_annot = self.n_annot
  252. tot_delete_vals = jknife.delete_values[
  253. :, 0:n_annot] # shape (n_blocks, n_annot)
  254. # shape (n_blocks, 1)
  255. tot_delete_vals = np.dot(tot_delete_vals, M.T) / Nbar
  256. return tot_delete_vals
  258. def _delete_vals_part(self, jknife, Nbar, M):
  259. '''Get delete values for partitioned h2 or gencov.'''
  260. n_annot = self.n_annot
  261. return jknife.delete_values[:, 0:n_annot] / Nbar
  263. def _coef(self, jknife, Nbar):
  264. '''Get coefficient estimates + cov from the jackknife.'''
  265. n_annot = self.n_annot
  266. coef = jknife.est[0, 0:n_annot] / Nbar
  267. coef_cov = jknife.jknife_cov[0:n_annot, 0:n_annot] / Nbar ** 2
  268. coef_se = np.sqrt(np.diag(coef_cov))
  269. return coef, coef_cov, coef_se
  271. def _cat(self, jknife, M, Nbar, coef, coef_cov):
  272. '''Convert coefficients to per-category h2 or gencov.'''
  273. cat = np.multiply(M, coef)
  274. cat_cov = np.multiply(np.dot(M.T, M), coef_cov)
  275. cat_se = np.sqrt(np.diag(cat_cov))
  276. return cat, cat_cov, cat_se
  278. def _tot(self, cat, cat_cov):
  279. '''Convert per-category h2 to total h2 or gencov.'''
  280. tot = np.sum(cat)
  281. tot_cov = np.sum(cat_cov)
  282. tot_se = np.sqrt(tot_cov)
  283. return tot, tot_cov, tot_se
  285. def _prop(self, jknife, M, Nbar, cat, tot):
  286. '''Convert total h2 and per-category h2 to per-category proportion h2 or gencov.'''
  287. n_annot = self.n_annot
  288. n_blocks = jknife.delete_values.shape[0]
  289. numer_delete_vals = np.multiply(
  290. M, jknife.delete_values[:, 0:n_annot]) / Nbar # (n_blocks, n_annot)
  291. denom_delete_vals = np.sum(
  292. numer_delete_vals, axis=1).reshape((n_blocks, 1))
  293. denom_delete_vals = np.dot(denom_delete_vals, np.ones((1, n_annot)))
  294. prop = jk.RatioJackknife(
  295. cat / tot, numer_delete_vals, denom_delete_vals)
  296. return prop.est, prop.jknife_cov, prop.jknife_se
  298. def _enrichment(self, M, M_tot, cat, tot):
  299. '''Compute proportion of SNPs per-category enrichment for h2 or gencov.'''
  300. M_prop = M / M_tot
  301. enrichment = np.divide(cat, M) / (tot / M_tot)
  302. return enrichment, M_prop
  304. def _intercept(self, jknife):
  305. '''Extract intercept and intercept SE from block jackknife.'''
  306. n_annot = self.n_annot
  307. intercept = jknife.est[0, n_annot]
  308. intercept_se = jknife.jknife_se[0, n_annot]
  309. return intercept, intercept_se
  311. def _combine_twostep_jknives(self, step1_jknife, step2_jknife, M_tot, c, Nbar=1):
  312. '''Combine free intercept and constrained intercept jackknives for --two-step.'''
  313. n_blocks, n_annot = step1_jknife.delete_values.shape
  314. n_annot -= 1
  315. if n_annot > 2:
  316. raise ValueError(
  317. 'twostep not yet implemented for partitioned LD Score.')
  319. step1_int, _ = self._intercept(step1_jknife)
  320. est = np.hstack(
  321. (step2_jknife.est, np.array(step1_int).reshape((1, 1))))
  322. delete_values = np.zeros((n_blocks, n_annot + 1))
  323. delete_values[:, n_annot] = step1_jknife.delete_values[:, n_annot]
  324. delete_values[:, 0:n_annot] = step2_jknife.delete_values -\
  325. c * (step1_jknife.delete_values[:, n_annot] -
  326. step1_int).reshape((n_blocks, n_annot)) # check this
  327. pseudovalues = jk.Jackknife.delete_values_to_pseudovalues(
  328. delete_values, est)
  329. jknife_est, jknife_var, jknife_se, jknife_cov = jk.Jackknife.jknife(
  330. pseudovalues)
  331. jknife = namedtuple('jknife',
  332. ['est', 'jknife_se', 'jknife_est', 'jknife_var', 'jknife_cov', 'delete_values'])
  333. return jknife(est, jknife_se, jknife_est, jknife_var, jknife_cov, delete_values)
  336. class Hsq(LD_Score_Regression):
  338. __null_intercept__ = 1
  340. def __init__(self, y, x, w, N, M, n_blocks=200, intercept=None, slow=False, twostep=None, old_weights=False):
  341. step1_ii = None
  342. if twostep is not None:
  343. step1_ii = y < twostep
  345. LD_Score_Regression.__init__(self, y, x, w, N, M, n_blocks, intercept=intercept,
  346. slow=slow, step1_ii=step1_ii, old_weights=old_weights)
  347. self.mean_chisq, self.lambda_gc = self._summarize_chisq(y)
  348. if not self.constrain_intercept:
  349. self.ratio, self.ratio_se = self._ratio(
  350. self.intercept, self.intercept_se, self.mean_chisq)
  352. def _update_func(self, x, ref_ld_tot, w_ld, N, M, Nbar, intercept=None, ii=None):
  353. '''
  354. Update function for IRWLS
  356. x is the output of np.linalg.lstsq.
  357. x[0] is the regression coefficients
  358. x[0].shape is (# of dimensions, 1)
  359. the last element of x[0] is the intercept.
  361. intercept is None --> free intercept
  362. intercept is not None --> constrained intercept
  363. '''
  364. hsq = M * x[0][0] / Nbar
  365. if intercept is None:
  366. intercept = max(x[0][1]) # divide by zero error if intercept < 0
  367. else:
  368. if ref_ld_tot.shape[1] > 1:
  369. raise ValueError(
  370. 'Design matrix has intercept column for constrained intercept regression!')
  372. ld = ref_ld_tot[:, 0].reshape(w_ld.shape) # remove intercept
  373. w = self.weights(ld, w_ld, N, M, hsq, intercept, ii)
  374. return w
  376. def _summarize_chisq(self, chisq):
  377. '''Compute mean chi^2 and lambda_GC.'''
  378. mean_chisq = np.mean(chisq)
  379. # median and matrix don't play nice
  380. lambda_gc = np.median(np.asarray(chisq)) / 0.4549
  381. return mean_chisq, lambda_gc
  383. def _ratio(self, intercept, intercept_se, mean_chisq):
  384. '''Compute ratio (intercept - 1) / (mean chi^2 -1 ).'''
  385. if mean_chisq > 1:
  386. ratio_se = intercept_se / (mean_chisq - 1)
  387. ratio = (intercept - 1) / (mean_chisq - 1)
  388. else:
  389. ratio = 'NA'
  390. ratio_se = 'NA'
  392. return ratio, ratio_se
  394. def _overlap_output(self, category_names, overlap_matrix, M_annot, M_tot, print_coefficients):
  395. '''LD Score regression summary for overlapping categories.'''
  396. overlap_matrix_prop = np.zeros([self.n_annot,self.n_annot])
  397. for i in range(self.n_annot):
  398. overlap_matrix_prop[i, :] = overlap_matrix[i, :] / M_annot
  400. prop_hsq_overlap = np.dot(
  401. overlap_matrix_prop, self.prop.T).reshape((1, self.n_annot))
  402. prop_hsq_overlap_var = np.diag(
  403. np.dot(np.dot(overlap_matrix_prop, self.prop_cov), overlap_matrix_prop.T))
  404. prop_hsq_overlap_se = np.sqrt(
  405. np.maximum(0, prop_hsq_overlap_var)).reshape((1, self.n_annot))
  406. one_d_convert = lambda x: np.array(x).reshape(np.prod(x.shape))
  407. prop_M_overlap = M_annot / M_tot
  408. enrichment = prop_hsq_overlap / prop_M_overlap
  409. enrichment_se = prop_hsq_overlap_se / prop_M_overlap
  410. overlap_matrix_diff = np.zeros([self.n_annot,self.n_annot])
  411. for i in range(self.n_annot):
  412. if not M_tot == M_annot[0,i]:
  413. overlap_matrix_diff[i, :] = overlap_matrix[i,:]/M_annot[0,i] - \
  414. (M_annot - overlap_matrix[i,:]) / (M_tot-M_annot[0,i])
  416. diff_est = np.dot(overlap_matrix_diff,self.coef)
  417. diff_cov = np.dot(np.dot(overlap_matrix_diff,self.coef_cov),overlap_matrix_diff.T)
  418. diff_se = np.sqrt(np.diag(diff_cov))
  419. diff_p = ['NA' if diff_se[i]==0 else 2*tdist.sf(abs(diff_est[i]/diff_se[i]),self.n_blocks) \
  420. for i in range(self.n_annot)]
  422. df = pd.DataFrame({
  423. 'Category': category_names,
  424. 'Prop._SNPs': one_d_convert(prop_M_overlap),
  425. 'Prop._h2': one_d_convert(prop_hsq_overlap),
  426. 'Prop._h2_std_error': one_d_convert(prop_hsq_overlap_se),
  427. 'Enrichment': one_d_convert(enrichment),
  428. 'Enrichment_std_error': one_d_convert(enrichment_se),
  429. 'Enrichment_p':diff_p,
  430. 'Coefficient': one_d_convert(self.coef),
  431. 'Coefficient_std_error': self.coef_se,
  432. 'Coefficient_z-score': one_d_convert(self.coef) / one_d_convert(self.coef_se)
  433. })
  434. if print_coefficients:
  435. df = df[['Category', 'Prop._SNPs', 'Prop._h2', 'Prop._h2_std_error',
  436. 'Enrichment','Enrichment_std_error', 'Enrichment_p',
  437. 'Coefficient', 'Coefficient_std_error','Coefficient_z-score']]
  438. else:
  439. df = df[['Category', 'Prop._SNPs', 'Prop._h2', 'Prop._h2_std_error',
  440. 'Enrichment','Enrichment_std_error', 'Enrichment_p']]
  441. return df
  444. def summary(self, ref_ld_colnames=None, P=None, K=None, overlap=False):
  445. '''Print summary of the LD Score Regression.'''
  446. if P is not None and K is not None:
  447. T = 'Liability'
  448. c = h2_obs_to_liab(1, P, K)
  449. else:
  450. T = 'Observed'
  451. c = 1
  453. out = ['Total ' + T + ' scale h2: ' +
  454. s(c * self.tot) + ' (' + s(c * self.tot_se) + ')']
  455. if self.n_annot > 1:
  456. if ref_ld_colnames is None:
  457. ref_ld_colnames = ['CAT_' + str(i)
  458. for i in xrange(self.n_annot)]
  460. out.append('Categories: ' + ' '.join(ref_ld_colnames))
  462. if not overlap:
  463. out.append(T + ' scale h2: ' + s(c * self.cat))
  464. out.append(T + ' scale h2 SE: ' + s(c * self.cat_se))
  465. out.append('Proportion of SNPs: ' + s(self.M_prop))
  466. out.append('Proportion of h2g: ' + s(self.prop))
  467. out.append('Enrichment: ' + s(self.enrichment))
  468. out.append('Coefficients: ' + s(self.coef))
  469. out.append('Coefficient SE: ' + s(self.coef_se))
  471. out.append('Lambda GC: ' + s(self.lambda_gc))
  472. out.append('Mean Chi^2: ' + s(self.mean_chisq))
  473. if self.constrain_intercept:
  474. out.append(
  475. 'Intercept: constrained to {C}'.format(C=s(self.intercept)))
  476. else:
  477. out.append(
  478. 'Intercept: ' + s(self.intercept) + ' (' + s(self.intercept_se) + ')')
  479. if self.mean_chisq > 1:
  480. if self.ratio < 0:
  481. out.append(
  482. 'Ratio < 0 (usually indicates GC correction).')
  483. else:
  484. out.append(
  485. 'Ratio: ' + s(self.ratio) + ' (' + s(self.ratio_se) + ')')
  486. else:
  487. out.append('Ratio: NA (mean chi^2 < 1)')
  489. return remove_brackets('\n'.join(out))
  491. def _update_weights(self, ld, w_ld, N, M, hsq, intercept, ii=None):
  492. if intercept is None:
  493. intercept = self.__null_intercept__
  495. return self.weights(ld, w_ld, N, M, hsq, intercept, ii)
  497. @classmethod
  498. def weights(cls, ld, w_ld, N, M, hsq, intercept=None, ii=None):
  499. '''
  500. Regression weights.
  502. Parameters
  503. ----------
  504. ld : np.matrix with shape (n_snp, 1)
  505. LD Scores (non-partitioned).
  506. w_ld : np.matrix with shape (n_snp, 1)
  507. LD Scores (non-partitioned) computed with sum r^2 taken over only those SNPs included
  508. in the regression.
  509. N : np.matrix of ints > 0 with shape (n_snp, 1)
  510. Number of individuals sampled for each SNP.
  511. M : float > 0
  512. Number of SNPs used for estimating LD Score (need not equal number of SNPs included in
  513. the regression).
  514. hsq : float in [0,1]
  515. Heritability estimate.
  517. Returns
  518. -------
  519. w : np.matrix with shape (n_snp, 1)
  520. Regression weights. Approx equal to reciprocal of conditional variance function.
  522. '''
  523. M = float(M)
  524. if intercept is None:
  525. intercept = 1
  527. hsq = max(hsq, 0.0)
  528. hsq = min(hsq, 1.0)
  529. ld = np.fmax(ld, 1.0)
  530. w_ld = np.fmax(w_ld, 1.0)
  531. c = hsq * N / M
  532. het_w = 1.0 / (2 * np.square(intercept + np.multiply(c, ld)))
  533. oc_w = 1.0 / w_ld
  534. w = np.multiply(het_w, oc_w)
  535. return w
  538. class Gencov(LD_Score_Regression):
  539. __null_intercept__ = 0
  541. def __init__(self, z1, z2, x, w, N1, N2, M, hsq1, hsq2, intercept_hsq1, intercept_hsq2,
  542. n_blocks=200, intercept_gencov=None, slow=False, twostep=None):
  543. self.intercept_hsq1 = intercept_hsq1
  544. self.intercept_hsq2 = intercept_hsq2
  545. self.hsq1 = hsq1
  546. self.hsq2 = hsq2
  547. self.N1 = N1
  548. self.N2 = N2
  549. y = z1 * z2
  550. step1_ii = None
  551. if twostep is not None:
  552. step1_ii = np.logical_and(z1**2 < twostep, z2**2 < twostep)
  554. LD_Score_Regression.__init__(self, y, x, w, np.sqrt(N1 * N2), M, n_blocks,
  555. intercept=intercept_gencov, slow=slow, step1_ii=step1_ii)
  556. self.p, self.z = p_z_norm(self.tot, self.tot_se)
  557. self.mean_z1z2 = np.mean(np.multiply(z1, z2))
  559. def summary(self, ref_ld_colnames, P=None, K=None):
  560. '''Print summary of the LD Score regression.'''
  561. out = []
  562. if P is not None and K is not None and\
  563. all((i is not None for i in P)) and all((i is not None for i in K)):
  564. T = 'Liability'
  565. c = gencov_obs_to_liab(1, P[0], P[1], K[0], K[1])
  566. else:
  567. T = 'Observed'
  568. c = 1
  570. out.append('Total ' + T + ' scale gencov: ' +
  571. s(c * self.tot) + ' (' + s(c * self.tot_se) + ')')
  572. if self.n_annot > 1:
  573. out.append('Categories: ' + str(' '.join(ref_ld_colnames)))
  574. out.append(T + ' scale gencov: ' + s(c * self.cat))
  575. out.append(T + ' scale gencov SE: ' + s(c * self.cat_se))
  576. out.append('Proportion of SNPs: ' + s(self.M_prop))
  577. out.append('Proportion of gencov: ' + s(self.prop))
  578. out.append('Enrichment: ' + s(self.enrichment))
  580. out.append('Mean z1*z2: ' + s(self.mean_z1z2))
  581. if self.constrain_intercept:
  582. out.append(
  583. 'Intercept: constrained to {C}'.format(C=s(self.intercept)))
  584. else:
  585. out.append(
  586. 'Intercept: ' + s(self.intercept) + ' (' + s(self.intercept_se) + ')')
  588. return remove_brackets('\n'.join(out))
  590. def _update_func(self, x, ref_ld_tot, w_ld, N, M, Nbar, intercept=None, ii=None):
  591. '''
  592. Update function for IRWLS
  593. x is the output of np.linalg.lstsq.
  594. x[0] is the regression coefficients
  595. x[0].shape is (# of dimensions, 1)
  596. the last element of x[0] is the intercept.
  598. '''
  599. rho_g = M * x[0][0] / Nbar
  600. if intercept is None: # if the regression includes an intercept
  601. intercept = x[0][1]
  603. # remove intercept if we have one
  604. ld = ref_ld_tot[:, 0].reshape(w_ld.shape)
  605. if ii is not None:
  606. N1 = self.N1[ii].reshape((w_ld.shape))
  607. N2 = self.N2[ii].reshape((w_ld.shape))
  608. else:
  609. N1 = self.N1
  610. N2 = self.N2
  612. return self.weights(ld, w_ld, N1, N2, np.sum(M), self.hsq1, self.hsq2, rho_g,
  613. intercept, self.intercept_hsq1, self.intercept_hsq2, ii)
  615. def _update_weights(self, ld, w_ld, sqrt_n1n2, M, rho_g, intercept, ii=None):
  616. '''Weight function with the same signature for Hsq and Gencov.'''
  617. w = self.weights(ld, w_ld, self.N1, self.N2, M, self.hsq1, self.hsq2, rho_g,
  618. intercept, self.intercept_hsq1, self.intercept_hsq2)
  619. return w
  621. @classmethod
  622. def weights(cls, ld, w_ld, N1, N2, M, h1, h2, rho_g, intercept_gencov=None,
  623. intercept_hsq1=None, intercept_hsq2=None, ii=None):
  624. '''
  625. Regression weights.
  627. Parameters
  628. ----------
  629. ld : np.matrix with shape (n_snp, 1)
  630. LD Scores (non-partitioned)
  631. w_ld : np.matrix with shape (n_snp, 1)
  632. LD Scores (non-partitioned) computed with sum r^2 taken over only those SNPs included
  633. in the regression.
  634. M : float > 0
  635. Number of SNPs used for estimating LD Score (need not equal number of SNPs included in
  636. the regression).
  637. N1, N2 : np.matrix of ints > 0 with shape (n_snp, 1)
  638. Number of individuals sampled for each SNP for each study.
  639. h1, h2 : float in [0,1]
  640. Heritability estimates for each study.
  641. rhog : float in [0,1]
  642. Genetic covariance estimate.
  643. intercept : float
  644. Genetic covariance intercept, on the z1*z2 scale (so should be Ns*rho/sqrt(N1*N2)).
  646. Returns
  647. -------
  648. w : np.matrix with shape (n_snp, 1)
  649. Regression weights. Approx equal to reciprocal of conditional variance function.
  651. '''
  652. M = float(M)
  653. if intercept_gencov is None:
  654. intercept_gencov = 0
  655. if intercept_hsq1 is None:
  656. intercept_hsq1 = 1
  657. if intercept_hsq2 is None:
  658. intercept_hsq2 = 1
  660. h1, h2 = max(h1, 0.0), max(h2, 0.0)
  661. h1, h2 = min(h1, 1.0), min(h2, 1.0)
  662. rho_g = min(rho_g, 1.0)
  663. rho_g = max(rho_g, -1.0)
  664. ld = np.fmax(ld, 1.0)
  665. w_ld = np.fmax(w_ld, 1.0)
  666. a = np.multiply(N1, h1 * ld) / M + intercept_hsq1
  667. b = np.multiply(N2, h2 * ld) / M + intercept_hsq2
  668. sqrt_n1n2 = np.sqrt(np.multiply(N1, N2))
  669. c = np.multiply(sqrt_n1n2, rho_g * ld) / M + intercept_gencov
  670. try:
  671. het_w = 1.0 / (np.multiply(a, b) + np.square(c))
  672. except FloatingPointError: # bizarre error; should never happen
  673. raise FloatingPointError('Why did you set hsq intercept <= 0?')
  675. oc_w = 1.0 / w_ld
  676. w = np.multiply(het_w, oc_w)
  677. return w
  680. class RG(object):
  682. def __init__(self, z1, z2, x, w, N1, N2, M, intercept_hsq1=None, intercept_hsq2=None,
  683. intercept_gencov=None, n_blocks=200, slow=False, twostep=None):
  684. self.intercept_gencov = intercept_gencov
  685. self._negative_hsq = None
  686. n_snp, n_annot = x.shape
  687. hsq1 = Hsq(np.square(z1), x, w, N1, M, n_blocks=n_blocks, intercept=intercept_hsq1,
  688. slow=slow, twostep=twostep)
  689. hsq2 = Hsq(np.square(z2), x, w, N2, M, n_blocks=n_blocks, intercept=intercept_hsq2,
  690. slow=slow, twostep=twostep)
  691. gencov = Gencov(z1, z2, x, w, N1, N2, M, hsq1.tot, hsq2.tot, hsq1.intercept,
  692. hsq2.intercept, n_blocks, intercept_gencov=intercept_gencov, slow=slow,
  693. twostep=twostep)
  694. gencov.N1 = None # save memory
  695. gencov.N2 = None
  696. self.hsq1, self.hsq2, self.gencov = hsq1, hsq2, gencov
  697. if (hsq1.tot <= 0 or hsq2.tot <= 0):
  698. self._negative_hsq = True
  699. self.rg_ratio = self.rg = self.rg_se = 'NA'
  700. self.p = self.z = 'NA'
  701. else:
  702. rg_ratio = np.array(
  703. gencov.tot / np.sqrt(hsq1.tot * hsq2.tot)).reshape((1, 1))
  704. denom_delete_values = np.sqrt(
  705. np.multiply(hsq1.tot_delete_values, hsq2.tot_delete_values))
  706. rg = jk.RatioJackknife(
  707. rg_ratio, gencov.tot_delete_values, denom_delete_values)
  708. self.rg_jknife = float(rg.jknife_est)
  709. self.rg_se = float(rg.jknife_se)
  710. self.rg_ratio = float(rg_ratio)
  711. self.p, self.z = p_z_norm(self.rg_ratio, self.rg_se)
  713. def summary(self, silly=False):
  714. '''Print output of Gencor object.'''
  715. out = []
  716. if self._negative_hsq:
  717. out.append('Genetic Correlation: nan (nan) (h2 out of bounds) ')
  718. out.append('Z-score: nan (nan) (h2 out of bounds)')
  719. out.append('P: nan (nan) (h2 out of bounds)')
  720. out.append('WARNING: One of the h2\'s was out of bounds.')
  721. out.append(
  722. 'This usually indicates a data-munging error ' +
  723. 'or that h2 or N is low.')
  724. elif (self.rg_ratio > 1.2 or self.rg_ratio < -1.2) and not silly:
  725. out.append('Genetic Correlation: nan (nan) (rg out of bounds) ')
  726. out.append('Z-score: nan (nan) (rg out of bounds)')
  727. out.append('P: nan (nan) (rg out of bounds)')
  728. out.append('WARNING: rg was out of bounds.')
  729. if self.intercept_gencov is None:
  730. out.append(
  731. 'This often means that h2 is not significantly ' +
  732. 'different from zero.')
  733. else:
  734. out.append(
  735. 'This often means that you have constrained' +
  736. ' the intercepts to the wrong values.')
  737. else:
  738. out.append(
  739. 'Genetic Correlation: ' + s(self.rg_ratio) +
  740. ' (' + s(self.rg_se) + ')')
  741. out.append('Z-score: ' + s(self.z))
  742. out.append('P: ' + s(self.p))
  743. return remove_brackets('\n'.join(out))

regressions.py at commit 2fdeeb3, under GPL-3.0 · at the source

Overview

Authors: Dipti Debnath1,2, Mohammad Housini2, Sanjeev Sariya3, Nicole R. Phillips1,2, Gita A. Pathak3,4, Robert C. Barber2
ORCID iDs: Dipti Debnath, Mohammad Housini, Sanjeev Sariya, Gita A. Pathak, Robert C. Barber
  1. Department of Microbiology, Immunology & Genetics, University of North Texas Health Science Center, 3500 Camp Bowie Blvd, Fort Worth, TX 76107, USA
  2. Institute for Translational Research, University of North Texas Health Science Center, 3500 Camp Bowie Blvd, Fort Worth, TX 76107, USA
  3. Institute for Genomic Health, Icahn School of Medicine at Mount Sinai, 1 Gustave L. Levy Pl, New York, NY 10029, USA
  4. Department of Genetics and Genomic Sciences, Icahn School of Medicine at Mount Sinai, 1 Gustave L. Levy Pl, New York, NY 10029, USA
Institutions: University of North Texas Health Science Center (United States); Icahn School of Medicine at Mount Sinai (United States)
Journal: International journal of molecular sciences, volume 27, issue 6, article 2900
Dates: received 28 January 2026; accepted 17 March 2026; published online 23 March 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.3390/ijms27062900 · PMID 41898760 · PMCID PMC13026669 · OpenAlex W7140146012
Open access: gold, a free copy (OpenAlex)
Status: code verified (how the code is found and verified)
Categories: genetics / omics (modality), human (organism), other condition (population), Alzheimer's / dementia (population), cellular / molecular (subfield)
Methods: Statistics
Keywords: apoptosis, APOE, local genetic correlation, genetically regulated gene expression, inverse comorbidity
MeSH: Alzheimer Disease*, Genetic Predisposition to Disease*, Neoplasms*, Gene Expression Regulation, Neoplastic, Genome-Wide Association Study, Humans, Polymorphism, Single Nucleotide (* major topic)
Topic: Cancer-related cognitive impairment studies (Oncology, Medicine), according to OpenAlex
Citations: cited by 2 papers (Europe PMC); 124 references in the paper

Abstract

Alzheimer’s disease (AD) and cancer are both age-related conditions, yet numerous large-scale epidemiological studies have consistently documented an inverse association, with individuals diagnosed with cancer exhibiting a reduced risk of AD and vice versa. Although this relationship has been replicated across diverse populations, its biological basis remains poorly understood. To address this gap, the present study applies a framework that integrates locus-level genetic correlation (rg) with genetically regulated gene expression to clarify the molecular factors contributing to the inverse epidemiological patterns observed between the two diseases. We used the largest available genome-wide association studies (GWAS) (Nmax = 448,150) to quantify local genetic correlations between AD and several age-associated cancers, including breast, prostate, lung, colorectal, melanoma, basal cell carcinoma, bladder, and endometrial cancer. Eight genomic regions showed significant negative local rg, at the 19q13.31–19q13.32 locus demonstrating strong negative correlations across multiple cancers, including breast, prostate, lung, melanoma, and endometrial cancer. To evaluate the contribution of genetically regulated gene expression, we conducted transcriptome-wide association studies (TWAS) using precomputed gene expression weights from cancer tissues (The Cancer Genome Atlas-TCGA), disease-agnostic tissues (Genotype-Tissue Expression-GTEx), and brain tissue (dorsolateral prefrontal cortex-DLPFC). For each AD–cancer pair, we prioritized genes that were nominally significant in both traits (p < 0.05) and exhibited inverse TWAS Z scores. This analysis identified 24 genes with opposite effect directions between AD and at least three cancer types. TWAS signals also aligned with local rg findings at the 19q13.31–19q13.32 region, suggesting that regulatory variation near this locus contributes to shared but opposing genetic effects beyond the canonical APOE signal. Across cancer types, genes inversely associated with AD converged on pathways involved in cell cycle regulation, apoptosis, DNA-damage response, and metabolic processes. These results support the hypothesis that biological mechanisms promoting proliferation and survival in cancer may oppose those contributing to neurodegeneration in AD.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 1 match between paragraphs and lines of code.

bulik/ldsc

License: GPL-3.0
State: the link answers, verified on 26 October 2026
Evidence: files inventoried
Commit: 2fdeeb3b44379408794154993dbd6101b8946b7e, 16 January 2026
Languages: Python (19), MATLAB (4), Perl (1), R (1)
Size: 1,093 files, 25 scripts
Software Heritage: archived
Found in: the text, “4.3. Global Genetic Correlation Analysis”
Holds: README, license file, environment (environment.yml, requirements.txt, setup.py), tests
Not found: CITATION.cff, continuous integration, documentation
Tools: NumPy (18 files), pandas (11 files), SciPy (5 files), BEDTools (2 files)
Availability: 2 checks, the latest on 26 October 2026: the link answers
  • 26 October 2026: the link answers
  • 26 September 2026: the link answers
27 files
At the source: github.com/bulik/ldsc

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 25 scripts, each with its path and the digest of its content;
  • 1 match between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data Availability Statement

The original contributions presented in this study are included in the article/Supplementary Materials. Further inquiries can be directed to the corresponding authors.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 30 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 6 authors, 5 keywords, 7 MeSH terms, 123 references.

Cite

This paper

Debnath, D., Housini, M., Sariya, S., Phillips, N. R., Pathak, G. A., & Barber, R. C. (2026). Locus- and Gene-Level Insights into the Inverse Association Between Alzheimer's Disease and Cancer. International journal of molecular sciences, 27(6), 2900. https://doi.org/10.3390/ijms27062900

BibTeX

@article{debnath2026locus,
author = {Debnath, Dipti and Housini, Mohammad and Sariya, Sanjeev and Phillips, Nicole R. and Pathak, Gita A. and Barber, Robert C.},
title = {{Locus- and Gene-Level Insights into the Inverse Association Between Alzheimer's Disease and Cancer}},
journal = {International journal of molecular sciences},
year = {2026},
month = mar,
volume = {27},
number = {6},
pages = {2900},
publisher = {Multidisciplinary Digital Publishing Institute (MDPI)},
issn = {1422-0067},
doi = {10.3390/ijms27062900},
url = {https://doi.org/10.3390/ijms27062900},
pmid = {41898760},
pmcid = {PMC13026669}
}

RIS

TY - JOUR
AU - Debnath, Dipti
AU - Housini, Mohammad
AU - Sariya, Sanjeev
AU - Phillips, Nicole R.
AU - Pathak, Gita A.
AU - Barber, Robert C.
TI - Locus- and Gene-Level Insights into the Inverse Association Between Alzheimer's Disease and Cancer
T2 - International journal of molecular sciences
J2 - Int J Mol Sci
PY - 2026
DA - 2026/03/23
VL - 27
IS - 6
SP - 2900
SN - 1422-0067
PB - Multidisciplinary Digital Publishing Institute (MDPI)
DO - 10.3390/ijms27062900
UR - https://doi.org/10.3390/ijms27062900
LA - en
ER -

CSL-JSON

{
"id": "10.3390/ijms27062900",
"type": "article-journal",
"title": "Locus- and Gene-Level Insights into the Inverse Association Between Alzheimer's Disease and Cancer",
"container-title": "International journal of molecular sciences",
"author": [
{
"family": "Debnath",
"given": "Dipti"
},
{
"family": "Housini",
"given": "Mohammad"
},
{
"family": "Sariya",
"given": "Sanjeev"
},
{
"family": "Phillips",
"given": "Nicole R."
},
{
"family": "Pathak",
"given": "Gita A."
},
{
"family": "Barber",
"given": "Robert C."
}
],
"container-title-short": "Int J Mol Sci",
"volume": "27",
"issue": "6",
"page": "2900",
"DOI": "10.3390/ijms27062900",
"PMID": "41898760",
"PMCID": "PMC13026669",
"ISSN": "1422-0067",
"publisher": "Multidisciplinary Digital Publishing Institute (MDPI)",
"URL": "https://doi.org/10.3390/ijms27062900",
"language": "en",
"issued": {
"date-parts": [
[
2026,
3,
23
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1371/journal.pcbi.1012356 [code]
Regression convolutional neural network models implicate peripheral immune regulatory variants in the predisposition to Alzheimer's disease.
Journal: PLoS computational biology
In common: BEDTools, pandas, SciPy, 1 other tool, Alzheimer's / dementia, 8 references
[2] doi:10.3389/fnmol.2022.948456 [code]
Cell type-specific histone acetylation profiling of Alzheimer's disease subjects and integration with genetics.
Journal: Frontiers in molecular neuroscience
In common: pandas, SciPy, NumPy, Alzheimer's / dementia, genetics / omics, cellular / molecular, 9 references
[3] doi:10.1016/j.isci.2023.107958 [code]
Identification of region-specific splicing QTLs in human hippocampal tissue and its distinctive role in brain disorders.
Journal: iScience
In common: BEDTools, pandas, SciPy, 1 other tool, genetics / omics, cellular / molecular, 7 references
[4] doi:10.1038/s42003-026-10030-4 [code]
Cell-type-aware transcriptome-wide association studies identify 91 independent risk genes for Alzheimer's disease dementia.
Journal: Communications biology
In common: Alzheimer's / dementia, genetics / omics, cellular / molecular, 9 references
[5] doi:10.1186/s13073-023-01248-6 [code]
DeepGAMI: deep biologically guided auxiliary learning for multimodal integration and imputation to improve genotype-phenotype prediction.
Journal: Genome medicine
In common: pandas, SciPy, NumPy, Alzheimer's / dementia, genetics / omics, cellular / molecular, 7 references
[6] doi:10.3390/ijms24021428 [code]
Molecular Landscape of Tourette's Disorder.
Journal: International journal of molecular sciences
In common: BEDTools, pandas, SciPy, 1 other tool, other condition, cellular / molecular, 6 references
[7] doi:10.1038/s41562-024-02051-y [code]
The shared genetic architecture and evolution of human language and musical rhythm.
Journal: Nature human behaviour
In common: BEDTools, pandas, SciPy, 1 other tool, genetics / omics, other condition, 6 references
[8] doi:10.1038/s41588-025-02189-z [code]
Genome-wide analyses identify 30 loci associated with obsessive-compulsive disorder.
Journal: Nature genetics
In common: BEDTools, pandas, SciPy, 1 other tool, genetics / omics, other condition, cellular / molecular, 6 references
[9] doi:10.1038/s41467-026-75193-4 [code]
Multi-ancestry gene expression models amplify transcriptome-wide association study discovery and validation.
Journal: Nature communications
In common: pandas, SciPy, NumPy, genetics / omics, cellular / molecular, 7 references
[10] doi:10.1038/s41467-024-52025-x [code]
Genetic, transcriptomic, histological, and biochemical analysis of progressive supranuclear palsy implicates glial activation and novel risk genes.
Journal: Nature communications
In common: pandas, NumPy, Alzheimer's / dementia, genetics / omics, cellular / molecular, 8 references

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Add the badge to the code's README

OSCR: the badge of this paper

If you are an author or a maintainer of this paper's code, you can add this badge. It says what OSCR knows of the paper (its code, its validated tracing map) and links the code to this page. Paste one of these snippets under the README's title. Keep the two marks around it, so that the badge can be replaced in place and is never added twice. OSCR changes nothing in your repository.

Markdown
<!-- oscr-badge:start -->
[![OSCR](https://openscicode.org/badge/doi/10.3390/ijms27062900.svg)](https://openscicode.org/badge/doi/10.3390/ijms27062900)
<!-- oscr-badge:end -->
reStructuredText
.. oscr-badge:start

.. image:: https://openscicode.org/badge/doi/10.3390/ijms27062900.svg
   :target: https://openscicode.org/badge/doi/10.3390/ijms27062900
   :alt: OSCR

.. oscr-badge:end
HTML
<!-- oscr-badge:start -->
<a href="https://openscicode.org/badge/doi/10.3390/ijms27062900"><img src="https://openscicode.org/badge/doi/10.3390/ijms27062900.svg" alt="OSCR"></a>
<!-- oscr-badge:end -->

On GitHub, open the README, choose Edit (the pencil), paste the snippet under the title, then "Commit changes".

If you choose to, you can have the badge kept up to date for you, with oscr badge --write from the command line or with the GitHub Action in your repository's workflows. Every other form and option (a card, AsciiDoc, LaTeX, shields.io) is on the badge page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: code errors and mismatches between the code and the paper about this paper and its code, from readers and its authors, beside the paper and its tracing map.

Report a mismatch between the code and the paper · Report a code error (signed in)

Reproductions: reports from readers who ran the authors' code and did not get the paper's result, what they ran, with which environment, commit and data.

Report a reproduction failure (signed in)

Versions of its code: the releases of its code tied to a version of this paper (the preprint, the accepted manuscript, the version of record…), each with the tracing map versioned with it.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.