Causal ASCEND infers which downstream molecular variables are causal ancestors of which, using an upstream layer that is known to precede them. This vignette shows how to run it on matrices, SummarizedExperiment and MultiAssayExperiment objects, how to read its output, and which settings matter.
causalASCEND 0.99.5
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
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.
causalASCEND works directly on Bioconductor data containers. ascend()
accepts a SummarizedExperiment for each layer, and ascendMultiAssay()
takes a MultiAssayExperiment and the names of its background and
foreground experiments, matching samples across them. Results are
returned as an S4 AscendResult object, whose accessors give the
ancestral relation as a matrix or an edge list, so it can be passed to
graph and network packages for further analysis.
Bioconductor already offers gene regulatory network inference, for example GENIE3 (tree-based regulator-to-target scores) and minet (mutual-information networks such as ARACNE, CLR and MRNET), and causal reasoning over a known network with CausalR. causalASCEND differs from these in two ways. It estimates causal ancestry (whether one gene acts upstream of another, directly or through intermediates) rather than association scores or direct edges, and it uses a second, upstream data layer, such as genotypes, to orient those relations. Unlike CausalR, it needs no prior network.
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("causalASCEND")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)
set.seed(1)
sim <- simulateTwoTier(n = 2000, nBackground = 20, nForeground = 15,
r2 = 0.7, sparsity = 0.9, pCross = 0.2, effect = 0.9)
dim(sim$background)
#> [1] 20 2000
dim(sim$foreground)
#> [1] 15 2000ascend() 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, 0 non-descendant, 86 unrelated, 0 unresolved
#> 3283 CI tests in 2 sweep(s)
#> accessors: ancestralMatrix(), ancestralEdges(), pairPValues(), ascendStats()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 x8 ancestor R1 1 4 5.101074e-50
#> 2 x3 x11 ancestor R1 1 5 4.953271e-43
#> 3 x3 x1 ancestor R1 1 2 5.160937e-04
#> 4 x3 x2 ancestor R1 1 3 5.815322e-07
#> 5 x3 x9 ancestor R1 1 3 3.443731e-06
#> 6 x3 x4 ancestor transitivity 1 NA 2.407359e-02
table(edges$relation, edges$rule)
#>
#> R1 transitivity
#> ancestor 15 4The 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 unrelated
#> 19 86The 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 NABecause the data are simulated, we can score the result against the true graph.
evaluateAncestral(res, trueAncestral(sim))
#> $tp
#> [1] 15
#>
#> $fp
#> [1] 4
#>
#> $fn
#> [1] 0
#>
#> $tn
#> [1] 191
#>
#> $unresolved
#> [1] 0
#>
#> $precision
#> [1] 0.7894737
#>
#> $recall
#> [1] 1
#>
#> $f1
#> [1] 0.8823529
#>
#> $orientationAccuracy
#> [1] 1
#>
#> $coverage
#> [1] 1ascend() 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':
#>
#> Filter, Find, Map, Position, Reduce, anyDuplicated, aperm, append,
#> as.data.frame, basename, cbind, colnames, dirname, do.call,
#> duplicated, eval, evalq, get, grep, grepl, is.unsorted, lapply,
#> mapply, match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
#> rank, rbind, 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':
#>
#> I, expand.grid, 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] TRUEFor 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] TRUEalpha 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
#> 15ascendStats() 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] 3283
#>
#> $nPairTests
#> [1] 105
#>
#> $nWitnessTests
#> [1] 772
#>
#> $nBlanketTests
#> [1] 2406
#>
#> $meanPairCondSize
#> [1] 11.90476
#>
#> $maxPairCondSize
#> [1] 19
#>
#> $nSweeps
#> [1] 2sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.5 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_GB LC_COLLATE=C
#> [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: America/New_York
#> 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.5
#> [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 R6_2.6.1
#> [10] XVector_0.53.0 S4Arrays_1.13.2 knitr_1.52
#> [13] DelayedArray_0.39.8 bookdown_0.48 bslib_0.12.0
#> [16] rlang_1.3.0 cachem_1.1.0 xfun_0.62
#> [19] sass_0.4.10 otel_0.3.0 SparseArray_1.13.4
#> [22] cli_3.6.6 digest_0.6.39 grid_4.6.1
#> [25] lifecycle_1.0.5 evaluate_1.0.5 abind_1.4-8
#> [28] rmarkdown_2.32 tools_4.6.1 htmltools_0.5.9