Multivariate genetic analysis reveals three distinct pathological dimensions in musculoskeletal disorders.

Code ↔ Paper

11 matches 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 11 matches

[1] § Methods › Multi-trait colocalization ↔ R/hyprcoloc.R, lines 911–958 · score 0.72 · HyPrColoc, posterior probability, causal variants, branch, bound, algorithm
[2] § Methods › Colocalization ↔ R/hyprcoloc.R, lines 911–958 · score 0.71 · posterior probabilities, prior probabilities, genomic region, causal variant, colocalization, signals
[3] § Methods › SMR analysis ↔ R/input_eqtl.r, lines 1–109 · score 0.70 · GTEx v8, sQTL, eQTL, tissue, GWAS, genes
[4] § Methods › Multivariate genetic modeling ↔ R/paLDSC.R, lines 1–62 · score 0.67 · sampling covariance matrix, genetic variance, components, multivariate, latent, Genomic
[5] § Methods › Multivariate genetic modeling ↔ R/commonfactor.R, lines 231–272 · score 0.66 · chi square, sampling covariance matrix, model parameters, CFI, fit, weighted
[6] § Methods › Spatial mapping of cells using gsMAP ↔ src/gsMap/report.py, lines 80–174 · score 0.65 · Cauchy combination, gsMap, gene expression, mapping, LDSC, latent
[7] § Methods › Heritability and genetic correlation ↔ R/input_eqtl.r, lines 314–396 · score 0.64 · binary phenotypes, genetic correlation matrix, sample overlap, prevalence, population, Heritability
[8] § Methods › Heritability and genetic correlation ↔ R/rgmodel.R, lines 1–35 · score 0.62 · liability scale, genetic correlation matrix, LD scores, Heritability, intercept, multivariate
[9] § Results › Three distinct genetic factors were identified in MSK phenotypes ↔ R/paLDSC.R, lines 1–62 · score 0.57 · genetic variance, genetic correlation matrix, genetic covariance, mass, dimensions, multivariate
[10] § Results › Three distinct genetic factors were identified in MSK phenotypes ↔ R/commonfactor.R, lines 314–403 · score 0.54 · model fit, correlation matrix, genetic covariance, AIC, SRMR, CFI
[11] § Methods › Spatial mapping of cells using gsMAP ↔ src/gsMap/run_all_mode.py, lines 177–241 · score 0.52 · Cauchy combination, gsMap, mapping, LDSC, correlation, gene

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 paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

R · 1,451 lines · 66 KB · no license · 2 matches

