OSCR

Cerebral Cortical Structural Variation and General Cognitive Ability: Evidence From Mendelian Randomization.

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] § Materials and Methods › Mendelian Randomization ↔ R/gsmr.R, lines 392–489 · score 0.59 · Mendelian Randomization, genome wide association, Linkage, GWAS, GSMR, instrumental

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

R · 542 lines · 26 KB · no license · 1 match

  1. # gsmr: A tool for GSMR and HEIDI analysis
  2. # The gsmr package Perform Generalized Summary-data-based Mendelian Randomization analysis (GSMR)
  3. # and HEterogeneity In Dependent Instruments analysis to remove pleiotropic outliers (HEIDI-outlier).
  4. # @author Zhihong Zhu <[email hidden]>
  5. # @author Zhili Zheng <[email hidden]>
  6. # @author Futao Zhang <[email hidden]>
  7. # @author Jian Yang <[email hidden]>
  9. eps = 1e-6;
  11. # ************************************************** #
  12. # check data is missing or not #
  13. # ************************************************** #
  14. check_element <- function(vals, argue) {
  15. vals = vals[which(is.finite(vals))]
  16. if(length(vals)==0) {
  17. stop("None of the ", argue, " is found. Please check.")
  18. }
  19. }
  21. # ************************************************** #
  22. # variance-covariane matrix of bXY #
  23. # ************************************************** #
  24. cov_bXY <- function(bzx, bzx_se, bzy, bzy_se, ldrho) {
  25. bXY = bzy / bzx
  26. zscoreZX = bzx / bzx_se
  27. nsnp = dim(ldrho)[1]
  28. covbXY = diag(nsnp)
  29. if(nsnp>1) {
  30. zszxij = zscoreZX%*%t(zscoreZX)
  31. bxyij = bXY%*%t(bXY)
  32. sezyij = bzy_se%*%t(bzy_se)
  33. bzxij = bzx%*%t(bzx)
  34. covbXY = ldrho*sezyij/bzxij + ldrho*bxyij/zszxij
  35. }
  36. return(covbXY)
  37. }
  39. # ************************************************** #
  40. # HEIDI test #
  41. # ************************************************** #
  42. #' @importFrom survey pchisqsum
  43. heidi <- function(bzx, bzx_se, bzx_pval, bzy, bzy_se, ldrho,
  44. heidi_thresh = pchisq(10, 1, lower.tail=F),
  45. nSNPs_thresh=10, maxid=integer(0) ) {
  47. remain_index = seq(1, length(bzx))
  48. m = length(remain_index)
  50. # recaculate the zscore_ZX for vector has been changed
  51. zscore_ZX = bzx / bzx_se
  52. bXY = bzy / bzx
  53. seSMR = sqrt( (bzy_se^2*bzx^2 + bzx_se^2*bzy^2) / bzx^4 )
  54. # remap the maxid according to filtered data
  55. maxid = which(remain_index==maxid)
  56. if(length(maxid) != 1) {
  57. stop("The top SNP for the HEIDI analysis is missing.")
  58. }
  59. # diff = bXY_top - bXY_-i
  60. dev = bXY[maxid] - bXY[-maxid]
  61. # v matrix
  62. covbXY = cov_bXY(bzx, bzx_se, bzy, bzy_se, ldrho)
  63. tmp1 = diag(covbXY)[maxid]
  64. tmp2 = covbXY[-maxid, -maxid]
  65. tmp3 = covbXY[maxid, -maxid]
  66. vdev = tmp1 + tmp2 - tmp3
  67. vdev = t(t(vdev) - tmp3)
  68. diag(vdev) = diag(vdev) + eps
  69. # variance of diff
  70. vardev = diag(covbXY)[-maxid] + tmp1 - 2*tmp3
  71. vardev = vardev + eps
  72. chisq_dev = dev^2 / vardev
  73. if(m>2) {
  74. # correlation matrix of diff
  75. corr_dev = diag(m-1)
  76. # more than 2 instruments
  77. for( i in 1 : (m-2) ) {
  78. for( j in (i+1) : (m-1) ) {
  79. corr_dev[i,j] = corr_dev[j,i] =
  80. vdev[i,j] / sqrt(vdev[i,i]*vdev[j,j])
  81. }
  82. }
  84. # estimate the p value
  85. lambda = eigen(corr_dev, symmetric=TRUE, only.values=TRUE)$values
  86. t = length(lambda)
  87. pHet = pchisqsum(sum(chisq_dev)+eps, df=rep(1,t), a=lambda, method="sadd", lower.tail=F)
  89. } else {
  90. pHet = pchisq(chisq_dev, 1, lower.tail=F)
  91. }
  93. return(list(pheidi=pHet, nsnps=m, used_index=remain_index))
  94. }
  96. # ************************************************** #
  97. # Iterations for HEIDI-ouliter #
  98. # ************************************************** #
  100. heidi_outlier_iter <- function(bzx, bzx_se, bzx_pval, bzy, bzy_se, ldrho, gwas_thresh, heidi_thresh, remain_index) {
  101. # remove pleiotropic instruments
  102. m = length(remain_index)
  103. # grab reference instrument
  104. bxy = bzy/bzx
  105. bxy_q = quantile(bxy, probs = seq(0, 1, 0.2))
  106. min_bxy = as.numeric(bxy_q[3]); max_bxy = as.numeric(bxy_q[4]);
  107. slctindx = which(bxy <= max_bxy & bxy >= min_bxy)
  108. if(length(slctindx)==0) {
  109. stop("The top SNP for the HEIDI-outlier analysis is missing. None SNPs are in the third quintile of the distribution of bxy.");
  110. }
  111. refid = slctindx[which.min(bzx_pval[slctindx])]
  112. pheidi = as.numeric()
  113. for( i in 1 : m ) {
  114. if( i==refid ) next
  115. heidi_result = heidi(bzx[c(refid,i)], bzx_se[c(refid,i)], bzx_pval[c(refid,i)],
  116. bzy[c(refid,i)], bzy_se[c(refid,i)],
  117. ldrho[c(refid,i), c(refid,i)],
  118. gwas_thresh, 2, 1)
  119. pheidi[i] = as.numeric(heidi_result$pheidi)
  120. }
  121. remain_index = remain_index[sort(c(refid,which(pheidi>=heidi_thresh)))]
  122. return(remain_index)
  123. }
  126. # ************************************************** #
  127. # Test identical elements #
  128. # ************************************************** #
  129. check_vec_elements_eq <- function(dim_vec){
  130. return(all(dim_vec == dim_vec[1]))
  131. }
  134. # ************************************************** #
  135. # LD pruning #
  136. # ************************************************** #
  137. # LD pruning, removing pairs of SNPs that LD r2 > threshold
  138. # return index: the index be used for further analysis
  139. snp_ld_prune = function(ldrho, ld_r2_thresh) {
  140. # initialization
  141. nsnp = dim(ldrho)[1]
  142. diag(ldrho) = 0
  143. ldrho[upper.tri(ldrho)]=0
  144. include_id = c(1:nsnp)
  146. # save the index which have ld r^2 > threshold
  147. indx = which(ldrho^2>ld_r2_thresh, arr.ind=T)
  148. if(length(indx)==0) return(NULL);
  149. indx1 = as.numeric(indx[,2])
  150. indx2 = as.numeric(indx[,1])
  151. slct_indx = c(indx1, indx2)
  152. indxbuf = unique(sort(slct_indx))
  154. # count how many SNPs in high LD
  155. nproc = length(indxbuf)
  156. n_slct_snp = as.numeric()
  157. for( i in 1 : nproc) {
  158. n_slct_snp[i] = length(which(slct_indx==indxbuf[i]))
  159. }
  161. # decide the index to remove
  162. nproc = length(indx1)
  163. if(nproc==0) return(NULL);
  164. for( i in 1 : nproc ) {
  165. n1 = n_slct_snp[which(indxbuf==indx1[i])]
  166. n2 = n_slct_snp[which(indxbuf==indx2[i])]
  167. if(n1 < n2) {
  168. t = indx1[i]; indx1[i] = indx2[i]; indx2[i] = t
  169. }
  170. }
  171. indx1 = unique(sort(indx1))
  172. return(indx1)
  173. }
  176. # ************************************************** #
  177. # Filter the summary data #
  178. # ************************************************** #
  179. # filter the input to se and zscore
  180. # parameters see HEIDI or SMR_Multi functions
  181. # return index: the index be used for further analysis
  182. # NOTICE: the parameters will be revised in place after run!!!!! TAKE CARE
  183. filter_summdat <- function(snp_id, bzx, bzx_se, bzx_pval, bzy, bzy_se, ldrho, n_ref, nsnps_thresh, pvalue_thresh, ld_r2_thresh, fdr_thresh){
  184. # make sure they are numeric numbers
  185. m = length(bzx)
  186. message("There are ", m, " SNPs in the dataset.")
  187. remain_index = seq(1,m)
  188. snpIDp = as.character(snp_id);
  189. bZXp = as.numeric(as.character(bzx))
  190. seZXp = as.numeric(as.character(bzx_se));
  191. pZXp = as.numeric(as.character(bzx_pval));
  192. bZYp = as.numeric(as.character(bzy))
  193. seZYp = as.numeric(as.character(bzy_se));
  194. ldrhop = matrix(as.numeric(as.character(unlist(ldrho))), m, m)
  196. # remove SNPs with missing value
  197. na_snps=c()
  198. indx = which(!is.finite(bZXp) | !is.finite(seZXp) | !is.finite(pZXp) | !is.finite(bZYp) | !is.finite(seZYp) |
  199. is.na(bZXp) | is.na(seZXp) | is.na(pZXp) | is.na(bZYp) | is.na(seZYp))
  200. if(length(indx)>0) {
  201. na_snps = snpIDp[indx];
  202. bZXp = bZXp[-indx]; seZXp = seZXp[-indx]; pZXp = pZXp[-indx];
  203. bZYp = bZYp[-indx]; seZYp = seZYp[-indx];
  204. ldrhop = ldrhop[-indx, -indx];
  205. remain_index = remain_index[-indx]
  206. snpIDp = snpIDp[-indx]
  207. warning(length(indx), " SNPs were removed due to missing estimates in the summary data.")
  208. }
  209. m = length(bZXp)
  210. if(m < nsnps_thresh) stop("At least ", nsnps_thresh, " SNPs are required. Note: this hard limit can be changed by the \"nsnps_thresh\".");
  212. # remove SNPs with very small SE
  213. indx = which(seZXp < eps | seZYp < eps )
  214. if(length(indx)>0) {
  215. na_snps = c(na_snps, snpIDp[indx]);
  216. bZXp = bZXp[-indx]; seZXp = seZXp[-indx]; pZXp = pZXp[-indx];
  217. bZYp = bZYp[-indx]; seZYp = seZYp[-indx];
  218. ldrhop = ldrhop[-indx, -indx];
  219. remain_index = remain_index[-indx]
  220. snpIDp = snpIDp[-indx]
  221. warning(length(indx), " SNPs were removed due to extremely small standard error. Please check that data.")
  222. }
  223. m = length(bZXp)
  224. if(m < nsnps_thresh) stop("At least ", nsnps_thresh, " SNPs are required. Note: this hard limit can be changed by the \"nsnps_thresh\".");
  225. # remove SNPs with missing LD
  226. indx = which(is.na(ldrhop[upper.tri(ldrhop)]))
  227. if(length(indx) > 0)
  228. stop("LD correlations between ", length(indx), " pairs of SNPs are missing. Please check the MAF of the SNPs and the missingness rate in the reference sample.")
  230. # z score of bzx
  231. weak_snps = c()
  232. indx = which(pZXp > pvalue_thresh)
  233. if(length(indx)>0) {
  234. weak_snps = snpIDp[indx];
  235. bZXp = bZXp[-indx]; seZXp = seZXp[-indx]; pZXp = pZXp[-indx];
  236. bZYp = bZYp[-indx]; seZYp = seZYp[-indx];
  237. ldrhop = ldrhop[-indx, -indx];
  238. remain_index = remain_index[-indx];
  239. snpIDp = snpIDp[-indx];
  240. warning(length(indx), " non-significant SNPs were removed.")
  241. }
  242. m = length(bZXp)
  243. if(m < nsnps_thresh) stop("At least ", nsnps_thresh, " SNPs are required. Note: this hard limit can be changed by the \"nsnps_thresh\".");
  245. # check LD r
  246. linkage_snps = c()
  247. indx = snp_ld_prune(ldrhop, ld_r2_thresh)
  248. if(length(indx) > 0) {
  249. linkage_snps = snpIDp[indx]
  250. bZXp = bZXp[-indx]; seZXp = seZXp[-indx]; pZXp = pZXp[-indx];
  251. bZYp = bZYp[-indx]; seZYp = seZYp[-indx];
  252. ldrhop = ldrhop[-indx, -indx];
  253. remain_index = remain_index[-indx]
  254. snpIDp = snpIDp[-indx]
  255. warning("There were SNPs in high LD. After LD pruning with a LD r2 threshold of ", ld_r2_thresh, ", ", length(indx), " SNPs were removed. Note: The threshold of LD can be changed by the \"ld_r2_thresh\".")
  256. }
  257. m = length(bZXp)
  258. if(m < nsnps_thresh) stop("At least ", nsnps_thresh, " SNPs are required. Note: this hard limit can be changed by the \"nsnps_thresh\".");
  260. # update LD correlation matrix
  261. var_rho = 1/n_ref
  262. pval_rho = pchisq(ldrhop[upper.tri(ldrhop)]^2/var_rho, 1, lower.tail=F)
  263. qval_rho = p.adjust(pval_rho, method = "fdr")
  264. qval_mat = matrix(0, m, m)
  265. qval_mat[upper.tri(qval_mat)] = qval_rho; qval_mat = t(qval_mat); qval_mat[upper.tri(qval_mat)] = qval_rho;
  266. rho_index = which(qval_mat >= fdr_thresh, arr.ind=T)
  267. ldrhop[rho_index] = 0
  269. message(length(remain_index), " SNPs were retained after filtering.")
  271. # replace the parameters in place
  272. eval.parent(substitute(bzx <- bZXp))
  273. eval.parent(substitute(bzy <- bZYp))
  274. eval.parent(substitute(bzx_se <- seZXp))
  275. eval.parent(substitute(bzy_se <- seZYp))
  276. eval.parent(substitute(bzx_pval <- pZXp))
  277. eval.parent(substitute(ldrho <- ldrhop))
  278. return(list(remain_index=remain_index, na_snps=na_snps, weak_snps=weak_snps, linkage_snps=linkage_snps))
  279. }
  281. # ************************************************** #
  282. # standardization of b and s.e. #
  283. # ************************************************** #
  284. #' @title Standardization of effect size and its standard error
  285. #' @description Standardization of SNP effect and its standard error using z-statistic, allele frequency and sample size
  286. #' @usage std_effect(snp_freq, b, se, n)
  287. #' @param snp_freq vector, allele frequencies
  288. #' @param b vector, SNP effects on risk factor
  289. #' @param se vector, standard errors of b
  290. #' @param n vector, per-SNP sample sizes for GWAS of the risk factor
  291. #' @examples
  292. #' data("gsmr")
  293. #' std_effects = std_effect(gsmr_data$a1_freq, gsmr_data$bzx, gsmr_data$bzx_se, gsmr_data$bzx_n)
  294. #'
  295. #' @return Standardised effect (b) and standard error (se)
  296. #' @export
  297. std_effect <- function(snp_freq, b, se, n) {
  298. # double check the counts
  299. # if data set is empty
  300. check_element(snp_freq,"Minor Allele Frequency")
  301. check_element(b,"Effect Size")
  302. check_element(se,"Standard Error")
  303. check_element(n,"Sample Size")
  305. message("std effect: ", length(b), " instruments loaded.")
  307. # length is different or not
  308. len_vec <- c(length(snp_freq),length(b),length(se),length(n))
  309. if (!check_vec_elements_eq(len_vec)){
  310. stop("Lengths of the input vectors are different. Please check.");
  311. }
  313. # check missing values
  314. indx = which(!is.finite(snp_freq) | !is.finite(b) | !is.finite(se) | !is.finite(n) | is.na(snp_freq) | is.na(b) | is.na(se) | is.na(n))
  315. if (length(indx)>0) {
  316. stop("There are ", length(indx), " SNPs with missing estimates in the summary data. Please check.");
  317. }
  319. # make sure they are numeric numbers
  320. snpfreq = as.numeric(as.character(snp_freq))
  321. b = as.numeric(as.character(b))
  322. se = as.numeric(as.character(se))
  323. n = as.numeric(as.character(n))
  325. zscore = b / se
  326. b_p = zscore / sqrt(2*snp_freq*(1-snp_freq)*(n+zscore^2))
  327. se_p = 1 / sqrt(2*snp_freq*(1-snp_freq)*(n+zscore^2))
  329. return(list(b=b_p,se=se_p))
  330. }
  332. # ************************************************** #
  333. # HEIDI-outlier analysis #
  334. # ************************************************** #
  335. #' @title HEIDI-outlier analysis
  336. #' @description An analysis to detect and eliminate from the analysis instruments that show significant pleiotropic effects on both risk factor and disease
  337. #' @usage heidi_outlier(bzx, bzx_se, bzx_pval, bzy, bzy_se, ldrho, snpid, n_ref, gwas_thresh=5e-8, heidi_outlier_thresh=0.01, nsnps_thresh=10, ld_fdr_thresh=0.05)
  338. #' @param bzx vector, SNP effects on risk factor
  339. #' @param bzx_se vector, standard errors of bzx
  340. #' @param bzx_pval vector, p values for bzx
  341. #' @param bzy vector, SNP effects on disease
  342. #' @param bzy_se vector, standard errors of bzy
  343. #' @param ldrho LD correlation matrix of the SNPs
  344. #' @param snpid genetic instruments
  345. #' @param n_ref sample size of the reference sample
  346. #' @param gwas_thresh threshold p-value to select instruments from GWAS for risk factor
  347. #' @param heidi_outlier_thresh threshold p-value to remove pleiotropic outliers (the default value is 0.01)
  348. #' @param nsnps_thresh the minimum number of instruments required for the GSMR analysis (we do not recommend users to set this number smaller than 10)
  349. #' @param ld_r2_thresh LD r2 threshold to remove SNPs in high LD
  350. #' @param ld_fdr_thresh FDR threshold to remove the chance correlations between SNP instruments
  351. #' @examples
  352. #' data("gsmr")
  353. #' filtered_index = heidi_outlier(gsmr_data$bzx, gsmr_data$bzx_se, gsmr_data$bzx_pval, gsmr_data$bzy, gsmr_data$bzy_se, ldrho, gsmr_data$SNP, n_ref, 5e-8, 0.01, 10, 0.1, 0.05)
  354. #'
  355. #' @return Retained index of genetic instruments, SNPs with missing values, with non-significant p-values and those in LD.
  356. #' @export
  357. heidi_outlier <- function(bzx, bzx_se, bzx_pval, bzy, bzy_se, ldrho, snpid,
  358. n_ref, gwas_thresh=5e-8, heidi_outlier_thresh=0.01, nsnps_thresh=10, ld_r2_thresh = 0.1, ld_fdr_thresh=0.05) {
  359. # Subset of LD r matrix
  360. len1 = length(Reduce(intersect, list(snpid, colnames(ldrho))))
  361. len2 = length(snpid)
  362. len_vec <- c(len1, len2)
  363. if (!check_vec_elements_eq(len_vec)){
  364. stop(paste(len2 - len1, " SNPs are missing in the LD correlation matrix. Please check.", sep=""))
  365. }
  366. ldrho = ldrho[snpid, snpid]
  367. # double check the counts
  368. len_vec <- c(length(bzx),length(bzx_se),length(bzx_pval),length(bzy),length(bzy_se),dim(ldrho)[1],dim(ldrho)[2])
  369. if (!check_vec_elements_eq(len_vec)){
  370. stop("Lengths of the input vectors are different. Please check.");
  371. }
  372. message("HEIDI-outlier: ", length(bzx), " instruments loaded.")
  373. # filter dataset
  374. resbuf <- filter_summdat(snpid, bzx, bzx_se, bzx_pval, bzy, bzy_se, ldrho, n_ref, nsnps_thresh, gwas_thresh, ld_r2_thresh, ld_fdr_thresh)
  375. remain_index <- resbuf$remain_index;
  376. na_snps <- resbuf$na_snps; weak_snps <- resbuf$weak_snps; linkage_snps <- resbuf$linkage_snps
  377. if(length(remain_index) < nsnps_thresh) {
  378. stop("Not enough SNPs for the HEIDI-outlier analysis. At least ", nsnps_thresh, " SNPs are required. Note: this hard limit can be changed by the \"nsnp_thresh\".");
  379. }
  380. # Perform HEIDI-outlier
  381. pleio_snps <- NULL
  382. remain_index_tmp <- remain_index
  383. remain_index <- heidi_outlier_iter(bzx, bzx_se, bzx_pval, bzy, bzy_se, ldrho, gwas_thresh, heidi_outlier_thresh, remain_index)
  384. if(length(remain_index) < length(remain_index_tmp)) {
  385. removed_index <- remain_index_tmp[-match(remain_index, remain_index_tmp)]
  386. pleio_snps <- snpid[removed_index]
  387. }
  388. message(length(remain_index), " SNPs were retained after the HEIDI-outlier analysis.")
  389. return(list(remain_index=remain_index, na_snps=na_snps, weak_snps=weak_snps, linkage_snps=linkage_snps, pleio_snps=pleio_snps))
  390. }
  392. # ************************************************** #
  393. # GSMR analysis #
  394. # ************************************************** #
  395. #' @title Generalized Summary-data-based Mendelian Randomization analysis
  396. #' @description GSMR (Generalised Summary-data-based Mendelian Randomisation) is a flexible and powerful approach that utilises multiple genetic instruments to test for causal association between a risk factor and disease using summary-level data from independent genome-wide association studies.
  397. #' @usage gsmr(bzx, bzx_se, bzx_pval, bzy, bzy_se, ldrho, snpid, heidi_outlier_flag=T, gwas_thresh=5e-8, heidi_outlier_thresh=0.01, nsnps_thresh=10)
  398. #' @param bzx vector, SNP effects on risk factor
  399. #' @param bzx_se vector, standard errors of bzx
  400. #' @param bzx_pval vector, p values for bzx
  401. #' @param bzy vector, SNP effects on disease
  402. #' @param bzy_se vector, standard errors of bzy
  403. #' @param ldrho LD correlation matrix of the SNPs
  404. #' @param snpid genetic instruments
  405. #' @param n_ref sample size of the reference sample
  406. #' @param heidi_outlier_flag flag for HEIDI-outlier analysis
  407. #' @param gwas_thresh threshold p-value to select instruments from GWAS for risk factor
  408. #' @param heidi_outlier_thresh HEIDI-outlier threshold
  409. #' @param nsnps_thresh the minimum number of instruments required for the GSMR analysis (we do not recommend users to set this number smaller than 10)
  410. #' @param ld_r2_thresh LD r2 threshold to remove SNPs in high LD
  411. #' @param ld_fdr_thresh FDR threshold to remove the chance correlations between SNP instruments
  412. #' @examples
  413. #' data("gsmr")
  414. #' gsmr_result = gsmr(gsmr_data$bzx, gsmr_data$bzx_se, gsmr_data$bzx_pval, gsmr_data$bzy, gsmr_data$bzy_se, ldrho, gsmr_data$SNP, n_ref, T, 5e-8, 0.01, 10, 0.1, 0.05)
  415. #'
  416. #' @return Estimate of causative effect of risk factor on disease (bxy), the corresponding standard error (bxy_se), p-value (bxy_pval), SNP index (used_index), SNPs with missing values, with non-significant p-values and those in LD.
  417. #' @export
  418. gsmr <- function(bzx, bzx_se, bzx_pval, bzy, bzy_se,
  419. ldrho, snpid, n_ref, heidi_outlier_flag=T, gwas_thresh=5e-8, heidi_outlier_thresh=0.01, nsnps_thresh=10, ld_r2_thresh=0.1, ld_fdr_thresh=0.05) {
  420. # subset of LD r matrix
  421. len1 = length(Reduce(intersect, list(snpid, colnames(ldrho))))
  422. len2 = length(snpid)
  423. len_vec <- c(len1, len2)
  424. if (!check_vec_elements_eq(len_vec)){
  425. stop(paste(len2 - len1, " SNPs are missing in the LD correlation matrix. Please check.", sep=""))
  426. }
  427. ldrho = ldrho[snpid, snpid]
  428. # double check the counts
  429. len_vec <- c(length(bzx),length(bzx_se),length(bzy),length(bzy_se),dim(ldrho)[1],dim(ldrho)[2])
  430. if (!check_vec_elements_eq(len_vec)){
  431. stop("Lengths of the input vectors are different. Please check.");
  432. }
  433. message("GSMR analysis: ", length(bzx), " instruments loaded.")
  434. # filter dataset
  435. resbuf <- filter_summdat(snpid, bzx, bzx_se, bzx_pval, bzy, bzy_se, ldrho, n_ref, nsnps_thresh, gwas_thresh, ld_r2_thresh, ld_fdr_thresh)
  436. remain_index<-resbuf$remain_index;
  437. na_snps<-resbuf$na_snps; weak_snps<-resbuf$weak_snps; linkage_snps<-resbuf$linkage_snps;
  438. if(length(remain_index) < nsnps_thresh) {
  439. stop("Not enough SNPs for the GSMR analysis. At least ", nsnps_thresh, " SNPs are required. Note: this hard limit can be changed by the \"nsnps_thresh\".");
  440. }
  441. pleio_snps=NULL;
  442. if(heidi_outlier_flag) {
  443. # Perform HEIDI-outlier
  444. remain_index2 <- heidi_outlier_iter(bzx, bzx_se, bzx_pval, bzy, bzy_se, ldrho, gwas_thresh, heidi_outlier_thresh, remain_index)
  445. if(length(remain_index2) < nsnps_thresh) {
  446. stop("Not enough SNPs for the GSMR analysis. At least ", nsnps_thresh, " are required. Note: this hard limit can be changed by \"nsnps_thresh\".");
  447. } else {
  448. message(length(remain_index2), " SNPs were retained after the HEIDI-outlier analysis.")
  449. }
  450. # Save pleiotropic SNPs
  451. if(length(remain_index2) < length(remain_index)) {
  452. removed_index <- remain_index[-match(remain_index2, remain_index)]
  453. pleio_snps <- snpid[removed_index]
  454. }
  455. # Update estimates
  456. remain_index2_tmp = remain_index2
  457. remain_index2 = match(remain_index2, remain_index)
  458. bzx = bzx[remain_index2]; bzx_se = bzx_se[remain_index2];
  459. bzy = bzy[remain_index2]; bzy_se = bzy_se[remain_index2];
  460. ldrho = ldrho[remain_index2,remain_index2];
  461. remain_index <- remain_index2_tmp
  462. }
  463. # do the SMR test with multiple instruments
  464. message("Computing the estimate of bxy at each instrument.")
  465. bXY = bzy/bzx
  467. message("Estimating the variance-covariance matrix for bxy.")
  468. covbXY = cov_bXY(bzx, bzx_se, bzy, bzy_se, ldrho)
  469. diag(covbXY) = diag(covbXY) + eps
  470. # Eigen decomposition
  471. resbuf = eigen(covbXY, symmetric=TRUE)
  472. eval = as.numeric(resbuf$values)
  473. evec = resbuf$vectors
  474. if(min(abs(eval)) < eps) {
  475. stop("The variance-covariance matrix for bxy is not invertible!");
  476. }
  477. covbXY_inv = evec%*%diag(1/eval)%*%t(evec)
  478. message("Estimating bxy using all the instruments.")
  479. vec_1 = rep(1, length(bzx))
  480. num_1_v_1 = as.numeric(solve(t(vec_1)%*%covbXY_inv%*%vec_1))
  481. vec_1_v = as.numeric(t(vec_1)%*%covbXY_inv)
  482. bXY_GLS = num_1_v_1*vec_1_v%*%bXY
  483. varbXY_GLS = num_1_v_1
  484. chisqbXY_GLS = bXY_GLS^2/varbXY_GLS
  485. pbXY_GLS = pchisq(chisqbXY_GLS, 1, lower.tail=F)
  486. message("GSMR analysis is completed.")
  487. return(list(bxy=bXY_GLS, bxy_se=sqrt(varbXY_GLS), bxy_pval=pbXY_GLS, used_index=remain_index,
  488. na_snps=na_snps, weak_snps=weak_snps, linkage_snps=linkage_snps, pleio_snps=pleio_snps))
  489. }
  491. # ************************************************** #
  492. # Bi-directional GSMR analysis #
  493. # ************************************************** #
  494. #' @title Bi-directional GSMR analysis
  495. #' @description Bi-directional GSMR analysis is composed of a forward-GSMR analysis and a reverse-GSMR analysis that uses SNPs associated with the disease (e.g. at < 5e-8) as the instruments to test for putative causal effect of the disease on the risk factor.
  496. #' @usage bi_gsmr(bzx, bzx_se, bzx_pval, bzy, bzy_se, bzy_pval, ldrho, snpid, heidi_outlier_flag=T, gwas_thresh=5e-8, heidi_outlier_thresh=0.01, nsnps_thresh=10)
  497. #' @param bzx vector, SNP effects on risk factor
  498. #' @param bzx_se vector, standard errors of bzx
  499. #' @param bzx_pval vector, p values for bzx
  500. #' @param bzy vector, SNP effects on disease
  501. #' @param bzy_se vector, standard errors of bzy
  502. #' @param bzy_pval vector, p values for bzy
  503. #' @param ldrho LD correlation matrix of the SNPs
  504. #' @param snpid genetic instruments
  505. #' @param n_ref sample size of the reference sample
  506. #' @param heidi_outlier_flag flag for HEIDI-outlier analysis
  507. #' @param gwas_thresh threshold p-value to select instruments from GWAS for risk factor
  508. #' @param heidi_outlier_thresh HEIDI-outlier threshold
  509. #' @param nsnps_thresh the minimum number of instruments required for the GSMR analysis (we do not recommend users to set this number smaller than 10)
  510. #' @param ld_r2_thresh LD r2 threshold to remove SNPs in high LD
  511. #' @param ld_fdr_thresh FDR threshold to remove the chance correlations between SNP instruments
  512. #' @examples
  513. #' data("gsmr")
  514. #' gsmr_result = bi_gsmr(gsmr_data$bzx, gsmr_data$bzx_se, gsmr_data$bzx_pval, gsmr_data$bzy, gsmr_data$bzy_se, gsmr_data$bzy_pval, ldrho, gsmr_data$SNP, n_ref, T, 5e-8, 0.01, 10, 0.1, 0.05)
  515. #'
  516. #' @return Estimate of causative effect of risk factor on disease (forward_bxy), the corresponding standard error (forward_bxy_se), p-value (forward_bxy_pval) and SNP index (forward_index), and estimate of causative effect of disease on risk factor (reverse_bxy), the corresponding standard error (reverse_bxy_se), p-value (reverse_bxy_pval), SNP index (reverse_index), SNPs with missing values, with non-significant p-values and those in LD.
  517. #' @export
  518. bi_gsmr <- function(bzx, bzx_se, bzx_pval, bzy, bzy_se, bzy_pval,
  519. ldrho, snpid, n_ref, heidi_outlier_flag=T, gwas_thresh=5e-8, heidi_outlier_thresh=0.01, nsnps_thresh=10, ld_r2_thresh=0.1, ld_fdr_thresh=0.05) {
  520. ## Forward GSMR
  521. message("Forward GSMR analysis...")
  522. gsmr_result=gsmr(bzx, bzx_se, bzx_pval, bzy, bzy_se, ldrho, snpid, n_ref, heidi_outlier_flag, gwas_thresh, heidi_outlier_thresh, nsnps_thresh, ld_r2_thresh, ld_fdr_thresh)
  523. bxy1 = gsmr_result$bxy; bxy1_se = gsmr_result$bxy_se; bxy1_pval = gsmr_result$bxy_pval;
  524. bxy1_index = gsmr_result$used_index;
  525. na_snps = gsmr_result$na_snps; weak_snps = gsmr_result$weak_snps; linkage_snps = gsmr_result$linkage_snps; pleio_snps = gsmr_result$pleio_snps;
  527. ## Reverse GSMR
  528. message("Reverse GSMR analysis...")
  529. gsmr_result=gsmr(bzy, bzy_se, bzy_pval, bzx, bzx_se, ldrho, snpid, n_ref, heidi_outlier_flag, gwas_thresh, heidi_outlier_thresh, nsnps_thresh, ld_r2_thresh, ld_fdr_thresh)
  530. bxy2 = gsmr_result$bxy; bxy2_se = gsmr_result$bxy_se; bxy2_pval = gsmr_result$bxy_pval;
  531. bxy2_index = gsmr_result$used_index;
  532. na_snps = c(na_snps, gsmr_result$na_snps);
  533. weak_snps = c(weak_snps, gsmr_result$weak_snps);
  534. linkage_snps = c(linkage_snps, gsmr_result$linkage_snps);
  535. pleio_snps = c(pleio_snps, gsmr_result$pleio_snps);
  536. return(list(forward_bxy=bxy1, forward_bxy_se=bxy1_se,
  537. forward_bxy_pval=bxy1_pval, forward_index=bxy1_index,
  538. reverse_bxy=bxy2, reverse_bxy_se=bxy2_se,
  539. reverse_bxy_pval=bxy2_pval, reverse_index=bxy2_index,
  540. na_snps=na_snps, weak_snps=weak_snps, linkage_snps=linkage_snps, pleio_snps=pleio_snps))
  541. }

