Introduction

Many multi-omic studies measure two layers of the same samples, where one layer is known to act before the other: genotypes or DNA methylation (the background) come before gene expression (the foreground). Causal ASCEND (Ancestral Scalable Causal discovEry via iNherited Descent; ASCEND for short) uses that order to infer the causal order among the foreground variables. For every pair of foreground variables it decides whether one is a causal ancestor of the other, whether they are causally unrelated, or leaves the pair unresolved.

Each conditional independence test conditions only on the nearest ancestors of the two variables, which ASCEND learns from the data and refines over several sweeps, rather than on the whole background layer. This keeps every test small, so the method scales to thousands of background variables.

ASCEND assumes that

no foreground variable causes a background variable, and
the data are Markov and faithful to a directed acyclic graph.

Its tests are Fisher-z partial-correlation tests, so it targets linear relations. Background variables that were not measured act as hidden confounders and lower its accuracy.

Installation

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("causalASCEND")

A first analysis

We simulate a two-tier system with 20 background and 15 foreground variables measured on 2,000 samples. Matrices have features in rows and samples in columns, as in other Bioconductor containers.

library(causalASCEND)
sim <- simulateTwoTier(n = 2000, nBackground = 20, nForeground = 15,
    r2 = 0.7, sparsity = 0.9, pCross = 0.2, effect = 0.9, seed = 1)
dim(sim$background)
#> [1]   20 2000
dim(sim$foreground)
#> [1]   15 2000

ascend() takes the two layers and returns an AscendResult.

res <- ascend(sim$background, sim$foreground)
res
#> AscendResult: 15 foreground and 20 background variables, 2000 samples
#>   105 pairs: 19 ancestor, 1 non-descendant, 85 unrelated, 0 unresolved
#>   3289 CI tests in 2 sweep(s)
#>   accessors: ancestralMatrix(), ancestralEdges(), pairPValues(), ascendStats()

Reading the result

ancestralEdges() lists the relations ASCEND committed to, with the rule that produced each one and its evidence.

edges <- ancestralEdges(res)
head(edges)
#>   from  to relation         rule sweep votes       pValue
#> 1   x3 x11 ancestor           R1     1     5 4.953271e-43
#> 2   x3  x2 ancestor           R1     1     3 5.815322e-07
#> 3   x3  x8 ancestor           R1     1     4 5.101074e-50
#> 4   x3  x9 ancestor           R1     1     3 3.443731e-06
#> 5   x3  x1 ancestor           R1     1     2 5.160937e-04
#> 6   x3  x4 ancestor transitivity     1    NA 2.407359e-02
table(edges$relation, edges$rule)
#>                 
#>                  R1 R2 transitivity
#>   ancestor       15  0            4
#>   non-descendant  0  1            0

The relations are:

ancestor: from is a causal ancestor of to (a direct or indirect cause);
non-descendant: from is not caused by to, but ASCEND could not show that it causes to;
unrelated and unresolved pairs are returned on request:

all <- ancestralEdges(res, include = c("ancestor", "non-descendant",
    "unrelated", "unresolved"))
table(all$relation)
#> 
#>       ancestor non-descendant      unrelated 
#>             19              1             85

The rules are R1 (deactivation, gives an ancestor relation), R2 (activation, gives a non-descendant relation), R3 (conditional independence, gives an unrelated call) and transitivity (if A causes B and B causes C, A causes C). votes counts the witnesses that agreed on an R1 or R2 call; pValue is the p-value of the pair’s last independence test, a continuous score that can rank pairs, for example to compute a precision-recall curve.

The full matrix is available too: entry [a, b] is 1 if a is an ancestor of b, 0.5 if a is not a descendant of b, 0 if the two are unrelated and NA if unresolved.

round(ancestralMatrix(res)[1:6, 1:6], 1)
#>     x3 x5 x6 x10 x12 x13
#> x3  NA  0  0   0   0   0
#> x5   0 NA  0   0   0   0
#> x6   0  0 NA   0   0   0
#> x10  0  0  0  NA   0   0
#> x12  0  0  0   0  NA   0
#> x13  0  0  0   0   0  NA

Because the data are simulated, we can score the result against the true graph.

evaluateAncestral(res, trueAncestral(sim))
#> $tp
#> [1] 15
#> 
#> $fp
#> [1] 5
#> 
#> $fn
#> [1] 0
#> 
#> $tn
#> [1] 190
#> 
#> $unresolved
#> [1] 0
#> 
#> $precision
#> [1] 0.75
#> 
#> $recall
#> [1] 1
#> 
#> $f1
#> [1] 0.8571429
#> 
#> $orientationAccuracy
#> [1] 1
#> 
#> $coverage
#> [1] 1

Bioconductor containers

ascend() also accepts SummarizedExperiment objects; choose the assay with backgroundAssay and foregroundAssay.