##########################################################
##### Prior #####
##########################################################
#' prior
#'
#' @param p.1 probability of one trait having a causal variant in the genetic region
#' @param gamma step probability
#' @param k number of traits
#' @export
prior <- function(p.1, gamma, k){
if(k==1){
p.1;
}else{
i=c(2:k);G=(1-gamma^(i-1));p.1*prod(G);
}
}
##########################################################
##### SNP scores #####
##########################################################
#' p.0
#'
#' @param j number of traits
#' @param Q number of SNPs
#' @export
p.0 <- function(j, Q){
(j-1)*Q - j*(j-1)/2 + 1;
}
#' ind.snp.score
#'
#' @param Q the number of traits
#' @param snp.scores vector of per snp contributions to the posterior probability of colocalization
#' @export
ind.snp.score <- function(Q, snp.scores){
loc = vector("numeric", Q-1);
res = vector("numeric", Q);
for(j in 1:Q){
i = 1;
while(i < j){
loc[i] = p.0(i, Q) + j - (i + 1);
i = i+1;
}
if(j<=Q-1){
loc[i : (Q-1)] = (p.0(j, Q)):(p.0(j+1, Q) - 1);
}
res[j] = sum(snp.scores[loc]);
}
return(res)
}
##########################################################
##### Reduce dimensions of datasets #####
##########################################################
#' samp.reduc
#'
#' @param Z matrix of Z-scores
#' @param W ratio matrix of prior standard deviation and observed standard errors
#' @param ld.matrix LD matrix
#' @param trait.cor correlation matrix between traits
#' @param sample.overlap matrix of sample overlap between traits
#' @param trts vector of traits
#' @param window.size size of window for 2CV testing
#' @param sentinel sentinel variant
#' @param Zsq matrix of Z-scores squared
#' @param Wsq matrix of W squared
#' @export
samp.reduc <- function(Z, W, ld.matrix, trait.cor, sample.overlap, trts, window.size = dim(Z)[1], sentinel = 0, Zsq, Wsq){
q = dim(Z)[1];
Z = Z[,trts];
W = W[,trts];
Zsq = Zsq[,trts];
Wsq = Wsq[,trts];
trait.cor = trait.cor[trts,trts];
sample.overlap = sample.overlap[trts,trts];
if(window.size < q){
if(sentinel == 0){
max.Z = which(abs(Z) == max(abs(Z)));
if(length(max.Z)>1){max.Z = sample(max.Z,1)};
row.Z = max.Z%%dim(Z)[1];
if(row.Z == 0){row.Z = dim(Z)[1]};
}else{ row.Z = sentinel;}
Q = window.size;
if(Q%%2 == 1){Q = Q-1;}
if(max(row.Z - Q/2,1) == 1){
row.col = 1:(Q+1);
Z = Z[row.col,];
W = W[row.col,];
Zsq = Zsq[row.col,];
Wsq = Wsq[row.col,];
ld.matrix = ld.matrix[row.col,row.col];
}else if(min(row.Z + Q/2,q)==q){
row.col = (q - Q):q;
Z = Z[row.col,];
W = W[row.col,];
Zsq = Zsq[row.col,];
Wsq = Wsq[row.col,];
ld.matrix = ld.matrix[row.col,row.col];
}else{
row.col = (row.Z-Q/2):(row.Z+Q/2);
Z = Z[row.col,];
W = W[row.col,];
Zsq = Zsq[row.col,];
Wsq = Wsq[row.col,];
ld.matrix = ld.matrix[row.col,row.col];
}
}else{row.col = 1:q}
return(list(Z, W, ld.matrix, trait.cor, sample.overlap, row.col, Zsq, Wsq))
}
#' colMax
#'
#' @param x data.frame or matrix
#' @export
colMax <- function(x){
apply(x, MARGIN = c(2), max)
}
################################################################################
##### Compute credible sets of snps for each clusetr of colocalized traits #####
################################################################################
#' cred.sets
#'
#' @param res list of results from hyprcoloc when setting "snpscores = TRUE"
#' @param value the sum of the probabilities of the snps included in the credible set
#' @export
cred.sets = function(res, value = 0.95){
clusters = which(! is.na(res[[1]]$traits));
results = vector('list', length(clusters));
if(length(clusters)>0){
count = 1;
for(i in clusters){
snp.scores = res[[2]][[i]];
snp.scores = snp.scores[order(-snp.scores)];
tmp.rmv = which(snp.scores < (1- value)/length(snp.scores));
tmp.lgt = length(snp.scores) - length(tmp.rmv);
snp.scores = snp.scores[- tmp.rmv];
if(tmp.lgt >1){
row.col = combn(tmp.lgt,2)
row = c(1:tmp.lgt, row.col[1,]);
col = c(1:tmp.lgt, row.col[2,]);
snp.cols = sparseMatrix(i = row, j = col, x = 1, dims = c(tmp.lgt,tmp.lgt));
wch.snps = snp.cols[,min(which(colSums(snp.scores*snp.cols) >= value))]==1;
results[[count]] = snp.scores[wch.snps];
}else{ results[[count]] = snp.scores[tmp.lgt];}
count = count + 1;
}
}
return(results)
}
###########################################################################################################################################
##### Perform a sensitivity analysis by varying the algorithm (regional and alignment) thresholds and coloclaization prior (prior.c) #####
###########################################################################################################################################
#' sensitivity.plot
#'
#' sensitivity.plot is a function which repeatedly calls the hyprcoloc function to compute a similarity matrix which illustrates how strongly clustered/colocalized pairs of traits are across different input thresholds and priors
#' @param effect.est matrix of snp regression coefficients (i.e. regression beta values) in the genomic region
#' @param effect.se matrix of standard errors associated with the beta values
#' @param binary.outcomes a binary vector of dimension the number of traits: 1 represents a binary trait 0 otherwise
#' @param trait.subset vector of trait names (or number) from the full trait list: used for trageted colocalization analysis in a region
#' @param trait.names vector of trait names corresponding to the columns in the effect.est matrix
#' @param snp.id vector of SNP IDs
#' @param ld.matrix LD matrix
#' @param trait.cor matrix of pairwise correlations between traits
#' @param sample.overlap matrix of pairwise sample overlap between traits
#' @param bb.alg branch and bound algorithm: TRUE, employ BB algorithm; FALSE, do not
#' @param bb.selection branch and bound algorithm type, e.g. regional or alignment selection
#' @param reg.steps regional step paramter
#' @param reg.thresh a vector of regional probability thresholds
#' @param align.thresh a vector of alignment probability thresholds
#' @param prior.1 prior probability of a SNP being associated with one trait
#' @param prior.c vector of conditional colocalization priors: where prior.c is the probability of a SNP being associated with an additional trait given that the SNP is associated with at least 1 other trait
#' @param prior.12 COLOC prior p12: prior probability of a SNP being associated with any two traits
#' @param uniform.priors uniform priors
#' @param ind.traits are the traits independent or to be treated as independent
#' @param equal.thresholds fix the regional and alignment thresholds to be equal
#' @param similarity.matrix To be viewed as a similarity matrix. Default FALSE.
#' @import pheatmap RColorBrewer
#' @export
sensitivity.plot = function(effect.est, effect.se, binary.outcomes = rep(0, dim(effect.est)[2]),
trait.subset = c(1:dim(effect.est)[2]), trait.names = c(1:dim(effect.est)[2]),
snp.id = c(1:dim(effect.est)[1]), ld.matrix = diag(1, dim(effect.est)[1], dim(effect.est)[1]),
trait.cor = diag(1, dim(effect.est)[2], dim(effect.est)[2]),
sample.overlap = matrix(rep(1,dim(effect.est)[2]^2), nrow = dim(effect.est)[2]),
bb.alg = TRUE, bb.selection = "regional",
reg.steps = 1, reg.thresh = c(0.6,0.7,0.8,0.9), align.thresh = c(0.6,0.7,0.8,0.9),
prior.1 = 1e-4, prior.c = c(0.02, 0.01, 0.005), prior.12 = NULL,
uniform.priors = FALSE,
ind.traits = TRUE, equal.thresholds = FALSE, similarity.matrix = FALSE){
m = dim(effect.est)[2];
snp.combin = function(x, y, vec){I = iterpc(x, y, labels = vec);return(getall(I)+0.0)};
sim.mat = diag(0,m);
if(!is.null(prior.12)){
prior.c = prior.12/(prior.1 + prior.12);
}
for(i in reg.thresh){
for(k in prior.c){
if(equal.thresholds){
j = i;
tmp.mat = diag(1,m);
res = hyprcoloc(effect.est, effect.se, binary.outcomes = binary.outcomes, trait.subset = trait.subset, trait.names = trait.names,
snp.id = snp.id, ld.matrix = ld.matrix, trait.cor = trait.cor, sample.overlap = sample.overlap, bb.alg = bb.alg, bb.selection = bb.selection,
reg.steps = reg.steps, reg.thresh = i, align.thresh = j,
prior.1 = prior.1, prior.c = k, uniform.priors = uniform.priors, ind.traits = ind.traits);
trt.clusts = res[[1]]$traits;
for(its in 1:length(trt.clusts)){
tmp.clust = unlist(strsplit(trt.clusts[its], split=", "));
if(tmp.clust[1]!="None"){
tmp.vec = which(trait.names %in% tmp.clust);
# coloc.pairs = t(snp.combin(m, 2, tmp.vec));
coloc.pairs = t(snp.combin(length(tmp.vec), 2, tmp.vec));
tmp.mat[t(coloc.pairs)] = 1;
}
}
sim.mat = sim.mat + tmp.mat + t(tmp.mat) - diag(1,m);
}else{
for(j in align.thresh){
tmp.mat = diag(1,m);
res = hyprcoloc(effect.est, effect.se, binary.outcomes = binary.outcomes, trait.subset = trait.subset, trait.names = trait.names,
snp.id = snp.id, ld.matrix = ld.matrix, trait.cor = trait.cor, sample.overlap = sample.overlap, bb.alg = bb.alg, bb.selection = bb.selection,
reg.steps = reg.steps, reg.thresh = i, align.thresh = j,
prior.1 = prior.1, prior.c = k, uniform.priors = uniform.priors, ind.traits = ind.traits);
trt.clusts = res[[1]]$traits;
for(its in 1:length(trt.clusts)){
tmp.clust = unlist(strsplit(trt.clusts[its], split=", "));
if(tmp.clust[1]!="None"){
tmp.vec = which(trait.names %in% tmp.clust);
# coloc.pairs = t(snp.combin(m, 2, tmp.vec));
coloc.pairs = t(snp.combin(length(tmp.vec), 2, tmp.vec));
tmp.mat[t(coloc.pairs)] = 1;
}
}
sim.mat = sim.mat + tmp.mat + t(tmp.mat) - diag(1,m);
}
}
}
}
sim.mat = sim.mat/length(reg.thresh)/length(align.thresh)/length(prior.c);
if(equal.thresholds){
sim.mat = sim.mat*length(align.thresh);
}
rownames(sim.mat) = trait.names;
colnames(sim.mat) = trait.names;
# annotation_row = data.frame(
# Clusters = factor(dta$cluster.class[match(rownames(smat), obs.names)])
# )
# rownames(annotation_row) <- rownames(smat)
breaksList = seq(0,1,by=0.02);
plot = pheatmap(
mat = sim.mat,
color = colorRampPalette((brewer.pal(n = 9, name = "Blues")))(length(breaksList)),
breaks = breaksList,
cluster_cols = FALSE,
cluster_rows = FALSE,
border_color = NA,
show_colnames = TRUE,
show_rownames = TRUE,
#annotation_col = annotation_row,
drop_levels = TRUE,
fontsize = 6,
main = "Prior sensitivity heatmap"
)
if(!similarity.matrix){
return(plot)}else{
return(list(plot, sim.mat))
}
}
##########################################################
##### Regional colocalisation (rapid) #####
##########################################################
#' rapid.reg
#'
#' @param Zsq matrix of Z-scores
#' @param Wsq ratio matrix of prior standard deviation and observed standard errors squared
#' @param prior.1 prior probability of a SNP being associated with one trait
#' @param prior.2 1 - prior probability of a SNP being associated with an additional trait given that the SNP is associated with at least 1 other trait
#' @param uniform.priors uniform priors
#' @export
rapid.reg <- function(Zsq, Wsq, prior.1, prior.2, uniform.priors){
m = dim(Zsq)[2];
Q = dim(Zsq)[1];
if(uniform.priors==T){
I.unif = 1;
}else{
I.unif = 0;
}
p.1 = prior.1;
p.1.m = p.1/Q;
prior.all = I.unif*p.1.m + (1-I.unif)*prior(prior.1, prior.2, k = m);
prior.sub = I.unif*p.1.m + (1-I.unif)*prior(prior.1, prior.2, k = m-1);
labf = 0.5*(log(Wsq) + (Zsq)*(1- Wsq));
sum.labf = rowSums(labf);
mx.labf = which.max(labf);
clc.max = which.max(sum.labf);
chi = sum.labf - sum.labf[clc.max];
exp.chi = exp(chi);
sum.exp.chi = sum(exp.chi);
sum.labf.subset = colSums(exp(chi - labf));
rm.trt = which.max(sum.labf.subset);
reg.prob = sum.exp.chi/(exp(-sum.labf[clc.max] - log(prior.all)) + sum.exp.chi + (prior.sub/prior.all)*sum(sum.labf.subset));
return(list(reg.prob, clc.max, rm.trt, exp.chi/sum.exp.chi))
}
##########################################################
##### Regional colocalisation #####
##########################################################
#' regional.ABF
#'
#' @param Z matrix of Z-scores
#' @param W ratio matrix of prior standard deviation and observed standard errors
#' @param snps.clc SNPs colocalisation
#' @param rho LD matrix
#' @param trait.cor correlation matrix between traits
#' @param sample.overlap matrix of sample overlap between traits
#' @param epsilon tolerance parameter
#' @param reg.thresh regional probability threshold
#' @param prior.1 prior probability of a SNP being associated with one trait
#' @param prior.2 1 - prior probability of a SNP being associated with an additional trait given that the SNP is associated with at least 1 other trait
#' @param prior.3 prior probability that a trait contains a second causal variant given it contains one already
#' @param prior.4 1 - prior probability that trait two co-localises with trait one given traits one and two already share a causal variant and trait one contains a second causal variant
#' @param flag flag variable
#' @param test.2 test for 2CV
#' @param reg.steps regional step paramter
#' @param cor.adj.priors correlation adjusted priors
#' @param uniform.priors uniform priors
#' @param branch.jump branch jump
#' @param Zsq matrix of Z-scores squared
#' @param Wsq matrix of W squared
#' @param ind.traits are the traits independent or to be treated as independent
#' @export
regional.ABF <- function(Z, W, snps.clc, rho, trait.cor, sample.overlap, epsilon, reg.thresh, prior.1, prior.2, prior.3, prior.4, flag, test.2, reg.steps, cor.adj.priors, uniform.priors, branch.jump, Zsq, Wsq, ind.traits){
m = dim(Z)[2];
Q = dim(Z)[1];
trait.cor =trait.cor*sample.overlap;
if(reg.steps>m){reg.steps = m;}
if(uniform.priors == T){
I.unif = 1;
}else{
I.unif = 0;
}
p.1 = prior.1;
p.1.m = p.1/Q;
p.2.m = p.1*2/(Q*(Q-1));
if(cor.adj.priors==T){
ave.cor = trait.cor[lower.tri(trait.cor)];
ave.cor = mean(abs(ave.cor));
prior.2 = min(prior.2, 1-ave.cor^2);
prior.4 = min(prior.4, 1-ave.cor^4);
}
if(branch.jump == T & reg.steps > 2){
jmp = 1;
}else{
jmp = 0;
}
ones=matrix(1,nrow=dim(snps.clc)[1],ncol=dim(snps.clc)[1]);
trait.cor = as.matrix(kronecker(trait.cor, ones));
log.sum.max.ABF.1 = sum.max.ABF.1 = log.sum.max.ABF.2 = sum.max.ABF.2 = mpfr(0,120);
log.sum.ABF.1 = sum.ABF.1 = sum.ABF.all.1 = log.sum.ABF.2 = sum.ABF.2 = sum.ABF.all.2 = mpfr(0,120);
log.all.traits.ABF.1 = log.all.traits.ABF.2 = mpfr(0,120);
log.max.ABF.1 = log.max.ABF.2 = mpfr(0,120)
mx.ABF.k = mpfr(0,120);
NAs = 0;
if(test.2 ==0){
df = data.frame(matrix(vector(), 0,10, dimnames=list(c(), c("log.sum.ABF.full", "log.sum.ABF.all.combs", "log.max.sum.ABF.all.traits", "log.max.SNP.ABF.full", "Post.ratio.2.vs.1and2", "Traits", "SNPs", "NaNs", "reg_bb_alg" , "reg_only_trts" ))), stringsAsFactors=F);
df[1,]$reg_bb_alg = FALSE;
reg.prob = 1;
i = m;
count = 0;
while(reg.prob > (1-jmp)*reg.thresh & i > 0){
clc.trt=combn(m,i);
count = count + dim(clc.trt)[2];
prior1 = I.unif*p.1.m + (1-I.unif)*prior(prior.1,prior.2, k = i);
prior2 = I.unif*p.2.m + (1-I.unif)*prior(prior.1,prior.2, k = i)*prior(prior.3, prior.4, k = i);
for(j in 1:dim(clc.trt)[2]){
trt.clc = clc.trt[,j] + 0.0;
if(flag==0){
if(ind.traits == T){
output = regional1ind(Zsq, Wsq, trt.clc);
}else{
output = regional1(Z, W, trt.clc, trait.cor, epsilon);
}
if(i == m){
max.ABF.1 = prior1*exp(mpfr(output[[1]][1], 120));
sum.ABF.1 = max.ABF.1*output[[1]][2];
clc.max.cvs.1 = output[[1]][3];
}else{
max.ABF.1 = prior1*exp(mpfr(output[[1]][1], 120));
sum.ABF.1 = max.ABF.1*output[[2]][1];
clc.max.cvs.1 = output[[3]][1];
if(mx.ABF.k < sum.ABF.1){
mx.ABF.k = sum.ABF.1;
df[1,]$reg_only_trts = toString(clc.trt[,j]);
}
}
if(!is.na(sum.ABF.1)){
sum.ABF.all.1 = sum.ABF.all.1 + sum.ABF.1;
if(sum.max.ABF.1<sum.ABF.1){
sum.max.ABF.1 = sum.ABF.1;
df[1,]$Traits = toString(clc.trt[,j]);
df[1,]$SNPs = toString(clc.max.cvs.1)
}
}else{NAs = NAs + 1;}
}else{
output = regional2(Z, W, snps.clc, trt.clc, rho, trait.cor, epsilon);
max.ABF.1 = prior1*exp(mpfr(output[[1]][1], 120));
sum.ABF.1 = max.ABF.1*output[[1]][2];
clc.max.cvs.1 = output[[1]][3];
if(!is.na(sum.ABF.1)){
sum.ABF.all.1 = sum.ABF.all.1 + sum.ABF.1;
if(sum.max.ABF.1<sum.ABF.1){
sum.max.ABF.1 = sum.ABF.1;
df[1,]$Traits = toString(clc.trt[,j]);
df[1,]$SNPs = toString(clc.max.cvs.1)
}
}else{NAs = NAs + 1;}
max.ABF.2 = prior2*exp(mpfr(output[[2]][1], 120));
sum.ABF.2 = max.ABF.2*output[[2]][2];
clc.max.cvs.2 = snps.clc[,output[[2]][3]];
if(!is.na(sum.ABF.2)){
sum.ABF.all.2 = sum.ABF.all.2 + sum.ABF.2;
if(sum.max.ABF.2<sum.ABF.2){
sum.max.ABF.2 = sum.ABF.2;
df[2,]$Traits = toString(clc.trt[,j]);
df[2,]$SNPs = toString(clc.max.cvs.2)
}
}else{NAs = NAs + 1;}
}
}
if(i==m){
log.max.ABF.1 = log(max.ABF.1)
sum.ABF.full.1 = sum.ABF.all.1;
log.sum.ABF.1=asNumeric(log(sum.ABF.1));
df[1,]$log.sum.ABF.full=log.sum.ABF.1;
df[1,]$log.max.SNP.ABF.full=asNumeric(log.max.ABF.1);
snp.scores.1 = output[[2]];
if(flag==1){
log.max.ABF.2 = log(max.ABF.2);
sum.ABF.full.2 = sum.ABF.all.2;
log.sum.ABF.2=asNumeric(log(sum.ABF.2));
df[2,]$log.sum.ABF.full=log.sum.ABF.2;
df[2,]$log.max.SNP.ABF.full=asNumeric(log.max.ABF.2);
snp.scores.1 = output[[3]];
snp.scores.2.tmp = output[[4]];
snp.scores.2 = ind.snp.score(Q, snp.scores.2.tmp);
}
}
if(flag==0){
reg.prob = sum.ABF.full.1/(sum.ABF.all.1);
}else{
reg.prob = (sum.ABF.full.1+sum.ABF.full.2)/(sum.ABF.all.1+sum.ABF.all.2);
}
i = i - 1;
if(reg.steps!=0 & i == (m-reg.steps-1)){i=0;}
}
if(jmp == 1 | sum.ABF.full.1/mx.ABF.k >= (1/min(9,m))*count){df[1,]$reg_bb_alg = TRUE;}
log.sum.max.ABF.1=asNumeric(log(sum.max.ABF.1));
log.all.traits.ABF.1=asNumeric(log(sum.ABF.all.1));
df[1,]$log.sum.ABF.all.combs=log.all.traits.ABF.1;
df[1,]$log.max.sum.ABF.all.traits=log.sum.max.ABF.1;
df[1,]$NaNs=NAs;
out.list = list(df, snp.scores.1);
if(flag==1){
log.sum.max.ABF.2=asNumeric(log(sum.max.ABF.2));
log.all.traits.ABF.2=asNumeric(log(sum.ABF.all.2));
df[2,]$log.sum.ABF.all.combs=log.all.traits.ABF.2;
df[2,]$log.max.sum.ABF.all.traits=log.sum.max.ABF.2;
df[2,]$NaNs=NAs;
df[2,]$Post.ratio.2.vs.1and2 = asNumeric(sum.ABF.all.2/(sum.ABF.all.1 + sum.ABF.all.2));
out.list = list(df, snp.scores.1, snp.scores.2);
}
}else{
prior1 = I.unif*p.1.m + (1-I.unif)*prior(prior.1,prior.2, k = m);
prior2 = I.unif*p.2.m + (1-I.unif)*prior(prior.1,prior.2, k = m)*prior(prior.3, prior.4, k = m);
output = regional2(Z, W, snps.clc, c(1:m)+0.0, rho, trait.cor, epsilon);
max.ABF.1 = prior1*exp(mpfr(output[[1]][1], 120));
sum.ABF.1 = max.ABF.1*output[[1]][2];
max.ABF.2 = prior2*exp(mpfr(output[[2]][1], 120));
sum.ABF.2 = max.ABF.2*output[[2]][2];
out.list = list(sum.ABF.2/(sum.ABF.2 + sum.ABF.1));
}
return(out.list)
}
##########################################################
##### Alignment (rapid) #####
##########################################################
#' rapid.align
#'
#' @param Zsq matrix of Z-scores
#' @param Wsq ratio matrix of prior standard deviation and observed standard errors squared
#' @param prior.1 prior probability of a SNP being associated with one trait
#' @param prior.2 1 - prior probability of a SNP being associated with an additional trait given that the SNP is associated with at least 1 other trait
#' @param uniform.priors uniform priors
#' @export
rapid.align <- function(Zsq, Wsq, prior.1, prior.2, uniform.priors){
m = dim(Zsq)[2];
Q = dim(Zsq)[1];
if(uniform.priors==T){
I.unif = 1;
}else{
I.unif = 0;
}
p.1 = prior.1;
p.1.m = p.1/Q;
p.m.1 = p.1/(m*Q*(Q-1));
if(m==2){
p.m.1 = 2*p.m.1;
cnst = 0.25;
}else{
cnst = 1
}
prior.all = I.unif*p.1.m + (1-I.unif)*prior(prior.1, prior.2, k = m);
prior.align= I.unif*p.m.1 + (1-I.unif)*prior.1*prior(prior.1, prior.2, k = m-1);
labf = 0.5*(log(Wsq) + (Zsq)*(1- Wsq));
sum.labf = rowSums(labf);
mx.labf = which.max(labf);
clc.max = which.max(sum.labf);
chi = sum.labf - sum.labf[clc.max];
sum.exp.chi = sum(exp(chi));
col.max.labf = colMax(labf);
exp.labf.std = exp(t(t(labf) - col.max.labf));
labf.tmp = t(colSums(exp.labf.std) - t(exp.labf.std));
sum.tmp = colSums(exp(t(t(chi - labf)+col.max.labf))*labf.tmp);
rm.trt = which.max(sum.tmp);
align.prob = sum.exp.chi/(sum.exp.chi + cnst*(prior.align/prior.all)*sum(sum.tmp));
return(list(align.prob, rm.trt))
}
##########################################################
##### Alignment #####
##########################################################
#' align.ABF.1
#'
#' @param Z matrix of Z-scores
#' @param W ratio matrix of prior standard deviation and observed standard errors
#' @param trait.cor correlation matrix between traits
#' @param sample.overlap matrix of sample overlap between traits
#' @param ld.matrix LD matrix
#' @param epsilon tolerance parameter
#' @param reg.res regional result
#' @param align.thresh alignment probability threshold
#' @param prior.1 prior probability of a SNP being associated with one trait
#' @param prior.2 1 - prior probability of a SNP being associated with an additional trait given that the SNP is associated with at least 1 other trait
#' @param cor.adj.priors correlation adjusted priors
#' @param uniform.priors uniform priors
#' @param Zsq matrix of Z-scores squared
#' @param Wsq matrix of W squared
#' @param ind.traits are the traits independent or to be treated as independent
#' @export
align.ABF.1 <- function(Z, W, trait.cor, sample.overlap, ld.matrix, epsilon, reg.res, align.thresh, prior.1, prior.2, cor.adj.priors, uniform.priors, Zsq, Wsq, ind.traits){
m = dim(Z)[2];
Q = dim(Z)[1];
traits = c(1:m);
trait.cor =trait.cor*sample.overlap;
if(uniform.priors == T){
I.unif = 1;
}else{
I.unif = 0;
}
p.1 = prior.1;
p.m.1 = p.1/(m*Q*(Q-1));
if(m==2){
p.m.1 = 2*p.m.1;
}
prior.3 = prior.1;
kappa = 2;
if(cor.adj.priors==T){
ave.cor = trait.cor[lower.tri(trait.cor)];
ave.cor = mean(abs(ave.cor));
prior.2 = min(prior.2, 1-ave.cor^2);
prior.3 = prior.1*(10^(-kappa*ave.cor^2));
p.m.1 = p.m.1*(1-ave.cor)^kappa;
}
max.ABF.overall = sum.max.ABF.co = mpfr(0,120);
log.sum.ABF.all.combs = sum.ABF = sum.ABF.co = mpfr(0,120);
NaNs = 0;
df = data.frame(matrix(vector(), 0,3, dimnames=list(c(), c("log.sum.ABF.all.combs", "Trait.no.clc", "NaNs"))), stringsAsFactors=F)
j = 1;
align.prob = 1;
trt.combn = sample(c(1:m));
prior.algn= I.unif*p.m.1 + (1-I.unif)*prior.3*prior(prior.1, prior.2, k = m-1);
if(m>2){
while(j <=m){
trt.no.clc = trt.combn[j] + 0.0;
trt.clc = traits[-trt.no.clc] + 0.0;
if(ind.traits==T){
output = align1ind(Zsq, Wsq, trt.clc, trt.no.clc);
}else{
output = align1(Z, W, 1, trt.clc, trt.no.clc, trait.cor, ld.matrix, epsilon);
}
max.ABF = prior.algn*exp(mpfr(output[[1]][1], 120));
sum.ABF.tmp = max.ABF*output[[2]][1];
if(!is.nan(sum.ABF.tmp)){
sum.ABF= sum.ABF + sum.ABF.tmp;
if(max.ABF.overall<sum.ABF.tmp){
max.ABF.overall = sum.ABF.tmp;
df[1,]$Trait.no.clc = toString(trt.no.clc);
}
}else{NaNs = NaNs + 1}
align.prob = (reg.res)/(reg.res + sum.ABF);
j = j+1;
}
}else{
trt.no.clc = 1.0;
trt.clc = traits[-trt.no.clc] + 0.0;
if(ind.traits == TRUE){
output = align1ind(Zsq, Wsq, trt.clc, trt.no.clc);
}else{
output = align1(Z, W, 1, trt.clc, trt.no.clc, trait.cor, ld.matrix, epsilon);
}
max.ABF = prior.algn*exp(mpfr(output[[1]][1], 120));
sum.ABF.tmp = max.ABF*output[[2]][1];
if(!is.nan(sum.ABF.tmp)){
sum.ABF= (sum.ABF + sum.ABF.tmp)/2;
df[1,]$Trait.no.clc = toString(trt.no.clc);
}else{NaNs = NaNs + 1}
}
log.sum.ABF.all.combs=asNumeric(log(sum.ABF));
df[1,]$log.sum.ABF.all.combs=log.sum.ABF.all.combs;
df[1,]$NaNs = NaNs;
return(df)
}
##########################################################
##### Alignment 1CV and 2CV #####
##########################################################
#' align.ABF.2
#'
#' @param Z matrix of Z-scores
#' @param W ratio matrix of prior standard deviation and observed standard errors
#' @param snps.clc SNPs colocalisation
#' @param trait.cor correlation matrix between traits
#' @param sample.overlap matrix of sample overlap between traits
#' @param ld.matrix LD matrix
#' @param epsilon tolerance parameter
#' @param reg.res regional result
#' @param align.thresh alignment probability threshold
#' @param prior.1 prior probability of a SNP being associated with one trait
#' @param prior.2 1 - prior probability of a SNP being associated with an additional trait given that the SNP is associated with at least 1 other trait
#' @param prior.3 prior probability that a trait contains a second causal variant given it contains one already
#' @param prior.4 1 - prior probability that trait two co-localises with trait one given traits one and two already share a causal variant and trait one contains a second causal variant
#' @param cor.adj.priors correlation adjusted priors
#' @param uniform.priors uniform priors
#' @export
align.ABF.2 <- function(Z, W, snps.clc, trait.cor, sample.overlap, ld.matrix, epsilon, reg.res, align.thresh, prior.1, prior.2, prior.3, prior.4, cor.adj.priors, uniform.priors){
m = dim(Z)[2];
Q = dim(Z)[1];
traits = c(1:m);
n.cv = dim(snps.clc)[1];
trait.cor =trait.cor*sample.overlap;
if(uniform.priors == T){
I.unif = 1;
}else{
I.unif = 0;
}
p.1 = prior.1;
p.m.1 = p.1/(m*Q*(Q-1));
p.1.m.2 = 2*p.1/(m*Q*(Q-1)*(Q-2));
p.2.m.1 = 2*p.1/(m*Q*(Q-1)*(Q-2));
p.2.m.2 = 4*p.1/(m*Q*(Q-1)*(Q-2)*(Q-3));
p.co.m.1.2 = p.1/(m*Q*(Q-1));
p.co.m.2.1 = p.1/(m*Q*(Q-1));
p.co.m.2.2 = p.1/(m*Q*(Q-1)*(Q-2));
if(m==2){
p.m.1 = 2*p.m.1;
p.2.m.2 = 2*p.2.m.2
p.co.m.1.2 = p.co.m.1.2/2;
p.co.m.2.1 = p.co.m.2.1/2;
p.co.m.2.2 = 2*p.co.m.2.2
}
prior.5 = prior.1
kappa = 10;
if(cor.adj.priors==T){
ave.cor = trait.cor[lower.tri(trait.cor)];
ave.cor = mean(abs(ave.cor));
prior.2 = min(prior.2, 1-ave.cor^2);
prior.3 = prior.3*(10^(-kappa*ave.cor^2));
prior.4 = min(prior.4, 1-ave.cor^4);
prior.5 = prior.1*(10^(-kappa*ave.cor^2));
p.m.1 = p.m.1*(1-ave.cor)^2;
p.1.m.2 = p.1.m.2*(1-ave.cor)^4;
p.2.m.1 = p.2.m.1/(1-ave.cor)^2;
p.2.m.2 = p.2.m.2*(1-ave.cor)^4;
p.co.m.1.2 = p.co.m.1.2*(1-ave.cor)^2
p.co.m.2.1 = p.co.m.2.1*(1-ave.cor)^2
p.co.m.2.2 = p.co.m.2.2*(1-ave.cor)^2
}
ones=matrix(1,nrow=dim(snps.clc)[1], ncol=dim(snps.clc)[1]);
trait.cor = as.matrix(kronecker(trait.cor, ones));
max.ABF.overall = sum.max.ABF.co = mpfr(0,120);
log.sum.ABF.align = sum.ABF = sum.ABF.co = mpfr(0,120);
NAs = 0;
df = data.frame(matrix(vector(), 0,7, dimnames=list(c(), c("log.sum.ABF.co.loc", "log.sum.ABF.align", "log.max.SNP.ABF.full", "Trait.no.clc", "Traits.clc", "SNPs.clc", "NaNs"))), stringsAsFactors=F)
j = 1;
align.prob = 1;
trt.combn = sample(c(1:m));
prior.no.11 = I.unif*p.m.1 + (1-I.unif)*prior.5*prior( prior.1, prior.2, k = m-1);
prior.no.12 = I.unif*p.1.m.2 + (1-I.unif)*(prior.5*prior.3)*prior(prior.1, prior.2, k = m-1);
prior.co.12 = I.unif*p.co.m.1.2 + (1-I.unif)*prior.3*prior(prior.1, prior.2, k = m);
prior.no.21 = I.unif*p.2.m.1 + (1-I.unif)*prior.5*prior( prior.1, prior.2, k = m-1)*prior(prior.3, prior.4, k = m-1);
prior.no.22 = I.unif*p.2.m.2 + (1-I.unif)*(prior.5*prior.3)*prior( prior.1, prior.2, k = m-1)*prior(prior.3, prior.4, k = m-1);
prior.co.21 = I.unif*p.co.m.2.1 + (1-I.unif)*prior(prior.1, prior.2, k = m)*prior(prior.3, prior.4, k = m-1);
prior.co.22 = I.unif*p.co.m.2.2 + (1-I.unif)*prior.3*prior(prior.1, prior.2, k = m)*prior(prior.3, prior.4, k = m-1);
while(align.prob > align.thresh & j <=m){
trt.no.clc = trt.combn[j] + 0.0;
trt.clc = traits[-trt.no.clc] + 0.0;
output_1 = align12(Z, W, snps.clc, trt.clc, trt.no.clc, ld.matrix, trait.cor, epsilon);
max.ABF = c(prior.no.11*exp(mpfr(output_1[[1]][1], 120)),prior.no.12*exp(mpfr(output_1[[2]][1], 120)));
sum.ABF.no = sum(max.ABF*c(output_1[[1]][2], output_1[[2]][2]));
if(!is.nan(sum.ABF.no)){
sum.ABF= sum.ABF + sum.ABF.no;
if(max.ABF.overall<sum.ABF.no){
max.ABF.overall = sum.ABF.no;
df[1,]$Trait.no.clc = toString(trt.no.clc);
}
}else{NAs = NAs + 1}
max.ABF.1 = prior.co.12*exp(mpfr(output_1[[3]][1], 120));
sum.ABF.1 = max.ABF.1*output_1[[3]][2];
clc.max.cvs.1 = output_1[[3]][3];
if(!is.na(sum.ABF.1)){
sum.ABF.co = sum.ABF.co + sum.ABF.1;
if(sum.max.ABF.co<sum.ABF.1){
sum.max.ABF.co = sum.ABF.1;
log.max.ABF.co = log(max(max.ABF.1))
df[1,]$Traits.clc = toString(trt.clc);
df[1,]$SNPs.clc = toString(clc.max.cvs.1)
}
}else{NAs = NAs + 1;}
output_2 = align2(Z, W, snps.clc, trt.clc, trt.no.clc, ld.matrix, trait.cor, epsilon);
max.ABF = c(prior.no.21*exp(mpfr(output_2[[1]][1], 120)),prior.no.22*exp(mpfr(output_2[[2]][1], 120)));
sum.ABF.no = sum(max.ABF*c(output_2[[1]][2], output_2[[2]][2]));
if(!is.nan(sum.ABF.no)){
sum.ABF= sum.ABF + sum.ABF.no;
if(max.ABF.overall<sum.ABF.no){
max.ABF.overall = sum.ABF.no;
df[1,]$Trait.no.clc = toString(trt.no.clc);
}
}else{NAs = NAs + 1}
max.ABF.2 = c(prior.co.21*exp(mpfr(output_2[[3]][1], 120)), prior.co.22*exp(mpfr(output_2[[4]][1], 120)));
sum.ABF.2 = sum(max.ABF.2*c(output_2[[3]][2], output_2[[4]][2]));
clc.max.cvs.2 = snps.clc[output_2[[2+which(max.ABF.2==max(max.ABF.2))]][3],output_2[[2+which(max.ABF.2==max(max.ABF.2))]][4]];
if(!is.na(sum.ABF.2)){
sum.ABF.co = sum.ABF.co + sum.ABF.2;
if(sum.max.ABF.co<sum.ABF.2){
sum.max.ABF.co = sum.ABF.2;
log.max.ABF.co = log(max(max.ABF.2));
df[1,]$Traits.clc = toString(trt.clc);
df[1,]$SNPs.clc = toString(clc.max.cvs.2)
}
}else{NAs = NAs + 1;}
align.prob = (reg.res + sum.ABF.co)/(reg.res + sum.ABF.co + sum.ABF);
j = j+1;
}
log.sum.ABF.align=asNumeric(log(sum.ABF));
df[1,]$log.sum.ABF.align=log.sum.ABF.align;
df[1,]$NaNs=NAs;
log.sum.ABF.co.loc = asNumeric(log(sum.ABF.co));
df[1,]$log.sum.ABF.co.loc=log.sum.ABF.co.loc;
df[1,]$log.max.SNP.ABF.full=asNumeric(log.max.ABF.co);
return(df)
}
##########################################################
##### HyPrColoc (rapid) #####
##########################################################
#' rapid.hyprcoloc
#'
#' @param Zsq matrix of Z-scores
#' @param Wsq ratio matrix of prior standard deviation and observed standard errors squared
#' @param prior.1 prior probability of a SNP being associated with one trait
#' @param prior.2 1 - prior probability of a SNP being associated with an additional trait given that the SNP is associated with at least 1 other trait
#' @param uniform.priors uniform priors
#' @export
rapid.hyprcoloc <- function(Zsq, Wsq, prior.1, prior.2, uniform.priors){
m = dim(Zsq)[2];
Q = dim(Zsq)[1];
if(uniform.priors==T){
I.unif = 1;
}else{
I.unif = 0;
}
p.1 = prior.1;
p.1.m = p.1/Q;
p.m.1 = p.1/(m*Q*(Q-1));
if(m==2){
p.m.1 = 2*p.m.1;
cnst = 0.25;
}else{
cnst = 1
}
prior.all = I.unif*p.1.m + (1-I.unif)*prior(prior.1, prior.2, k = m);
prior.sub = I.unif*p.1.m + (1-I.unif)*prior(prior.1, prior.2, k = m-1);
prior.align= I.unif*p.m.1 + (1-I.unif)*prior.1*prior(prior.1, prior.2, k = m-1);
labf = 0.5*(log(Wsq) + (Zsq)*(1- Wsq));
sum.labf = rowSums(labf);
mx.labf = which.max(labf);
clc.max = which.max(sum.labf);
chi = sum.labf - sum.labf[clc.max];
exp.chi = exp(chi);
sum.exp.chi = sum(exp.chi);
sum.labf.subset = colSums(exp(chi - labf));
rm.trt = which.max(sum.labf.subset);
reg.prob = sum.exp.chi/(exp(-sum.labf[clc.max] - log(prior.all)) + sum.exp.chi + (prior.sub/prior.all)*sum(sum.labf.subset));
col.max.labf = colMax(labf);
exp.labf.std = exp(t(t(labf) - col.max.labf));
labf.tmp = t(colSums(exp.labf.std) - t(exp.labf.std));
sum.tmp = colSums(exp(t(t(chi - labf)+col.max.labf))*labf.tmp);
align.prob = sum.exp.chi/(sum.exp.chi + cnst*(prior.align/prior.all)*sum(sum.tmp));
return(list(reg.prob, reg.prob*align.prob, clc.max, exp.chi/sum.exp.chi))
}
##########################################################
##### HyPrColoc #####
##########################################################
#' HyPrColoc
#'
#' hyprcoloc is a function used to identify clusters of colocalized traits and candidate causal SNPs in genomic regions
#' @param effect.est matrix of snp regression coefficients (i.e. regression beta values) in the genomic region
#' @param effect.se matrix of standard errors associated with the beta values
#' @param binary.outcomes a binary vector of dimension the number of traits: 1 represents a binary trait 0 otherwise
#' @param trait.subset vector of trait names (or number) from the full trait list: used for trageted colocalization analysis in a region
#' @param trait.names vector of trait names corresponding to the columns in the effect.est matrix
#' @param snp.id vector of SNP IDs
#' @param ld.matrix LD matrix
#' @param trait.cor matrix of pairwise correlations between traits
#' @param sample.overlap matrix of pairwise sample overlap between traits
#' @param bb.alg branch and bound algorithm: TRUE, employ BB algorithm; FALSE, do not
#' @param bb.selection branch and bound algorithm type, e.g. regional or alignment selection
#' @param reg.steps regional step paramter
#' @param reg.thresh threshold probability beyond which traits are believed to share a regional association signal
#' @param align.thresh threshold probability beyond which traits are believed to align at a single causal variant
#' @param prior.1 prior probability of a SNP being associated with one trait
#' @param prior.c conditional colocalization prior: probability of a SNP being associated with an additional trait given that the SNP is associated with at least 1 other trait
#' @param prior.12 COLOC prior p12: prior probability of a SNP being associated with any two traits
#' @param sensitivity perform senstivity analysis
#' @param sense.1 first sensitivity analysis
#' @param sense.2 second sensitivity analysis
#' @param uniform.priors uniform priors
#' @param ind.traits are the traits independent or to be treated as independent
#' @param snpscores output estimated posterior probability explained each SNP
#' @return A data.frame of HyPrColoc results: each row is a cluster of colocalized traits or is coded NA (if no colocalization is identified)
#' @return If snpscores=TRUE: additionally returns a list of posterior probability explained by each SNPs and for each cluster of colocalized traits identified
#' @import compiler Rmpfr iterpc Matrix
#' @importFrom Rcpp evalCpp
#' @useDynLib hyprcoloc
#' @author Christopher N Foley <[email hidden]> and James R Staley <[email hidden]>
#' @examples
#' # Regression coefficients and standard errors from ten GWAS studies
#' # (Traits 1-5, 6-8 & 9-10 are the clusters of colocalized traits)
#' betas <- hyprcoloc::test.betas
#' head(betas)
#' ses <- hyprcoloc::test.ses
#' head(ses)
#'
#' # Trait names and SNP IDs
#' traits <- paste0("T", 1:10)
#' rsid <- rownames(betas)
#'
#' # Colocalisation analyses
#' results <- hyprcoloc(betas, ses, trait.names=traits, snp.id=rsid)
#' @export
hyprcoloc <- function(effect.est, effect.se, binary.outcomes = rep(0, dim(effect.est)[2]),
trait.subset = c(1:dim(effect.est)[2]), trait.names = c(1:dim(effect.est)[2]),
snp.id = c(1:dim(effect.est)[1]), ld.matrix = diag(1, dim(effect.est)[1], dim(effect.est)[1]),
trait.cor = diag(1, dim(effect.est)[2], dim(effect.est)[2]),
sample.overlap = matrix(rep(1,dim(effect.est)[2]^2) , nrow = dim(effect.est)[2]),
bb.alg = TRUE, bb.selection = "regional",
reg.steps = 1, reg.thresh = "default", align.thresh = "default",
prior.1 = 1e-4, prior.c = 0.02, prior.12 = NULL,
sensitivity = FALSE, sense.1 = 1, sense.2 = 2,
uniform.priors = FALSE, ind.traits = FALSE, snpscores=FALSE){
if(any(is.na(effect.est))) stop("there are missing values in effect.est")
if(any(is.na(effect.se))) stop("there are missing values in effect.se")
if(any(effect.se==0)) stop("there are zero values in effect.se")
if(any(is.na(binary.outcomes))) stop("there are missing values in binary.outcomes")
if(any(!(binary.outcomes %in% c(0,1)))) stop("there are missing values in binary.outcomes")
if(any(is.na(trait.subset))) stop("there are missing values in trait.subset")
if(any(is.na(trait.names))) stop("there are missing values in trait.names")
if(any(is.na(snp.id))) stop("there are missing values in snp.id")
if(any(is.na(ld.matrix))) stop("there are missing values in ld.matrix")
if(any(is.na(trait.cor))) stop("there are missing values in trait.cor")
if(any(is.na(sample.overlap))) stop("there are missing values in sample.overlap")
if(!is.null(prior.12)){
prior.c = prior.12/(prior.1 + prior.12);
}
prior.2 = 1-prior.c;
n.cvs = 1;
test.2 = F;
sentinel = 0;
epsilon = 0;
window.size = dim(effect.est)[1];
prior.3 = 1e-3;
prior.4 = 0.995;
reg.tol = 0.699;
cor.adj.priors = F;
branch.jump = F;
Z = effect.est/effect.se;
m = dim(Z)[2];
Q = dim(Z)[1];
if(uniform.priors==F & (reg.thresh == "default" | align.thresh == "default")){
if(reg.thresh == "default"){reg.thresh = 0.5; reg.tol = 0.499;}
if(align.thresh == "default"){align.thresh = 0.5;}
}else if(uniform.priors==T & (reg.thresh == "default" | align.thresh == "default")){
if(reg.thresh == "default"){reg.thresh =0.7;}
if(align.thresh == "default"){align.thresh = 0.7;}
}
if(reg.tol >= reg.thresh){
reg.tol = reg.thresh-0.001;
}
if(bb.alg == T & uniform.priors==F & (reg.thresh < 0.5 | align.thresh < 0.5)){stop("Do not run branch and bound algorithm with reg.thresh or align.thresh < 0.5 when using non-uniform priors")};
if(bb.alg == T & uniform.priors==T & (reg.thresh < 0.6 | align.thresh < 0.6)){stop("Do not run branch and bound algorithm with reg.thresh or align.thresh < 0.6 when using uniform priors")};
if(! n.cvs %in% c(1,2)){stop("Current version of HyPrMTMC is limited to assessment of a maximum of 2 CVs in a region: set n.cvs to either 1 or 2")};
if(n.cvs == 2 & window.size >= 4E2){print("Computation likely to take a long time. Consider utilising the window.size and sentinel variables to focus assessment within a specified window of a sentinel SNP, default is the lead SNP across all traits")};
if(sensitivity == T){bb.alg = F; print("Performing a regional probability sensitivity assessment: posterior probability of co-localisation will not be computed")};
if(! reg.steps %in% 0:m){stop("reg.steps parameter must be an integer between 0 and number of traits m")};
if(! sense.1 %in% 0:m | ! sense.2 %in% 0:m ){stop("Sensitivity parameters must be intergers between 0 and number of traits m")};
if(sense.1 >= sense.2){stop("Sensitivity parameter sense.2 must be greater than sensitivity parameter sense.1")};
if(window.size > Q){stop("Window size cannot be larger than the number of SNPs Q in the region")}
if((length(trait.subset)<m & typeof(trait.subset)!="integer") & length(trait.names) < m){stop("When using character names in 'trait.subset' the variable 'trait.names' must be of length m, containing names for all traits considered")}
if(bb.selection == "reg.only" & reg.thresh < 0.9){warning("Posterior evaluation and the BB algorithm are based on values of the regional statistic only.\n Consider setting a larger value for the regional threshold, e.g. > 0.9, to avoid difficulties in interpreting any regional association signals.");}
if(class(trait.subset)=="character" & class(trait.names)=="character"){
tmp.traits = c(1:m);
traits = tmp.traits[tolower(trait.names) %in% tolower(trait.subset)];
trait.names = trait.names[traits];
}else if(class(trait.subset)=="integer"){
traits = trait.subset;
trait.names = trait.names[traits];
}else{stop("trait.subset must be character or intger valued")}
if(n.cvs==1 & test.2==F){
W = 0.15/effect.se;
W[, which(binary.outcomes==1)] = 0.2/(effect.se[, which(binary.outcomes==1)]);
}else if(n.cvs==1 & test.2==T){
W = 0.15/effect.se;
W[, which(binary.outcomes==1)] = 0.2/(effect.se[, which(binary.outcomes==1)]);
W.nudge = matrix(0.15 + runif(m*Q, -0.015, 0.015), nrow=Q)/effect.se;
W.nudge[, which(binary.outcomes==1)] = matrix(0.2 + runif(Q*length(which(binary.outcomes==1)), -0.02, 0.02), nrow=Q)/(effect.se[, which(binary.outcomes==1)]);
}else{
if(epsilon == 0){epsilon = 0.01;}
W = matrix(0.15 + runif(m*Q, -0.015, 0.015), nrow=Q)/effect.se;
W[, which(binary.outcomes==1)] = matrix(0.2 + runif(Q*length(which(binary.outcomes==1)), -0.02, 0.02), nrow=Q)/(effect.se[, which(binary.outcomes==1)]);
}
if(n.cvs==1 & (all(trait.cor[lower.tri(trait.cor)] == 0) | ind.traits == TRUE)){
ind.traits = TRUE;
Zsq = Z^2;
Wsq = 1/(1+ W^2);
if(reg.steps == 1 & test.2 == FALSE){rapid = TRUE}else{rapid = FALSE};
}else{
Zsq = Wsq = sparseMatrix(i=1, j=1, dims=c(Q,m));
rapid = FALSE;
}
if(isTRUE(test.2)){
test.2 = 1;
}else{
test.2 = 0;
}
snp.combin = function(x, y, vec){I = iterpc(x, y, labels = vec);return(getall(I)+0.0)};
snp.scores = NA;
if(class(sentinel)=="character"){sentinel = which(snp.id == sentinel);}
tmp.vars = samp.reduc(Z, W, ld.matrix, trait.cor, sample.overlap, traits, window.size, sentinel, Zsq = Zsq, Wsq = Wsq);
Z = tmp.vars[[1]];
W = tmp.vars[[2]];
Zsq = tmp.vars[[7]];
Wsq = tmp.vars[[8]]
ld.matrix = tmp.vars[[3]];
trait.cor = tmp.vars[[4]];
sample.overlap = tmp.vars[[5]];
m = dim(Z)[2];
Q = dim(Z)[1];
snp.id = snp.id[tmp.vars[[6]]];
traits = 1:length(traits);
cor.trts = identical(trait.cor, diag(m));
if(branch.jump == T & cor.trts == F){branch.jump = F};
snp.scores = list();
bb.selection = tolower(bb.selection);
if(n.cvs==1){
if(rapid == F){
snps = c(1:Q);
snps.clc = as.matrix(t(snp.combin(Q, 2, snps)));
flag = 0;
}
if(bb.alg==T){
df = data.frame(matrix(vector(), 0,9, dimnames=list(c(), c("Algorithm_iteration", "Putatively_co_localised_traits", "HyPr_MTC_posterior", "Regional_probability", "Putatively_causal_SNP", "Posterior_explained_by_SNP", "Posterior_ratio_2CVs_vs_1_or_2CVs", "Dropped_trait", "NaNs"))), stringsAsFactors=F)
i=m;
count = 1;
reg.prob = align.prob = 0;
while(i>1){
trts=traits;
while(reg.prob <= reg.thresh | align.prob <= align.thresh){
reg.prob = reg.only = 0;
if(bb.selection == "align"){
tmp.vars = samp.reduc(Z, W, ld.matrix, trait.cor, sample.overlap, trts, Zsq = Zsq, Wsq = Wsq);
if(rapid == F){
reg.result = regional.ABF(tmp.vars[[1]], tmp.vars[[2]], snps.clc, tmp.vars[[3]], tmp.vars[[4]], tmp.vars[[5]], epsilon, reg.thresh, prior.1, prior.2, prior.3, prior.4, flag, 0, reg.steps, cor.adj.priors, uniform.priors, branch.jump, tmp.vars[[7]], tmp.vars[[8]], ind.traits);
reg.res = reg.result[[1]];
reg.full = exp(mpfr(reg.res$log.sum.ABF.full,120));
reg.prob = reg.full/(1 + exp(mpfr(reg.res$log.sum.ABF.all.combs,120)));
}else{
reg.res = rapid.reg(tmp.vars[[7]], tmp.vars[[8]], prior.1, prior.2, uniform.priors)
reg.prob = reg.res[[1]];
}
}else{
while(reg.prob <= reg.thresh){
tmp.vars = samp.reduc(Z, W, ld.matrix, trait.cor, sample.overlap, trts, Zsq = Zsq, Wsq = Wsq);
if(rapid == F){
reg.result = regional.ABF(tmp.vars[[1]], tmp.vars[[2]], snps.clc, tmp.vars[[3]], tmp.vars[[4]], tmp.vars[[5]], epsilon, reg.thresh, prior.1, prior.2, prior.3, prior.4, flag, 0, reg.steps, cor.adj.priors, uniform.priors, branch.jump, tmp.vars[[7]], tmp.vars[[8]], ind.traits);
reg.res = reg.result[[1]];
trts.tmp = as.numeric(unlist(strsplit(reg.res$Traits, split=", ")));
reg.full = exp(mpfr(reg.res$log.sum.ABF.full,120));
reg.prob = reg.full/(1 + exp(mpfr(reg.res$log.sum.ABF.all.combs,120)));
if(reg.prob <= reg.thresh & length(trts) > length(trts.tmp)){
trts = trts[trts.tmp];
}else if(reg.prob <= reg.thresh & length(trts) == length(trts.tmp) & reg.prob < reg.tol){
reg.only.trts = as.numeric(unlist(strsplit(reg.res$reg_only_trts, split=", ")));
trts = trts[reg.only.trts];
}else if(reg.prob <= reg.thresh & length(trts) == length(trts.tmp) & reg.prob >= reg.tol){
break;
}else{
if(length(trts) != length(trts.tmp)){reg.prob = 0;trts = trts[trts.tmp];}
}
}else{
reg.res = rapid.reg(tmp.vars[[7]], tmp.vars[[8]], prior.1, prior.2, uniform.priors)
reg.prob = reg.res[[1]];
if(reg.prob <= reg.thresh){trts = trts[-reg.res[[3]]]}
}
if(length(trts)<=1){reg.only = 1; break;}
}
}
if(length(trts)<=1){break;}
if(bb.selection != "reg.only"){
if(ind.traits == F){
align.res = align.ABF.1(tmp.vars[[1]], tmp.vars[[2]], tmp.vars[[4]], tmp.vars[[5]], tmp.vars[[3]], epsilon, reg.full, align.thresh, prior.1, prior.2, cor.adj.priors, uniform.priors, tmp.vars[[7]], tmp.vars[[8]], ind.traits);
align.prob = asNumeric(reg.full/(reg.full + exp(mpfr(align.res$log.sum.ABF.all.combs, 120))));
if(align.prob <= align.thresh | reg.prob <= reg.thresh){
drop.trait = as.numeric(unlist(strsplit(align.res$Trait.no.clc, split=", ")))
trts = trts[-drop.trait];
}
}else{
align.res = rapid.align(tmp.vars[[7]], tmp.vars[[8]], prior.1, prior.2, uniform.priors);
align.prob = align.res[[1]];
if(align.prob <= align.thresh | reg.prob <= reg.thresh){trts = trts[-align.res[[2]]];}
}
}else{
reg.only = 1; break;
}
if(length(trts)<=1)break
}
df[count,]$Regional_probability = asNumeric(reg.prob);
if(length(trts)>=2){
df[count,]$Algorithm_iteration = count;
df[count,]$Putatively_co_localised_traits = toString(trait.names[trts]);
df[count,]$HyPr_MTC_posterior = asNumeric(reg.prob*align.prob);
if(test.2 == 1){
tmp.vars = samp.reduc(Z, W.nudge, ld.matrix, trait.cor, sample.overlap, trts, Zsq = Zsq, Wsq = Wsq);
reg.tmp = regional.ABF(tmp.vars[[1]], tmp.vars[[2]], snps.clc, tmp.vars[[3]], tmp.vars[[4]], tmp.vars[[5]], epsilon, reg.thresh, prior.1, prior.2, prior.3, prior.4, flag, 1, reg.steps, cor.adj.priors, uniform.priors=FALSE, branch.jump, tmp.vars[[7]], tmp.vars[[8]], ind.traits);
df[count,]$Posterior_ratio_2CVs_vs_1_or_2CVs = asNumeric(reg.tmp[[1]]);
}
if(rapid == F){
if(bb.selection != "reg.only" & ind.traits == F){
df[count,]$NaNs = reg.res$NaNs + align.res$NaNs;
}else{df[count,]$NaNs = reg.res$NaNs;}
df[count,]$Putatively_causal_SNP = toString(snp.id[as.numeric(reg.res$SNPs)]);
df[count,]$Posterior_explained_by_SNP = exp(reg.res$log.max.SNP.ABF.full-reg.res$log.sum.ABF.full);
snp.scores[[count]] = reg.result[[2]];
}else{
df[count,]$Putatively_causal_SNP = toString(snp.id[reg.res[[2]]]);
df[count,]$Posterior_explained_by_SNP = reg.res[[4]][reg.res[[2]]];
snp.scores[[count]] = reg.res[[4]];
}
}else{
df[count,]$Algorithm_iteration = count;
df[count,]$Putatively_co_localised_traits = "None";
df[count,]$Dropped_trait = toString(trait.names[trts]);
if(reg.only == 1 & rapid == F){df[count,]$NaNs = reg.res$NaNs;
}else if(rapid == F){df[count,]$NaNs = reg.res$NaNs + align.res$NaNs;}
}
traits = traits[! traits %in% trts];
reg.prob = align.prob = reg.only = 0;
i = length(traits);
count = count + 1;
}
}else{
if(sensitivity == F){
df = data.frame(matrix(vector(), 0,7, dimnames=list(c(), c("Traits", "HyPr_MTC_posterior", "Regional_probability", "Putatively_causal_SNP", "Posterior_explained_by_SNP", "Posterior_ratio_2CVs_vs_1_or_2CVs", "NaNs"))), stringsAsFactors=F)
df[1,]$Traits = toString(trait.names);
tmp.vars = samp.reduc(Z, W, ld.matrix, trait.cor, sample.overlap, traits, Zsq = Zsq, Wsq = Wsq);
if(rapid == F){
reg.result = regional.ABF(tmp.vars[[1]], tmp.vars[[2]], snps.clc, tmp.vars[[3]], tmp.vars[[4]], tmp.vars[[5]], epsilon, reg.thresh, prior.1, prior.2, prior.3, prior.4, flag, 0, reg.steps, cor.adj.priors, uniform.priors, branch.jump, tmp.vars[[7]], tmp.vars[[8]], ind.traits);
reg.res = reg.result[[1]];
reg.full = exp(mpfr(reg.res$log.sum.ABF.full,120));
reg.prob = asNumeric(reg.full/(1 + exp(mpfr(reg.res$log.sum.ABF.all.combs,120))));
df[1,]$Regional_probability = reg.prob;
if(reg.prob < reg.thresh*align.thresh){
df[1,]$NaNs = reg.res$NaNs;
}else{
if(ind.traits == F){
align.res = align.ABF.1(tmp.vars[[1]], tmp.vars[[2]], tmp.vars[[4]], tmp.vars[[5]], tmp.vars[[3]], epsilon, reg.full, align.thresh, prior.1, prior.2, cor.adj.priors, uniform.priors, tmp.vars[[7]], tmp.vars[[8]], ind.traits);
align.prob = asNumeric(reg.full/(reg.full + exp(mpfr(align.res$log.sum.ABF.all.combs, 120))));
df[1,]$NaNs = reg.res$NaNs + align.res$NaNs;
}else{
align.res = rapid.align(tmp.vars[[7]], tmp.vars[[8]], prior.1, prior.2, uniform.priors);
align.prob = align.res[[1]];
df[1,]$NaNs = reg.res$NaNs;
}
hypr_post = reg.prob*align.prob;
df[1,]$HyPr_MTC_posterior = hypr_post;
if(hypr_post>=(reg.thresh*align.thresh)){
df[1,]$Putatively_causal_SNP = toString(snp.id[as.numeric(reg.res$SNPs)])
df[1,]$Posterior_explained_by_SNP = exp(reg.res$log.max.SNP.ABF.full-reg.res$log.sum.ABF.full);
if(test.2 == 1){
reg.tmp = regional.ABF(tmp.vars[[1]], tmp.vars[[2]], snps.clc, tmp.vars[[3]], tmp.vars[[4]], tmp.vars[[5]], epsilon, reg.thresh, prior.1, prior.2, prior.3, prior.4, flag, 1, reg.steps, cor.adj.priors, uniform.priors = FALSE, branch.jump, tmp.vars[[7]], tmp.vars[[8]], ind.traits);
df[1,]$Posterior_ratio_2CVs_vs_1_or_2CVs = asNumeric(reg.tmp[[1]]);
}
snp.scores = reg.result[[2]];
}
}
}else{
rapid.result = rapid.hyprcoloc(tmp.vars[[7]], tmp.vars[[8]], prior.1, prior.2, uniform.priors);
df[1,]$Regional_probability = rapid.result[[1]];
df[1,]$HyPr_MTC_posterior = rapid.result[[2]];
df[1,]$Putatively_causal_SNP = toString(snp.id[rapid.result[[3]]]);
df[1,]$Posterior_explained_by_SNP = rapid.result[[4]][rapid.result[[3]]];
snp.scores = rapid.result[[4]];
}
}else{
df = data.frame(matrix(vector(), 0,4, dimnames=list(c(), c("Traits", "Regional_Pr_1", "Regional_Pr_2", "Relative_diff"))), stringsAsFactors=F)
df[1,]$Traits = toString(trait.names);
tmp.vars = samp.reduc(Z, W, ld.matrix, trait.cor, sample.overlap, traits, Zsq = Zsq, Wsq = Wsq);
reg.steps = sense.1 ;
reg.result = regional.ABF(tmp.vars[[1]], tmp.vars[[2]], snps.clc, tmp.vars[[3]], tmp.vars[[4]], tmp.vars[[5]], epsilon, reg.thresh, prior.1, prior.2, prior.3, prior.4, flag, 0, reg.steps, cor.adj.priors, uniform.priors, branch.jump, tmp.vars[[7]], tmp.vars[[8]], ind.traits);
reg.res = reg.result[[1]];
reg.full = exp(mpfr(reg.res$log.sum.ABF.full,120));
reg.prob.1 = asNumeric(reg.full/(1 + exp(mpfr(reg.res$log.sum.ABF.all.combs,120))));
df[1,]$Regional_Pr_1 = reg.prob.1;
reg.steps = sense.2;
reg.result = regional.ABF(tmp.vars[[1]], tmp.vars[[2]], snps.clc, tmp.vars[[3]], tmp.vars[[4]], tmp.vars[[5]], epsilon, reg.thresh, prior.1, prior.2, prior.3, prior.4, flag, 0, reg.steps, cor.adj.priors, uniform.priors, branch.jump, tmp.vars[[7]], tmp.vars[[8]], ind.traits);
reg.res = reg.result[[1]];
reg.full = exp(mpfr(reg.res$log.sum.ABF.full,120));
reg.prob.2 = asNumeric(reg.full/(1 + exp(mpfr(reg.res$log.sum.ABF.all.combs,120))));
df[1,]$Regional_Pr_2 = reg.prob.2;
df[1,]$Relative_diff = abs(reg.prob.1-reg.prob.2)/reg.prob.2;
}
}
}else{
flag = 1;
test.2 = 0;
snps = c(1:Q);
snps.clc = as.matrix(t(snp.combin(Q, 2, snps)));
if(bb.alg==T){
NaNs = 0;
i=m;
count = 1;
reg.prob = align.prob = 0;
df = data.frame(matrix(vector(), 0,8, dimnames=list(c(), c("Algorithm_iteration", "Putatively_co_localised_traits", "HyPr_MTC_posterior", "Regional_probability", "Putatively_causal_SNPs", "Posterior_explained_by_SNPs", "Dropped_trait", "NaNs"))), stringsAsFactors=F)
while(i>1){
trts=traits;
tmp.vars = samp.reduc(Z, W, ld.matrix, trait.cor, sample.overlap, trts, Zsq = Zsq, Wsq = Wsq)
while(reg.prob <= reg.thresh | align.prob <= align.thresh){
reg.prob = 0;
while(reg.prob <= reg.thresh){
reg.full = reg.all = reg.only = 0;
reg.result = regional.ABF(tmp.vars[[1]], tmp.vars[[2]], snps.clc, tmp.vars[[3]], tmp.vars[[4]], tmp.vars[[5]], epsilon, reg.thresh, prior.1, prior.2, prior.3, prior.4, flag, test.2, reg.steps, cor.adj.priors, uniform.priors, branch.jump, tmp.vars[[7]], tmp.vars[[8]], ind.traits);
reg.res = reg.result[[1]];
reg.full = sum(exp(mpfr(reg.res$log.sum.ABF.full,120)));
reg.all = sum(exp(mpfr(reg.res$log.sum.ABF.all.combs,120)));
trts.tmp = as.numeric(unlist(strsplit(reg.res$Traits[1], split=", ")));
tmp.mx = which(reg.res$log.sum.ABF.full==max(reg.res$log.sum.ABF.full));
clc.snp = reg.res$SNPs[tmp.mx];
reg.max.snp = reg.res$log.max.SNP.ABF.full[tmp.mx];
reg.max.full = reg.res$log.sum.ABF.full[tmp.mx];
reg.prob = asNumeric(reg.full/(1+reg.all));
if(reg.prob < reg.thresh){
if(length(trts) == length(trts.tmp)){
trts = sample(trts, (length(trts)-1));
}else{trts = trts[trts.tmp];}
if(length(trts)>1){
tmp.vars = samp.reduc(Z, W, ld.matrix, trait.cor, sample.overlap, trts, Zsq = Zsq, Wsq = Wsq);
}
}else{
if(length(trts) != length(trts.tmp)){
reg.prob = 0;
trts = trts[trts.tmp];
tmp.vars = samp.reduc(Z, W, ld.matrix, trait.cor, sample.overlap, trts, Zsq = Zsq, Wsq = Wsq);
}
}
if(length(trts)<=1){reg.only = 1; break;}
}
if(reg.only == 1){break;}
align.res = align.ABF.2(tmp.vars[[1]], tmp.vars[[2]], snps.clc, tmp.vars[[4]], tmp.vars[[5]], tmp.vars[[3]], epsilon, reg.full, align.thresh, prior.1, prior.2, prior.3, prior.4, cor.adj.priors, uniform.priors);
align.tot = exp(mpfr(align.res$log.sum.ABF.align,120));
reg.full.align = sum(exp(mpfr(align.res$log.sum.ABF.co.loc,120)));
if(reg.full.align > reg.full){
clc.snp = align.res$SNPs.clc;
reg.max.snp = exp(mpfr(align.res$log.max.SNP.ABF.full, 120));
reg.max.full = reg.full.align;
tmp.mx = 1;
}
reg.full = reg.full + reg.full.align;
drop.trait = as.numeric(unlist(strsplit(align.res$Trait.no.clc, split=", ")));
align.prob = asNumeric(reg.full/(reg.full + align.tot));
if(align.prob < align.thresh){
trts = trts[-drop.trait];
}
if(length(trts)<=1)break
}
if(reg.only !=1){reg.prob = asNumeric(reg.full/(reg.full.align + reg.all));}
df[count,]$Regional_probability = reg.prob;
if(length(trts)>=2){
NaNs = sum(reg.res$NaNs) + align.res$NaNs
hypr_post = reg.prob*align.prob;
df[count,]$Algorithm_iteration = count;
df[count,]$Putatively_co_localised_traits = toString(trait.names[trts]);
df[count,]$HyPr_MTC_posterior = hypr_post;
df[count,]$Putatively_causal_SNPs = toString(snp.id[as.numeric(unlist(strsplit(clc.snp, split=", ")))]);
df[count,]$Posterior_explained_by_SNPs = asNumeric(reg.max.snp/reg.max.full);
df[count,]$NaNs = NaNs;
snp.scores[[count]] = reg.result[[tmp.mx+1]];
}else{
df[count,]$Algorithm_iteration = count;
df[count,]$Dropped_trait =toString(trait.names[trts]);
if(reg.only == 1){df[count,]$NaNs = sum(reg.res$NaNs);
}else{df[count,]$NaNs = NaNs;}
}
traits = traits[! traits %in% trts];
reg.prob = align.prob = reg.only = 0;
i = length(traits);
count = count + 1;
}
}else{
if(sensitivity == F){
df = data.frame(matrix(vector(), 0,6, dimnames=list(c(), c("Traits", "HyPr_MTC_posterior", "Regional_probability", "Putatively_causal_SNPs", "Posterior_explained_by_SNPs", "NaNs"))), stringsAsFactors=F)
reg.full = reg.all = NaNs = 0;
align.tot = drop.trait = vector("numeric", window.size+2);
df[1,]$Traits = toString(trait.names);
reg.result = regional.ABF(tmp.vars[[1]], tmp.vars[[2]], snps.clc, tmp.vars[[3]], tmp.vars[[4]], tmp.vars[[5]], epsilon, reg.thresh, prior.1, prior.2, prior.3, prior.4, flag, test.2, reg.steps, cor.adj.priors, uniform.priors, branch.jump, tmp.vars[[7]], tmp.vars[[8]], ind.traits);
reg.res = reg.result[[1]];
reg.full = sum(exp(mpfr(reg.res$log.sum.ABF.full,120)));
reg.all = sum(exp(mpfr(reg.res$log.sum.ABF.all.combs,120)));
reg.prob = asNumeric(reg.full/(1 + reg.all));
df[1,]$Regional_probability = reg.prob;
NaNs = sum(reg.res$NaNs)
if(reg.prob < reg.thresh*align.thresh){
df[1,]$NaNs = NaNs;
}else{
trts.tmp = as.numeric(unlist(strsplit(reg.res$Traits[1], split=", ")));
tmp.mx = which(reg.res$log.sum.ABF.full==max(reg.res$log.sum.ABF.full));
clc.snp = reg.res$SNPs[tmp.mx];
reg.max.snp = reg.res$log.max.SNP.ABF.full[tmp.mx];
reg.max.full = reg.res$log.sum.ABF.full[tmp.mx];
align.res = align.ABF.2(tmp.vars[[1]], tmp.vars[[2]], snps.clc, tmp.vars[[4]], tmp.vars[[5]], tmp.vars[[3]], epsilon, reg.full, align.thresh, prior.1, prior.2, prior.3, prior.4, cor.adj.priors, uniform.priors);
align.tot = exp(mpfr(align.res$log.sum.ABF.align,120));
reg.full.align = sum(exp(mpfr(align.res$log.sum.ABF.co.loc,120)));
if(reg.full.align > reg.full){
clc.snp = align.res$SNPs.clc;
reg.max.full = reg.full.align;
tmp.mx = 1;
}
reg.full = reg.full + reg.full.align;
drop.trait = as.numeric(unlist(strsplit(align.res$Trait.no.clc, split=", ")));
NaNs = NaNs + align.res$NaNs;
reg.prob = asNumeric(reg.full/(1 + reg.full.align + reg.all));
align.prob = asNumeric(reg.full/(1 + reg.full + align.tot));
hypr_post = reg.prob*align.prob;
df[1,]$Regional_probability = reg.prob;
df[1,]$HyPr_MTC_posterior = hypr_post;
df[1,]$NaNs = NaNs;
if(hypr_post>=(reg.thresh*align.thresh)){
df[1,]$Putatively_causal_SNPs = toString(snp.id[as.numeric(unlist(strsplit(clc.snp, split=", ")))]);
df[1,]$Posterior_explained_by_SNPs = asNumeric(reg.max.snp/reg.max.full);
snp.scores = reg.result[[tmp.mx + 1]];
}
}
}else{
df = data.frame(matrix(vector(), 0,4, dimnames=list(c(), c("Traits", "Regional_Pr_1", "Regional_Pr_2", "Relative_diff"))), stringsAsFactors=F)
df[1,]$Traits = toString(trait.names);
tmp.vars = samp.reduc(Z, W, ld.matrix, trait.cor, sample.overlap, traits, Zsq = Zsq, Wsq = Wsq);
reg.steps = sense.1;
reg.result = regional.ABF(tmp.vars[[1]], tmp.vars[[2]], snps.clc, tmp.vars[[3]], tmp.vars[[4]], tmp.vars[[5]], epsilon, reg.thresh, prior.1, prior.2, prior.3, prior.4, flag, test.2, reg.steps, cor.adj.priors, uniform.priors, branch.jump, tmp.vars[[7]], tmp.vars[[8]], ind.traits);
reg.res = reg.result[[1]];
reg.full = sum(exp(mpfr(reg.res$log.sum.ABF.full,120)));
reg.all = sum(exp(mpfr(reg.res$log.sum.ABF.all.combs,120)));
reg.prob.1 = asNumeric(reg.full/(1 + reg.all));
df[1,]$Regional_Pr_1 = reg.prob.1;
reg.steps = sense.2;
reg.result = regional.ABF(tmp.vars[[1]], tmp.vars[[2]], snps.clc, tmp.vars[[3]], tmp.vars[[4]], tmp.vars[[5]], epsilon, reg.thresh, prior.1, prior.2, prior.3, prior.4, flag, test.2, reg.steps, cor.adj.priors, uniform.priors, branch.jump, tmp.vars[[7]], tmp.vars[[8]], ind.traits);
reg.res = reg.result[[1]];
reg.full = sum(exp(mpfr(reg.res$log.sum.ABF.full,120)));
reg.all = sum(exp(mpfr(reg.res$log.sum.ABF.all.combs,120)));
reg.prob.2 = asNumeric(reg.full/(1 + reg.all));
df[1,]$Regional_Pr_2 = reg.prob.2;
df[1,]$Relative_diff = abs(reg.prob.1-reg.prob.2)/reg.prob.2;
}
}
}
if(length(which(df$NaNs!=0)) > 1){print("Results include NAs. Possible the Z score correlation matrix or adjusted prior matrix has a high condition number somewhere: consider increasing epsilon by a small amount");}
names(df)[names(df)=="Algorithm_iteration"] <- "iteration"
names(df)[names(df)=="Putatively_co_localised_traits"] <- "traits"
names(df)[names(df)=="Traits"] <- "traits"
names(df)[names(df)=="HyPr_MTC_posterior"] <- "posterior_prob"
names(df)[names(df)=="Regional_probability"] <- "regional_prob"
names(df)[names(df)=="Posterior_ratio_2CVs_vs_1_or_2CVs"] <- "posterior_ratio_2cvs_vs_1_or_2cvs"
names(df)[names(df)=="Putatively_causal_SNP"] <- "candidate_snp"
names(df)[names(df)=="Putatively_causal_SNPs"] <- "candidate_snp"
names(df)[names(df)=="Posterior_explained_by_SNP"] <- "posterior_explained_by_snp"
names(df)[names(df)=="Posterior_explained_by_SNPs"] <- "posterior_explained_by_snp"
names(df)[names(df)=="Dropped_trait"] <- "dropped_trait"
names(df)[names(df)=="Regional_Pr_1"] <- "regional_prob_1"
names(df)[names(df)=="Regional_Pr_2"] <- "regional_prob_2"
names(df)[names(df)=="Relative_diff"] <- "relative_diff"
df$NaNs <- NULL
rownames(df) <- NULL
if(test.2==F){df$posterior_ratio_2cvs_vs_1_or_2cvs <- NULL}
if(any(!is.na(df$posterior_prob))){df$posterior_prob[!is.na(df$posterior_prob)] <- round(as.numeric(df$posterior_prob[!is.na(df$posterior_prob)]),4); df$regional_prob[!is.na(df$regional_prob)] <- round(as.numeric(df$regional_prob[!is.na(df$regional_prob)]),4); df$posterior_explained_by_snp[!is.na(df$posterior_explained_by_snp)] <- round(as.numeric(df$posterior_explained_by_snp[!is.na(df$posterior_explained_by_snp)]),4)}
if(any(!is.na(df$regional_prob))){df$regional_prob[!is.na(df$regional_prob)] <- round(as.numeric(df$regional_prob[!is.na(df$regional_prob)]),4); df$posterior_explained_by_snp[!is.na(df$posterior_explained_by_snp)] <- round(as.numeric(df$posterior_explained_by_snp[!is.na(df$posterior_explained_by_snp)]),4)}
if(length(df$posterior_ratio_2cvs_vs_1_or_2cvs)>0){if(any(!is.na(df$posterior_ratio_2cvs_vs_1_or_2cvs))){df$posterior_ratio_2cvs_vs_1_or_2cvs[!is.na(df$posterior_ratio_2cvs_vs_1_or_2cvs)] <- round(as.numeric(df$posterior_ratio_2cvs_vs_1_or_2cvs[!is.na(df$posterior_ratio_2cvs_vs_1_or_2cvs)]),4)}}
if(sum(is.na(snp.scores))==0 & length(snp.scores)>0 & snpscores==T){results = list(results=df, snpscores=snp.scores)}else{results = list(results=df)};
class(results) <- "hyprcoloc";
return(results)
}
#' Print HyPrColoc
#'
#' print method for class "hyprcoloc"
#' @param x an object of class "hyprcoloc"
#' @author Christopher N Foley (University of Cambridge) <[email hidden]> and James R Staley (University of Bristol) <[email hidden]>
#' @param ... Other arguments passed to or from other methods
#'
#' @export print.hyprcoloc
#' @export
print.hyprcoloc <- function(x, ...){
cat("\nCall: \nhyprcoloc")
cat("\n\nResults:\n")
print(x$results)
cat("\n")
}
# prevents R check throwing errors about unassigned global variables
utils::globalVariables(c("colorRampPalette", "combn", "runif"))