gsmr.R at commit 53f5034, no license · at the source

Overview

Authors: Chun‐Ju Chou1,2, Mark Fiecas3, Elisabetta C. del Re4,5, Eero Vuoksimaa6, Chi‐Hua Chen1
ORCID iDs: Chun‐Ju Chou, Mark Fiecas, Elisabetta C. del Re, Eero Vuoksimaa, Chi‐Hua Chen
  1. Department of Radiology University of California, San Diego San Diego California USA
  2. Department of Bioengineering University of California, San Diego San Diego California USA
  3. Division of Biostatistics University of Minnesota School of Public Health Minneapolis Minnesota USA
  4. Department of Psychiatry Harvard Medical School Boston Massachusetts USA
  5. College of Science Northeastern University Boston Massachusetts USA
  6. Institute for Molecular Medicine Finland (FIMM), Helsinki Institute of Life Science, University of Helsinki Helsinki Finland
Institutions: University of California San Diego (United States); University of Minnesota (United States); Northeastern University (United States); Harvard University (United States); University of Helsinki (Finland); Institute for Molecular Medicine Finland (Finland); HiLIFE – Elämäntieteiden Instituutti
Journal: Human brain mapping, volume 47, issue 13, article e70635
Dates: received 25 February 2026; accepted 16 August 2026; published online 30 August 2026; in print September 2026
Type: Research article · Language: English
License: CC BY-NC
Identifiers: DOI 10.1002/hbm.70635 · PMID 42669588 · PMCID PMC13526639 · OpenAlex W7204736479
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: genetics / omics (modality), structural MRI / diffusion (modality), human (organism), cognitive (subfield)
Methods: Statistics, Smoothing, state filtering, decompositions
Keywords: cortical morphology, general cognitive ability, genetically‐informed parcellation, genome‐wide association study, Mendelian randomization
MeSH: Aptitude*, Cerebral Cortex*, Cognition*, Individuality*, Aged, Brain Cortical Thickness, Female, Genome-Wide Association Study, Humans, Magnetic Resonance Imaging, Male, Mendelian Randomization Analysis, Middle Aged, UK Biobank (* major topic)
Topic: Cognitive Abilities and Testing (Experimental and Cognitive Psychology, Psychology), according to OpenAlex
Funding: NIH (R01MH132783); Sigrid Jusélius Foundation and the Research Council of Finland (314639, 345988)
Citations: not cited yet (Europe PMC); 55 references in the paper

