Introduction

Cell-cell communication is usually inferred from the co-expression of a ligand in one cell type and its receptor in another. Co-expression, however, is not evidence of a causal signal: shared regulators, tissue composition and technical factors all induce correlations between cell types without any signalling taking place. MRCCC implements MR-CCC, a Bayesian Mendelian randomization framework that uses cis-eQTL genotypes as instrumental variables for ligand expression in a sender cell type and for receptor expression in a receiver cell type, and asks whether the genetically predicted ligand level, modulated by the genetically predicted receptor level, shifts the activity of a downstream pathway in the receiver.

The unit of analysis is a ligand-receptor-pathway triplet. For each triplet the package returns a posterior inclusion probability (PIP) that summarises the evidence for a causal signal, together with effect estimates, credible intervals, posterior sign probabilities, the receptor level at which the ligand effect reverses sign, and the Monte Carlo and convergence diagnostics needed to trust these numbers. Across many triplets, the PIPs can be thresholded with Bayesian false discovery rate control.

library(MRCCC)

The model

For donor $i$, let $X_i$ be ligand expression in the sender cell type, $Z_i$ receptor expression in the receiver cell type, and $Y_i$ pathway activity in the receiver cell type. Let $G_i$ and $H_i$ be the cis-eQTL genotypes that instrument the ligand and the receptor, and $V_i$ a vector of shared donor covariates. The first stage predicts each exposure from its instruments, $X_i^* = E[X_i \mid G_i, V_i]$ and $Z_i^* = E[Z_i \mid H_i, V_i]$. The second stage models the outcome as \[ Y_i = \mu + \gamma\,(\beta_X X_i^* + \beta_{XZ} X_i^* Z_i^*) + \beta_Z Z_i^* + \alpha_Y^\top V_i + \varepsilon_{Yi}, \] where $\gamma \in \{0, 1\}$ is a spike-and-slab inclusion indicator on the causal block $(\beta_X, \beta_{XZ})$. Because $X^*$ and $Z^*$ depend only on genotype and the observed covariates, they are unaffected by the unmeasured confounders that distort a naive regression of $Y$ on $X$ and $Z$. All three stages are estimated jointly by a blocked Gibbs sampler.

The posterior inclusion probability $\mathrm{PIP} = P(\gamma = 1 \mid \text{data})$ is the communication score. The receptor-modulated ligand effect is $\beta_X + \beta_{XZ} z$, which changes sign at $\tau = -\beta_X / \beta_{XZ}$: when $\tau$ lies inside the observed receptor range, the ligand activates the pathway in some donors and suppresses it in others. The model requires $X$, $Z$ and $Y$, and every column of $G$, $H$ and $V$, to be centred, because neither first stage has an intercept; mr_ccc() does this internally and reports effects on a standardised scale (per standard deviation of $X$ and $Z$, in standard deviations of $Y$) alongside the raw scale.

Simulating a triplet

simulate_mrccc() generates one triplet from the working model with an unmeasured confounder that loads on the ligand, the receptor and the pathway. The default generates 300 donors, five instruments per gene, three covariates, and a positive interaction.

sim <- simulate_mrccc(n = 300, beta_X = 0.3, beta_XZ = 0.3, beta_Z = 0.5,
                      conf_strength = 0.7, seed = 1)
str(sim[c("X", "Z", "Y", "G", "H", "V")])
#> List of 6
#>  $ X: num [1:300, 1] -2.28 2.15 -1.6 -1.74 1.53 ...
#>  $ Z: num [1:300, 1] 2.07 0.288 0.844 1.701 -2.609 ...
#>  $ Y: num [1:300, 1] -1.413 1.894 -0.349 1.717 -2.794 ...
#>  $ G: num [1:300, 1:5] -0.626 0.184 -0.836 1.595 0.33 ...
#>  $ H: num [1:300, 1:5] 0.85 -0.925 0.894 -0.941 0.539 ...
#>  $ V: num [1:300, 1:3] 0.739 0.387 1.296 -0.804 -1.603 ...
sim$truth
#> $beta_X
#> [1] 0.3
#> 
#> $beta_XZ
#> [1] 0.3
#> 
#> $beta_Z
#> [1] 0.5
#> 
#> $pi_val
#> [1] 0.5
#> 
#> $conf_strength
#> [1] 0.7
#> 
#> $tau
#> [1] -1