hyprcoloc.R at commit 26ea595, no license · at the source

Overview

Authors: Weiming Gong^1,2, Yantong Guo^1,2, Xiubin Sun^1,2, Shukang Wang^1,2, Fuzhong Xue^1,2,3, Wei Zhao⁴, Xiang Zhou^5,6, Lu Liu^1,2, Zhongshang Yuan^1,2

ORCID iDs: Fuzhong Xue, Wei Zhao, Xiang Zhou, Lu Liu, Zhongshang Yuan

Department of Biostatistics, School of Public Health, Cheeloo College of Medicine, Shandong University,Jinan, China
Institute for Medical Dataology, Cheeloo College of Medicine, Shandong University,Jinan, China
Qilu Hospital, Cheeloo College of Medicine, Shandong University,Jinan, China
Department of Pathogenic Biology, Key Laboratory of Infection and Immunity of Shandong Province, and Key Laboratory for Experimental Teratology of the Chinese Ministry of Education, School of Basic Medical Science, Cheeloo College of Medicine, Shandong University,Jinan, China
Department of Biostatistics, University of Michigan,Ann Arbor, MI USA
Center for Statistical Genetics, University of Michigan,Ann Arbor, MI USA

Institutions: Shandong University (China); University of Michigan (United States)

Journal: Nature communications, volume 17, issue 1, article 5482