Abstract

Understanding the cortical architecture underlying individual differences in general cognitive ability (GCA) remains a central question in cognitive neuroscience. Prior work has established associations between global brain size and GCA, yet the regional effects and directionality of these relationships remain debated. Using a genetically informed cortical parcellation in 11,289 UK Biobank participants, we examined associations between cortical surface area (SA), cortical thickness (CT), and GCA measured via verbal–numerical reasoning. Total SA showed a robust positive association with GCA. At the regional level, dorsolateral prefrontal and superior temporal SA exhibited the strongest positive associations, which persisted after adjustment for global SA. In contrast, CT showed comparatively modest associations. Using Mendelian randomization (MR) with genome‐wide significant genetic instruments, we observed evidence consistent with a bidirectional relationship between total SA and GCA. At the regional level, dorsolateral prefrontal and temporal SA demonstrated evidence of MR‐inferred directional effects on GCA, while GCA showed evidence of MR‐inferred directional effects on total SA and perisylvian thickness. These findings support a polyregional SA architecture underlying GCA, with prominent contributions from prefrontal and temporal association cortices. Our results refine global brain–GCA models and highlight the value of genetically informed parcellation for identifying regional cortical contributions.

Reproduced under the paper's license (CC BY-NC), 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.

jianyang-lab/gsmr