Fitting the model

mr_ccc() validates the inputs, centres $X$, $Z$, $Y$ and the columns of $G$, $H$ and $V$, runs the Gibbs sampler, and pools the retained draws. The chain length below is kept short so that this document builds quickly; the section on chain length discusses the settings appropriate for reported results.

fit <- mr_ccc(sim$X, sim$Z, sim$Y, sim$G, sim$H, sim$V,
              n_iter = 4000, burn_in = 500, seed = 1)
#> The Monte Carlo uncertainty of the PIP could not be estimated (too few retained draws per chain, or a constant gamma chain), so the borderline check was skipped.
fit
#> MR-CCC fit
#>   Chains: 4; retained draws per chain: 3500 (n_iter = 4000, burn_in = 500, thin = 1)
#>   Posterior inclusion probability: 1  (MCSE NA, ESS NA)
#>   First-stage F: X = 62.8 (5 instruments), Z = 55.5 (5 instruments)
#>   Standardised posterior means: beta_X = 0.202, beta_XZ = 0.485, beta_Z = 0.345

summary() prints the PIP with its Monte Carlo standard error, the standardised effects with credible intervals, and the convergence diagnostics, and returns the assembled table invisibly.

tab <- summary(fit)
#> MR-CCC fit summary
#>   Chains: 4; retained draws per chain: 3500
#>   Posterior inclusion probability: 1  (MCSE NA, ESS NA)
#>   Discovery threshold: 0.5
#>   First-stage F: X = 62.8, Z = 55.5
#> 
#> Standardised effects (posterior mean and 95% credible interval):
#>   beta_X   0.202  [0.0706, 0.338]
#>   beta_XZ  0.485  [0.31, 0.671]
#>   beta_Z   0.345  [0.207, 0.488]
#> 
#> Gelman-Rubin R-hat (NA unless n_chains >= 2):
#>   beta_X   R-hat = 1
#>   beta_XZ  R-hat = 1
#>   gamma    R-hat = NA
#>   loglik   R-hat = 1
tab
#>   parameter  raw_mean  std_mean  std_lower std_upper     rhat
#> 1    beta_X 0.2577800 0.2023687 0.07062601 0.3382584 0.999891
#> 2   beta_XZ 0.3385379 0.4853670 0.30957803 0.6710859 1.000091
#> 3    beta_Z 0.4299449 0.3452156 0.20695710 0.4878785       NA
#> 4       pip 1.0000000 1.0000000         NA        NA       NA

Input from a SummarizedExperiment

Donor-level (pseudo-bulk) expression is commonly stored in a SummarizedExperiment, with genes as rows and donors (samples) as columns. se_row() extracts one gene from such an object as a named numeric vector that can be passed to mr_ccc() as X, Z or Y. A SummarizedExperiment that contains exactly one gene can also be passed to mr_ccc() directly, in which case its first assay row is used. The chunk below stores the simulated ligand expression in a one-gene SummarizedExperiment and refits the model from it, with a short chain to keep this document fast.

m <- rbind(LIG1 = as.numeric(sim$X))
colnames(m) <- paste0("donor", seq_len(nrow(sim$G)))
se <- SummarizedExperiment::SummarizedExperiment(assays = list(logcounts = m))
x_se <- se_row(se, "LIG1")
fit_se <- mr_ccc(x_se, sim$Z, sim$Y, sim$G, sim$H, sim$V,
                 n_iter = 1500, burn_in = 300, seed = 1)
#> Warning: Gelman-Rubin R-hat exceeds 1.01 for beta_X, beta_XZ, gamma, loglik:
#> the chains have not mixed. Increase 'n_iter' and 'burn_in'.
fit_se$pip
#> [1] 0.8327083

Posterior summaries

Credible intervals