Dates: received 7 May 2025; accepted 9 April 2026; published online 21 April 2026

Type: Research article · Language: English

License: CC BY-NC-ND

Identifiers: DOI 10.1038/s41467-026-72164-7 · PMID 42014408 · PMCID PMC13284391 · OpenAlex W7155065151

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: genetics / omics (modality), human (organism), mouse (organism)

Methods: Statistics, Preprocessing

Keywords: Genetic association study, Musculoskeletal abnormalities, Computational biology and bioinformatics

MeSH: Genetic Predisposition to Disease*, Musculoskeletal Diseases*, Animals, Autoimmune Diseases, Female, Genome-Wide Association Study, Humans, Mice, Multifactorial Inheritance, Multivariate Analysis, Polymorphism, Single Nucleotide (* major topic)

Topic: Genetic Associations and Epidemiology (Genetics, Biochemistry, Genetics and Molecular Biology), according to OpenAlex

Funding: National Natural Science Foundation of China (National Science Foundation of China) (82373686); Shandong Provincial Natural Science Foundation (ZR2024JQ029); Taishan Scholar Project of Shandong Province (tsqn202211025)

Citations: cited by 1 paper (Europe PMC); 126 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repositories

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

Zenodo 19012445

License: CC-BY-4.0

State: the link answers, verified on 29 September 2026