License: none: the authors keep all their rights
State: the link answers, verified on 26 September 2026
Evidence: files inventoried
Commit: 53f5034f4ebde2d149930a924449b57f6a36cb85, 13 September 2023
Languages: R (5)
Size: 18 files, 5 scripts
Software Heritage: not archived
Found in: “Data Availability Statement”
Holds: README, environment (DESCRIPTION), documentation, 2 notebooks
Not found: license file, CITATION.cff, tests, continuous integration
Availability: 1 check, the latest on 26 September 2026: the link answers
  • 26 September 2026: the link answers
6 files
At the source: github.com/JianYang-Lab/gsmr/releases

The paper's code and data availability statement is in the Data section.

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;
  • 5 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

This study utilizes individual‐level genetic and imaging data from the UK Biobank (https://www.ukbiobank.ac.uk/). The genome‐wide association data for cortical regions were provided from our previously published studies, which can be accessed via the GWAS Catalog (https://www.ebi.ac.uk/gwas/publications/35113692 and https://www.ebi.ac.uk/gwas/publications/36893272). The genome‐wide association data for intelligence were obtained directly from the authors upon request, as a means of avoiding overlap with UKB samples, and the original source was from the Psychiatric Genomics Consortium (PGC) (https://pgc.unc.edu/) or the GWAS Catalog (https://www.ebi.ac.uk/gwas/publications/29942086). GSMR method and code are publicly available via the GCTA software repository on GitHub (https://github.com/JianYang‐Lab/gsmr/releases (https://github.com/JianYang-Lab/gsmr/releases)).

Reproduced under the paper's license (CC BY-NC), 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, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 5 authors, 5 keywords, 14 MeSH terms, 2 funders, 55 references.

Cite

This paper

Chou, C., Fiecas, M., del Re, E. C., Vuoksimaa, E., & Chen, C. (2026). Cerebral Cortical Structural Variation and General Cognitive Ability: Evidence From Mendelian Randomization. Human brain mapping, 47(13), e70635. https://doi.org/10.1002/hbm.70635

BibTeX

@article{chou2026cerebral,
author = {Chou, Chun‐Ju and Fiecas, Mark and del Re, Elisabetta C. and Vuoksimaa, Eero and Chen, Chi‐Hua},
title = {{Cerebral Cortical Structural Variation and General Cognitive Ability: Evidence From Mendelian Randomization}},
journal = {Human brain mapping},
year = {2026},
month = sep,
volume = {47},
number = {13},
pages = {e70635},
publisher = {Wiley},
issn = {1065-9471},
doi = {10.1002/hbm.70635},
url = {https://doi.org/10.1002/hbm.70635},
pmid = {42669588},
pmcid = {PMC13526639}
}

RIS

TY - JOUR
AU - Chou, Chun‐Ju
AU - Fiecas, Mark
AU - del Re, Elisabetta C.
AU - Vuoksimaa, Eero
AU - Chen, Chi‐Hua
TI - Cerebral Cortical Structural Variation and General Cognitive Ability: Evidence From Mendelian Randomization
T2 - Human brain mapping
J2 - Hum Brain Mapp
PY - 2026
DA - 2026/09/01
VL - 47
IS - 13
SP - e70635
SN - 1065-9471
PB - Wiley
DO - 10.1002/hbm.70635
UR - https://doi.org/10.1002/hbm.70635
LA - en
ER -

CSL-JSON

{
"id": "10.1002/hbm.70635",
"type": "article-journal",
"title": "Cerebral Cortical Structural Variation and General Cognitive Ability: Evidence From Mendelian Randomization",
"container-title": "Human brain mapping",
"author": [
{
"family": "Chou",
"given": "Chun‐Ju"
},
{
"family": "Fiecas",
"given": "Mark"
},
{
"family": "del Re",
"given": "Elisabetta C."
},
{
"family": "Vuoksimaa",
"given": "Eero"
},
{
"family": "Chen",
"given": "Chi‐Hua"
}
],
"container-title-short": "Hum Brain Mapp",
"volume": "47",
"issue": "13",
"page": "e70635",
"DOI": "10.1002/hbm.70635",
"PMID": "42669588",
"PMCID": "PMC13526639",
"ISSN": "1065-9471",
"publisher": "Wiley",
"URL": "https://doi.org/10.1002/hbm.70635",
"language": "en",
"issued": {
"date-parts": [
[
2026,
9,
1
]
]
}
}

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.1038/s42003-026-10131-0 [code]
Shared genetic architecture between the topology of brain white matter structural connectome and fluid intelligence.
Journal: Communications biology
In common: structural MRI / diffusion, genetics / omics, 8 references
[2] doi:10.1371/journal.pcbi.1014422 [code]
Deciphering cell type-specific causal genetic effects on brain imaging-derived phenotypes and disorders with single-cell Mendelian randomization.
Journal: PLoS computational biology
In common: genetics / omics, 7 references
[3] doi:10.1038/s41467-026-74274-8 [code]
Regional sex differences in human cortical anatomy vary in their morphometric bases and overlap with sex chromosomal and gonadal influences.
Journal: Nature communications
In common: structural MRI / diffusion, 7 references
[4] doi:10.1038/s41467-026-73714-9 [code]
The genetic architecture of cortical similarity networks.
Journal: Nature communications
In common: structural MRI / diffusion, genetics / omics, 5 references
[5] doi:10.1038/s42003-026-09956-6 [code]
Linking changes in sulcal morphometry to cognitive development from childhood to adolescence.
Journal: Communications biology
In common: structural MRI / diffusion, 5 references
[6] doi:10.3390/life16081356
Regional Brain Volume Variation Across Adulthood: A Cross-Sectional MRI Analysis of Age, Sex, and Hemispheric Asymmetry.
Journal: Life (Basel, Switzerland)
In common: structural MRI / diffusion, 4 references
[7] doi:10.1038/s43856-026-01510-z [code]
Mapping genetic convergence across brain structure, mental health, and cardiometabolic disease.
Journal: Communications medicine
In common: genetics / omics, 4 references
[8] doi:10.1038/s41467-026-71738-9 [code]
Genetic landscape of adult executive function reveals a cell-type-specific developmental origin.
Journal: Nature communications
In common: cognitive, genetics / omics, 4 references
[9] doi:10.1162/imag.a.1235 [code]
Intracranial volume: To adjust or not to adjust? It is not a matter of if, but how.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: structural MRI / diffusion, 4 references
[10] doi:10.1038/s41380-026-03641-0
Cerebral cortical alterations in adolescent early-onset psychosis: a surface-based morphometry mega-analysis.
Journal: Molecular psychiatry
In common: structural MRI / diffusion, 4 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.

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: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

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