credible_intervals(fit)
#>   parameter      lower    median     upper
#> 1    beta_X 0.07062601 0.2011509 0.3382584
#> 2   beta_XZ 0.30957803 0.4838821 0.6710859
#> 3    beta_Z 0.20695710 0.3444264 0.4878785
credible_intervals(fit, level = 0.9, standardized = FALSE)
#>   parameter     lower    median     upper
#> 1    beta_X 0.1189013 0.2562288 0.4021894
#> 2   beta_XZ 0.2336778 0.3375022 0.4473497
#> 3    beta_Z 0.2856495 0.4289619 0.5782621

Sign probabilities of the interaction

The prior on $(\beta_X, \beta_{XZ})$ is a spike-and-slab mixture, so the marginal posterior of $\beta_{XZ}$ mixes draws taken while $\gamma = 0$ (concentrated near zero) with draws taken while $\gamma = 1$. The marginal sign probabilities combine the evidence that there is an effect with the evidence about its direction; the conditional probabilities describe the direction given that the ligand acts on the pathway. Both are informative.

sign_probabilities(fit)
#>   pip p_pos_marginal p_neg_marginal p_pos_conditional p_neg_conditional
#> 1   1              1              0                 1                 0
#>   n_conditional
#> 1         14000

Sign-reversal threshold

sign_reversal() computes $\tau = -\beta_X / \beta_{XZ}$ from paired draws and reports it in standard-deviation units of receptor expression relative to the mean. The element p_in_range is the posterior probability that $\tau$ lies inside the observed receptor range, which determines whether the reversal is realised in the donor population.

sr <- sign_reversal(fit, Z_observed = sim$Z)
sr[c("tau_median", "tau_lower", "tau_upper", "Z_range_std", "p_in_range")]
#> $tau_median
#> [1] -0.4163952
#> 
#> $tau_lower
#> [1] -0.8180225
#> 
#> $tau_upper
#> [1] -0.1422925
#> 
#> $Z_range_std
#> [1] -2.766094  2.656162
#> 
#> $p_in_range
#> [1] 1

In this simulation the true threshold is $Z = -0.3 / 0.3 = -1$, which corresponds to -0.54 standard deviations from the mean receptor expression.

Effect curve

plot_effect_curve() draws the posterior mean of $\beta_X + \beta_{XZ} z$ over the observed receptor range with a pointwise credible band built from paired draws.

plot_effect_curve(fit, Z_observed = sim$Z)

Receptor-modulated ligand effect with a 95 percent pointwise credible band.

Choosing the chain length

The PIP is the mean of an autocorrelated binary chain, so its Monte Carlo standard error (MCSE) is estimated by non-overlapping batch means within each chain and combined across chains in quadrature. The effective sample size (ESS) is the number of independent draws that would give the same precision.

c(pip = fit$pip, mcse = fit$pip_mcse, ess = fit$pip_ess)
#>  pip mcse  ess 
#>    1   NA   NA

When the PIP lies within two MCSEs of the discovery threshold (0.5 by default), mr_ccc() issues a warning: the discovery status of the triplet is not resolved at the current chain length, and either n_iter or n_chains should be increased. The warning is a property of the Monte Carlo approximation, not of the data; a longer run narrows the MCSE until the status is resolved or the PIP is shown to sit genuinely at the threshold.

The short chain used in this document is adequate for illustration only. For reported results, substantially longer chains should be used; the analysis in the paper pools four chains of 100,000 iterations.

Chains alternate their starting value of $\gamma$ and initialise $(\beta_X, \beta_{XZ})$ from a normal distribution with standard deviation init_scale, so they begin at dispersed points. The Gelman-Rubin statistic $\hat R$ is then reported for $\beta_X$, $\beta_{XZ}$, $\gamma$ and the joint log-likelihood in fit$diagnostics; it is NA if a single chain was run, because $\hat R$ compares the variance between chains with the variance within them and there is nothing to compare. This is why the default is n_chains = 4.

Pooling chains costs little more than lengthening one. Per-chain Monte Carlo standard errors combine in quadrature, so $C$ chains of n_iter iterations, each discarding burn_in, carry approximately the same standard error as one chain that retains $C \times$ (n_iter $-$ burn_in) draws. The additional cost is the $C - 1$ extra burn-in periods, and several chains also yield $\hat R$.