Evidence: files inventoried

Size: 4 files

Software Heritage: not checked

Found in: the references

Holds: README

Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation

Availability: 1 check, the latest on 29 September 2026: the link answers (HTTP 200)

29 September 2026: the link answers (HTTP 200)

1 file

README.txt, Text, 23 lines

At the source:

GenomicSEM/GenomicSEM

License: GPL-3.0

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: 6b65ca5db39fdade08b0d811477be1cdd57b5039, 26 August 2026

Languages: R (32)

Size: 82 files, 32 scripts

Software Heritage: not archived

Found in: “Code availability”

Holds: README, license file, environment (DESCRIPTION), documentation

Not found: CITATION.cff, tests, continuous integration

Tools: tidyverse (4 files), data.table (3 files), ggplot2 (2 files), ggpubr (2 files), lavaan (2 files), psych (1 file)

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

34 files

R/addGenes.R, R, 143 lines
R/addSNPs.R, R, 267 lines
R/commonfactor.R, R, 403 lines, 2 matches
R/commonfactorGWAS.R, R, 276 lines
R/commonfactorGWAS_main.R, R, 216 lines
R/enrich.R, R, 611 lines
R/hdl.R, R, 524 lines
R/indexS.R, R, 51 lines
R/ldsc.R, R, 583 lines
R/localSRMD.R, R, 22 lines
R/multiGene.R, R, 349 lines
R/multiSNP.R, R, 364 lines
R/munge.R, R, 109 lines
R/munge_main.R, R, 138 lines
R/paLDSC.R, R, 380 lines, 2 matches
R/qtrait.r, R, 738 lines
R/read_fusion.R, R, 99 lines
R/rgmodel.R, R, 1,063 lines, 1 match
R/s_ldsc.R, R, 851 lines
R/simLDSC.R, R, 225 lines
R/subSV.R, R, 77 lines
R/summaryGLS.R, R, 57 lines
R/summaryGLSbands.R, R, 271 lines
R/sumstats.R, R, 142 lines
R/sumstats_main.R, R, 226 lines
R/userGWAS.R, R, 490 lines
R/userGWAS_main.R, R, 403 lines
R/userGWASa.r, R, 388 lines
R/usermodel.R, R, 685 lines
R/utils.R, R, 226 lines
R/utils_sanitychecks.R, R, 67 lines
R/write.model.R, R, 108 lines
LICENSE, License, 621 lines
README.md, Text, 99 lines