library(SummarizedExperiment)
#> Loading required package: MatrixGenerics
#> Loading required package: matrixStats
#> 
#> Attaching package: 'MatrixGenerics'
#> The following objects are masked from 'package:matrixStats':
#> 
#>     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
#>     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
#>     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
#>     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
#>     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
#>     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
#>     colWeightedMeans, colWeightedMedians, colWeightedSds,
#>     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
#>     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
#>     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
#>     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
#>     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
#>     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
#>     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
#>     rowWeightedSds, rowWeightedVars
#> Loading required package: GenomicRanges
#> Loading required package: stats4
#> Loading required package: BiocGenerics
#> Loading required package: generics
#> 
#> Attaching package: 'generics'
#> The following objects are masked from 'package:base':
#> 
#>     as.difftime, as.factor, as.ordered, intersect, is.element, setdiff,
#>     setequal, union
#> 
#> Attaching package: 'BiocGenerics'
#> The following objects are masked from 'package:stats':
#> 
#>     IQR, mad, sd, var, xtabs
#> The following object is masked from 'package:utils':
#> 
#>     data
#> The following objects are masked from 'package:base':
#> 
#>     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
#>     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
#>     get, grep, grepl, is.unsorted, lapply, Map, mapply, match, mget,
#>     order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
#>     rbind, Reduce, rownames, sapply, saveRDS, scale, sequence, table,
#>     tapply, transform, unique, unsplit, which.max, which.min
#> Loading required package: S4Vectors
#> 
#> Attaching package: 'S4Vectors'
#> The following object is masked from 'package:utils':
#> 
#>     findMatches
#> The following objects are masked from 'package:base':
#> 
#>     expand.grid, I, unname
#> Loading required package: IRanges
#> Loading required package: Seqinfo
#> Loading required package: Biobase
#> Welcome to Bioconductor
#> 
#>     Vignettes contain introductory material; view with
#>     'browseVignettes()'. To cite Bioconductor, see
#>     'citation("Biobase")', and for packages 'citation("pkgname")'.
#> 
#> Attaching package: 'Biobase'
#> The following object is masked from 'package:MatrixGenerics':
#> 
#>     rowMedians
#> The following objects are masked from 'package:matrixStats':
#> 
#>     anyMissing, rowMedians
geno <- SummarizedExperiment(assays = list(dosage = sim$background))
expr <- SummarizedExperiment(assays = list(logcounts = sim$foreground))
res_se <- ascend(geno, expr, backgroundAssay = "dosage",
    foregroundAssay = "logcounts")
identical(ancestralMatrix(res_se), ancestralMatrix(res))
#> [1] TRUE

For a MultiAssayExperiment, name the two experiments. ascendMultiAssay() matches samples through the sample map and keeps those present in both.

library(MultiAssayExperiment)
mae <- MultiAssayExperiment(experiments = list(snp = sim$background,
    rna = sim$foreground))
res_mae <- ascendMultiAssay(mae, background = "snp", foreground = "rna")
identical(ancestralMatrix(res_mae), ancestralMatrix(res))
#> [1] TRUE

Settings

alpha and alphaMB are the significance levels of the orientation rules and of the nearest-ancestor searches (both 0.05).
fdr = TRUE applies a Benjamini-Hochberg correction to the pairwise tests of each sweep.
minVotes is the number of agreeing witnesses needed before a direction is committed. Raising it gives fewer, more reliable directions.
prescreen limits each foreground variable’s background candidates to those marginally associated with it (Benjamini-Hochberg threshold 0.3). This bounds the cost when the background is large; set it to 1 to keep every background variable.

strict <- ascend(sim$background, sim$foreground, minVotes = 2)
table(ancestralEdges(strict)$relation)
#> 
#> ancestor 
#>       15

Cost

ascendStats() counts the conditional independence tests by type and reports the size of the conditioning sets. Most tests are spent finding nearest ancestors, and the pairwise conditioning sets stay far smaller than the background layer.

st <- ascendStats(res)
st[c("nTests", "nPairTests", "nWitnessTests", "nBlanketTests",
    "meanPairCondSize", "maxPairCondSize", "nSweeps")]
#> $nTests
#> [1] 3289
#> 
#> $nPairTests
#> [1] 105
#> 
#> $nWitnessTests
#> [1] 772
#> 
#> $nBlanketTests
#> [1] 2412
#> 
#> $meanPairCondSize
#> [1] 11.90476
#> 
#> $maxPairCondSize
#> [1] 19
#> 
#> $nSweeps
#> [1] 2

Limitations

The two-tier order must be known. Mislabelling a few variables does little harm, but feedback from the foreground to the background does.
Unmeasured background variables act as hidden confounders.
The Fisher-z test targets linear relations. ASCEND’s rules use only the outcomes of the tests, so other conditional independence tests can be used in principle; this version implements Fisher-z.
In dense graphs, orientation errors are mostly reversals; directions are most reliable in sparse networks.

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] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] MultiAssayExperiment_1.39.1 SummarizedExperiment_1.43.0
#>  [3] Biobase_2.73.2              GenomicRanges_1.65.4       
#>  [5] Seqinfo_1.3.2               IRanges_2.47.5             
#>  [7] S4Vectors_0.51.10           BiocGenerics_0.59.12       
#>  [9] generics_0.1.4              MatrixGenerics_1.25.0      
#> [11] matrixStats_1.5.0           causalASCEND_0.99.0        
#> [13] BiocStyle_2.41.0           
#> 
#> loaded via a namespace (and not attached):
#>  [1] Matrix_1.7-6        jsonlite_2.0.0      compiler_4.6.1     
#>  [4] BiocManager_1.30.27 jquerylib_0.1.4     yaml_2.3.12        
#>  [7] fastmap_1.2.0       lattice_0.23-1      XVector_0.53.0     
#> [10] R6_2.6.1            S4Arrays_1.13.2     knitr_1.52         
#> [13] DelayedArray_0.39.8 maketools_1.3.2     bslib_0.12.0       
#> [16] rlang_1.3.0         cachem_1.1.0        xfun_0.61          
#> [19] sass_0.4.10         sys_3.4.3           otel_0.2.0         
#> [22] SparseArray_1.13.4  cli_3.6.6           digest_0.6.39      
#> [25] grid_4.6.1          lifecycle_1.0.5     evaluate_1.0.5     
#> [28] buildtools_1.0.0    abind_1.4-8         rmarkdown_2.32     
#> [31] tools_4.6.1         htmltools_0.5.9

Two-tier causal discovery with causalASCEND