fit2 <- mr_ccc(sim$X, sim$Z, sim$Y, sim$G, sim$H, sim$V,
               n_iter = 2000, burn_in = 500, n_chains = 2, seed = 2)
#> The Monte Carlo uncertainty of the PIP could not be estimated (too few retained draws per chain, or a constant gamma chain), so the borderline check was skipped.
fit2$diagnostics
#>   parameter      rhat
#> 1    beta_X 1.0003549
#> 2   beta_XZ 0.9996861
#> 3     gamma        NA
#> 4    loglik 1.0021389

Screening many triplets

In a full analysis mr_ccc() is applied to every candidate triplet and the PIPs are collected in a vector. For a declared set $S$ the Bayesian false discovery rate is \[ \mathrm{bFDR}(S) = |S|^{-1} \sum_{j \in S} (1 - \mathrm{PIP}_j), \] the posterior expected proportion of false discoveries. bayes_fdr() returns the bFDR of each nested top-$k$ set ranked by PIP, and bfdr_threshold_set() selects the largest set whose bFDR does not exceed a target level.

pips <- c(0.98, 0.95, 0.91, 0.84, 0.72, 0.61, 0.48, 0.30, 0.12, 0.05)
data.frame(pip = pips, bfdr = bayes_fdr(pips),
           selected_10 = bfdr_threshold_set(pips, alpha = 0.10),
           selected_20 = bfdr_threshold_set(pips, alpha = 0.20))
#>     pip       bfdr selected_10 selected_20
#> 1  0.98 0.02000000        TRUE        TRUE
#> 2  0.95 0.03500000        TRUE        TRUE
#> 3  0.91 0.05333333        TRUE        TRUE
#> 4  0.84 0.08000000        TRUE        TRUE
#> 5  0.72 0.12000000       FALSE        TRUE
#> 6  0.61 0.16500000       FALSE        TRUE
#> 7  0.48 0.21571429       FALSE       FALSE
#> 8  0.30 0.27625000       FALSE       FALSE
#> 9  0.12 0.34333333       FALSE       FALSE
#> 10 0.05 0.40400000       FALSE       FALSE

Session information

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] MRCCC_0.99.4     BiocStyle_2.41.0
#> 
#> loaded via a namespace (and not attached):
#>  [1] sass_0.4.10                 generics_0.1.4             
#>  [3] SparseArray_1.13.4          lattice_0.23-1             
#>  [5] digest_0.6.39               evaluate_1.0.5             
#>  [7] grid_4.6.1                  RColorBrewer_1.1-3         
#>  [9] fastmap_1.2.0               jsonlite_2.0.0             
#> [11] Matrix_1.7-6                BiocManager_1.30.27        
#> [13] scales_1.4.0                jquerylib_0.1.4            
#> [15] abind_1.4-8                 cli_3.6.6                  
#> [17] rlang_1.3.0                 XVector_0.53.0             
#> [19] Biobase_2.73.2              withr_3.0.3                
#> [21] cachem_1.1.0                DelayedArray_0.39.8        
#> [23] yaml_2.3.12                 otel_0.2.0                 
#> [25] S4Arrays_1.13.2             tools_4.6.1                
#> [27] ggplot2_4.0.3               SummarizedExperiment_1.43.0
#> [29] BiocGenerics_0.59.12        buildtools_1.0.0           
#> [31] vctrs_0.7.3                 R6_2.6.1                   
#> [33] matrixStats_1.5.0           stats4_4.6.1               
#> [35] lifecycle_1.0.5             Seqinfo_1.3.2              
#> [37] S4Vectors_0.51.10           IRanges_2.47.5             
#> [39] bslib_0.12.0                gtable_0.3.6               
#> [41] glue_1.8.1                  Rcpp_1.1.2                 
#> [43] xfun_0.61                   GenomicRanges_1.65.4       
#> [45] sys_3.4.3                   MatrixGenerics_1.25.0      
#> [47] knitr_1.52                  farver_2.1.2               
#> [49] htmltools_0.5.9             labeling_0.4.3             
#> [51] rmarkdown_2.32              maketools_1.3.2            
#> [53] compiler_4.6.1              S7_0.2.2

Introduction to MRCCC