At the source: github.com/GenomicSEM/GenomicSEM

josefin-werme/LAVA

License: none: the authors keep all their rights

State: the link answers, verified on 29 September 2026

Evidence: files inventoried

Commit: e729a245f7b6923967a96804fbf5246eadf2d6c6, 28 July 2026

Languages: R (12), C/C++ (5), C++ (3)

Size: 56 files, 20 scripts

Software Heritage: not archived

Found in: “Code availability”

Holds: README, environment (DESCRIPTION), documentation, 2 notebooks

Not found: license file, CITATION.cff, tests, continuous integration

Tools: data.table (2 files)

Availability: 1 check, the latest on 29 September 2026: the link answers

29 September 2026: the link answers

21 files

R/alignment.R, R, 65 lines
R/analysis_functions.R, R, 411 lines
R/binary_processing.R, R, 88 lines
R/confidence_intervals.R, R, 109 lines
R/cpp11.R, R, 17 lines
R/input_eqtl.r, R, 398 lines, 2 matches
R/input_processing.R, R, 692 lines
R/on_load.R, R, 5 lines
R/partial_corrs.R, R, 64 lines
R/wishart.R, R, 155 lines
README.Rmd, R, 432 lines
src/cpp11.cpp, C++, 51 lines
src/global.h, C/C++, 48 lines
src/load_ld.cpp, C++, 346 lines
src/load_ld.h, C/C++, 75 lines
src/load_ld.tpl.h, C/C++, 65 lines
src/load_plink.cpp, C++, 226 lines
src/load_plink.h, C/C++, 66 lines
src/utils.h, C/C++, 139 lines
vignettes/sample_overlap.Rmd, R, 66 lines
README.md, Text, 741 lines

At the source: github.com/josefin-werme/LAVA

bulik/ldsc

License: GPL-3.0

State: the link answers, verified on 26 September 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: “Code availability”

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: 1 check, the latest on 26 September 2026: the link answers

26 September 2026: the link answers

27 files

ContinuousAnnotations/quantile_M.pl, Perl, 241 lines
ContinuousAnnotations/quantile_h2g.r, R, 76 lines
ldsc.py, Python, 660 lines
ldscore/__init__.py, Python, 1 line
ldscore/irwls.py, Python, 196 lines
ldscore/jackknife.py, Python, 514 lines
ldscore/ldscore.py, Python, 415 lines
ldscore/parse.py, Python, 292 lines
ldscore/regressions.py, Python, 743 lines
ldscore/sumstats.py, Python, 581 lines
make_annot.py, Python, 56 lines
munge_sumstats.py, Python, 745 lines
setup.py, Python, 20 lines
test/parse_test/test.l2.M, MATLAB, 1 line
test/parse_test/test1.l2.M, MATLAB, 1 line
test/parse_test/test2.l2.M, MATLAB, 1 line
test/parse_test/test_bad.l2.M, MATLAB, 1 line
test/simulate.py, Python, 81 lines
test/test_irwls.py, Python, 69 lines
test/test_jackknife.py, Python, 267 lines
test/test_ldscore.py, Python, 111 lines
test/test_munge_sumstats.py, Python, 358 lines
test/test_parse.py, Python, 129 lines
test/test_regressions.py, Python, 342 lines
test/test_sumstats.py, Python, 487 lines
LICENSE, License, 675 lines
README.md, Text, 122 lines

At the source: github.com/bulik/ldsc

cnfoley/hyprcoloc

License: none: the authors keep all their rights

State: the link answers, verified on 29 September 2026

Evidence: files inventoried

Commit: 26ea5953a46b3e204dfa8eadd202f746244afa13, 12 February 2021

Languages: C++ (8), R (3)

Size: 45 files, 11 scripts

Software Heritage: archived

Found in: “Code availability”

Holds: README, environment (DESCRIPTION), documentation, 1 notebook

Not found: license file, CITATION.cff, tests, continuous integration

Tools: pheatmap (1 file)

Availability: 1 check, the latest on 29 September 2026: the link answers

29 September 2026: the link answers

12 files

R/RcppExports.R, R, 31 lines
R/hyprcoloc.R, R, 1,451 lines, 2 matches
src/RcppExports.cpp, C++, 137 lines
src/align1.cpp, C++, 150 lines
src/align12.cpp, C++, 255 lines
src/align1ind.cpp, C++, 97 lines
src/align2.cpp, C++, 335 lines
src/regional1.cpp, C++, 112 lines
src/regional1ind.cpp, C++, 73 lines
src/regional2.cpp, C++, 218 lines
vignettes/hyprcoloc.Rmd, R, 453 lines
README.md, Text, 44 lines

At the source: github.com/cnfoley/hyprcoloc

chr1swallace.github.io/coloc

License: none: the authors keep all their rights

State: the link answers, verified on 29 September 2026

Evidence: the link answers

Software Heritage: not checked

Found in: “Code availability”

Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation

Availability: 1 check, the latest on 29 September 2026: the link answers (HTTP 200)

29 September 2026: the link answers (HTTP 200)

At the source: chr1swallace.github.io/coloc/

JianYang-Lab/gsMap

License: MIT

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: 7fbdcf80cb3e35a341cf1a955bb4809f3e5e6e02, 19 August 2026

Languages: Python (37)

Size: 95 files, 37 scripts

Software Heritage: not archived

Found in: “Code availability”

Holds: README, license file, environment (docker-compose.yml, Dockerfile, pyproject.toml, docs/requirements.txt, gsmap_vibe/pyproject.toml, visualization_web_docs/requirements.txt, src/gsMap/setup.py), tests, continuous integration, documentation

Not found: CITATION.cff

Tools: pandas (14 files), NumPy (13 files), SciPy (10 files), PyTorch (6 files), Scanpy (6 files), scikit-learn (3 files), anndata (2 files), Plotly (2 files), Numba (1 file), PyTorch Geometric (1 file)

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

39 files

docs/source/conf.py, Python, 53 lines
gsmap_vibe/gsmap_mcp/__init__.py, Python, 1 line
gsmap_vibe/gsmap_mcp/__main__.py, Python, 3 lines
gsmap_vibe/gsmap_mcp/docs.py, Python, 19 lines
gsmap_vibe/gsmap_mcp/server.py, Python, 276 lines
gsmap_vibe/gsmap_mcp/tools.py, Python, 178 lines
src/gsMap/GNN/__init__.py, Python, 1 line
src/gsMap/GNN/adjacency_matrix.py, Python, 73 lines
src/gsMap/GNN/model.py, Python, 92 lines
src/gsMap/GNN/train.py, Python, 85 lines
src/gsMap/__init__.py, Python, 5 lines
src/gsMap/__main__.py, Python, 4 lines
src/gsMap/cauchy_combination_test.py, Python, 144 lines
src/gsMap/config.py, Python, 1,494 lines
src/gsMap/create_slice_mean.py, Python, 169 lines
src/gsMap/diagnosis.py, Python, 365 lines
src/gsMap/find_latent_representation.py, Python, 157 lines
src/gsMap/format_sumstats.py, Python, 445 lines
src/gsMap/generate_ldscore.py, Python, 1,018 lines
src/gsMap/latent_to_gene.py, Python, 322 lines
src/gsMap/main.py, Python, 46 lines
src/gsMap/report.py, Python, 174 lines, 1 match
src/gsMap/run_all_mode.py, Python, 241 lines, 1 match
src/gsMap/setup.py, Python, 5 lines
src/gsMap/spatial_ldsc_multiple_sumstats.py, Python, 439 lines
src/gsMap/utils/__init__.py, Python, 1 line
src/gsMap/utils/generate_r2_matrix.py, Python, 610 lines
src/gsMap/utils/jackknife.py, Python, 518 lines
src/gsMap/utils/manhattan_plot.py, Python, 620 lines
src/gsMap/utils/regression_read.py, Python, 175 lines
src/gsMap/utils/torch_utils.py, Python, 23 lines
src/gsMap/visualize.py, Python, 217 lines
tests/conftest.py, Python, 469 lines
tests/test_advanced_usage.py, Python, 183 lines
tests/test_cli.py, Python, 69 lines
tests/test_docs_cli_parsing.py, Python, 125 lines
visualization_web_docs/source/conf.py, Python, 56 lines
LICENSE, License, 21 lines
README.md, Text, 130 lines

At the source: github.com/JianYang-Lab/gsMap

Code availability statement

The paper has a code availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

it points to the authors' code: chr1swallace.github.io/coloc, bulik/ldsc, cnfoley/hyprcoloc, GenomicSEM/GenomicSEM, JianYang-Lab/gsMap, josefin-werme/LAVA

Read it in the paper: doi.org/10.1038/s41467-026-72164-7.

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:

7 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
125 scripts, each with its path and the digest of its content;
11 matches 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.

Code and data availability statement

The paper has a code and data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

it points to the authors' code: chr1swallace.github.io/coloc, bulik/ldsc, cnfoley/hyprcoloc, GenomicSEM/GenomicSEM, JianYang-Lab/gsMap, josefin-werme/LAVA

Read it in the paper: doi.org/10.1038/s41467-026-72164-7.

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, 29 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 9 authors, 3 keywords, 11 MeSH terms, 3 funders, 126 references.

Cite

This paper

Gong, W., Guo, Y., Sun, X., Wang, S., Xue, F., Zhao, W., Zhou, X., Liu, L., & Yuan, Z. (2026). Multivariate genetic analysis reveals three distinct pathological dimensions in musculoskeletal disorders. Nature communications, 17(1), 5482. https://doi.org/10.1038/s41467-026-72164-7

BibTeX

@article{gong2026multivariate, author = {Gong, Weiming and Guo, Yantong and Sun, Xiubin and Wang, Shukang and Xue, Fuzhong and Zhao, Wei and Zhou, Xiang and Liu, Lu and Yuan, Zhongshang}, title = {{Multivariate genetic analysis reveals three distinct pathological dimensions in musculoskeletal disorders}}, journal = {Nature communications}, year = {2026}, month = apr, volume = {17}, number = {1}, pages = {5482}, publisher = {Nature Publishing Group}, issn = {2041-1723}, doi = {10.1038/s41467-026-72164-7}, url = {https://doi.org/10.1038/s41467-026-72164-7}, pmid = {42014408}, pmcid = {PMC13284391} }

RIS

TY - JOUR AU - Gong, Weiming AU - Guo, Yantong AU - Sun, Xiubin AU - Wang, Shukang AU - Xue, Fuzhong AU - Zhao, Wei AU - Zhou, Xiang AU - Liu, Lu AU - Yuan, Zhongshang TI - Multivariate genetic analysis reveals three distinct pathological dimensions in musculoskeletal disorders T2 - Nature communications J2 - Nat Commun PY - 2026 DA - 2026/04/21 VL - 17 IS - 1 SP - 5482 SN - 2041-1723 PB - Nature Publishing Group DO - 10.1038/s41467-026-72164-7 UR - https://doi.org/10.1038/s41467-026-72164-7 LA - en ER -

CSL-JSON

{ "id": "10.1038/s41467-026-72164-7", "type": "article-journal", "title": "Multivariate genetic analysis reveals three distinct pathological dimensions in musculoskeletal disorders", "container-title": "Nature communications", "author": [ { "family": "Gong", "given": "Weiming" }, { "family": "Guo", "given": "Yantong" }, { "family": "Sun", "given": "Xiubin" }, { "family": "Wang", "given": "Shukang" }, { "family": "Xue", "given": "Fuzhong" }, { "family": "Zhao", "given": "Wei" }, { "family": "Zhou", "given": "Xiang" }, { "family": "Liu", "given": "Lu" }, { "family": "Yuan", "given": "Zhongshang" } ], "container-title-short": "Nat Commun", "volume": "17", "issue": "1", "page": "5482", "DOI": "10.1038/s41467-026-72164-7", "PMID": "42014408", "PMCID": "PMC13284391", "ISSN": "2041-1723", "publisher": "Nature Publishing Group", "URL": "https://doi.org/10.1038/s41467-026-72164-7", "language": "en", "issued": { "date-parts": [ [ 2026, 4, 21 ] ] } }

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/s41467-026-73428-y [code]: Regional heterogeneity in phenotypic and genetic associations between bone and brain in humans.
Journal: Nature communications
In common: lavaan, BEDTools, psych, 7 other tools, genetics / omics, 9 references
[2] doi:10.1038/s43856-026-01510-z [code]: Mapping genetic convergence across brain structure, mental health, and cardiometabolic disease.
Journal: Communications medicine
In common: lavaan, BEDTools, psych, 8 other tools, genetics / omics, 4 references
[3] doi:10.1038/s41467-026-71803-3 [code]: Charting the transition from in vitro gliogenesis to the in vivo maturation of human glial progenitor cells transplanted into the hypomyelinated mouse brain.
Journal: Nature communications
In common: BEDTools, anndata, Numba, 12 other tools, genetics / omics, mouse
[4] doi:10.1038/s41562-026-02486-5 [code]: Genome-wide association studies of infant and toddler temperament in European and multi-ancestry populations.
Journal: Nature human behaviour
In common: psych, data.table, ggplot2, 4 other tools, genetics / omics, 8 references
[5] doi:10.1038/s41592-026-03194-8 [code]: Beyond benchmarking: an expert-guided consensus approach to spatially aware clustering.
Journal: Nature methods
In common: PyTorch Geometric, anndata, Scanpy, 10 other tools, genetics / omics, 1 reference
[6] doi:10.1186/s13059-026-04177-w [code]: Genomic sequence evolution underlying human neocortical interareal diversification.
Journal: Genome biology
In common: BEDTools, anndata, Scanpy, 10 other tools, genetics / omics, mouse
[7] doi:10.1016/j.isci.2026.116055 [code]: Mapping the transcriptional diversity of calcium signaling in the mouse and human brain.
Journal: iScience
In common: PyTorch Geometric, anndata, Numba, 10 other tools, genetics / omics, mouse
[8] doi:10.1016/j.xcrm.2026.102766 [code]: A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: anndata, Numba, Scanpy, 10 other tools, genetics / omics
[9] doi:10.1038/s41467-026-71682-8 [code]: GWAS meta-analysis of cerebrospinal fluid Alzheimer's biomarkers reveals loci regulating lipids, brain volume and autophagy.
Journal: Nature communications
In common: BEDTools, data.table, ggplot2, 5 other tools, genetics / omics, 6 references
[10] doi:10.1002/imt2.70163 [code]: Spatial multi-omics unveils sphingolipid metabolic reprogramming within the retinal pathological niche.
Journal: iMeta
In common: anndata, Scanpy, pheatmap, 9 other tools, genetics / omics, mouse, 1 reference

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.

Multivariate genetic analysis reveals three distinct pathological dimensions in musculoskeletal disorders.

Code ↔ Paper

Paper

The authors' code

Overview

Abstract

Repositories

Zenodo 19012445

GenomicSEM/GenomicSEM

josefin-werme/LAVA

bulik/ldsc

cnfoley/hyprcoloc

chr1swallace.github.io/coloc

JianYang-Lab/gsMap

Code availability statement

Tracing map

Data

Code and data availability statement

Versions

Version 1, 29 September 2026: the first record

Cite

This paper

Similar papers

Contribute

Claim this paper

Correct its record

Validate its tracing map

Add the badge to its README

Request its removal

Discussion, reproductions, activity