| Type: | Package |
| Title: | Cell Explanatory Power with Inflection Point (CEP-IP) with CMRF |
| Version: | 0.1.0 |
| Description: | Implements the cell explanatory power with inflection point (CEP-IP) framework for single-cell co-expression analysis. A generalized additive model (GAM) fitted by restricted maximum likelihood (REML) with a penalized residual sum of squares (PRSS) search over the basis dimension relates a gene of interest to the mean of its partner genes. Cells are ranked by their explanatory power, and the inflection point is the first negative-to-positive zero crossing of a residual GAM. Optional extensions are the Spearman-Kendall cell-level score (SKCS) and the copula monotone recovery filter (CMRF), which recovers cells near the module through a tolerance band and a repeated stratified k-fold cross-validation gate. Downstream helpers give per-donor module fractions, paired donor-level differential expression and donor-level association models. CEP-IP is described in Wong (2026) <doi:10.1016/j.cmpb.2026.109372>. |
| License: | MIT + file LICENSE |
| Encoding: | UTF-8 |
| LazyData: | true |
| Depends: | R (≥ 4.1.0) |
| Imports: | mgcv, stats, utils, methods, grDevices, grid, parallel, Matrix, ggplot2 (≥ 3.4.0), gridExtra, openxlsx |
| Suggests: | Hmisc, energy, clinfun, pcaPP, testthat (≥ 3.0.0), knitr, rmarkdown |
| VignetteBuilder: | knitr |
| Config/testthat/edition: | 3 |
| RoxygenNote: | 7.3.1 |
| URL: | https://github.com/kahkengwong/cepip |
| BugReports: | https://github.com/kahkengwong/cepip/issues |
| NeedsCompilation: | no |
| Packaged: | 2026-10-19 23:24:38 UTC; root |
| Author: | Kah Keng Wong |
| Maintainer: | Kah Keng Wong <kahkeng3@gmail.com> |
| Repository: | CRAN |
| Date/Publication: | 2026-10-29 19:30:19 UTC |
cepip: CEP-IP with CMRF for single-cell co-expression
Description
The main function is [cepip()]. The steps can also be run one by one: [cepip_screen()], [cepip_composite()], [cepip_fit()], [cepip_ep()], [cepip_ip()], [cepip_skcs()], [cepip_band()], [cepip_cmrf()] and [cepip_validate()], then [cepip_donor()], [cepip_deg()] and [cepip_associate()]. [cepip_save()] writes the tables and plots of one analysis to a folder.
Author(s)
Maintainer: Kah Keng Wong kahkeng3@gmail.com (ORCID)
See Also
Useful links:
Run CEP-IP, optionally with SKCS and CMRF
Description
Runs the whole method on one population of cells.
Usage
cepip(
x,
y,
donor = NULL,
id = NULL,
skcs = FALSE,
cmrf = FALSE,
validate = cmrf,
ip_k = "prss",
ip_axis = c("y", "x"),
band = NULL,
kappa = 0.5,
bw = 0.1,
y_spine = NULL,
workers = 1L,
repeats = 40L,
tests = TRUE,
verbose = TRUE
)
Arguments
x |
Expression of the gene of interest (GOI), one value per cell. |
y |
Composite expression of the partner genes (see [cepip_composite()]). |
donor |
Optional donor label per cell. |
id |
Optional cell identifiers (for example barcodes). |
skcs |
Add the SKCS step. |
cmrf |
Add the CMRF steps (implies 'skcs = TRUE'). |
validate |
Run Validation1 after CMRF. |
ip_k |
k of the GAM f: '"prss"' (default) or a fixed integer. The published SEA-AD analysis used 10. |
ip_axis |
'"y"' or '"x"', the axis that splits Pre-IP from Post-IP cells. |
band |
Optional 'cepip_band' object to use instead of building the band. |
kappa, bw |
Band settings (see [cepip_band()]). |
y_spine |
Optional composite values for the spine fit (see [cepip_band()]). |
workers |
Number of parallel workers (GAM k searches and the RSk-fCV gate). |
repeats |
Number of RSk-fCV repeats of the gate (see [cepip_cmrf()]). |
tests |
Run the CMRF class tests. |
verbose |
Print progress messages. |
Details
* Always: the CEP-IP GAM with the PRSS k search ([cepip_fit()]), the EP and TREP cells ([cepip_ep()]) and the inflection point with its IPRS ([cepip_ip()]). With 'skcs = FALSE' and 'cmrf = FALSE' this is the published CEP-IP, except that cells are labelled Pre-IP or Post-IP on the response axis ('ip_axis = "y"'); set 'ip_axis = "x"' for the published split on the GOI axis. * 'skcs = TRUE': the Spearman-Kendall cell-level scores, the SKCS-core cells and the SKCS-core fit ([cepip_skcs()]). * 'cmrf = TRUE': also SKCS, the CMRF band ([cepip_band()]), the RSk-fCV gate ([cepip_cmrf()]) and, by default, Validation1 ([cepip_validate()]). CMRF always uses the CEP-IP inflection point.
When no IP is found, the run stops after the IP step and the object reports why; no error is raised.
Value
An object of class 'cepip'.
Examples
d <- cepip_toy
res <- cepip(d$x, d$y, donor = d$donor, id = d$id, ip_k = 6)
res
# With SKCS, CMRF and Validation1 (a short gate of 5 repeats)
res2 <- cepip(d$x, d$y, donor = d$donor, id = d$id, cmrf = TRUE, repeats = 5,
tests = FALSE)
res2
Donor-level association with outcomes
Description
For each outcome, fits 'outcome ~ predictor + covariates' and 'outcome ~ covariates' by ordinary least squares on the donors with complete data, after z-scoring every variable. It reports the standardized beta of the predictor, its p value, and delta R2, the R2 of the full model minus that of the covariate-only model. An outcome is never used as its own covariate. The BH false discovery rate is applied across all outcomes passed.
Usage
cepip_associate(
data,
predictor,
outcomes,
covariates = character(0),
min_donors = 20L
)
Arguments
data |
Donor table (one row per donor). |
predictor |
Column name of the predictor. |
outcomes |
Column names of the outcomes. |
covariates |
Column names of the covariates. |
min_donors |
Minimum number of complete donors for a model. |
Details
In the SEA-AD analysis the predictor was the Module Post-IP percentage per donor, the covariates were age at death, sex, PMI and APOE e4 dose, and the BH family was 15 outcomes (nine clinical and pathology outcomes plus six control variables). Outcomes and covariates depend on the data available.
Value
A data frame with one row per outcome.
Examples
d <- cepip_toy
res <- cepip(d$x, d$y, donor = d$donor, id = d$id, verbose = FALSE)
dt <- cepip_donor(res)
set.seed(1)
dt$age <- rnorm(nrow(dt), 80, 5)
dt$score <- 0.05 * dt$pct_PostIP + rnorm(nrow(dt))
cepip_associate(dt, "pct_PostIP", outcomes = "score", covariates = "age",
min_donors = 10)
CMRF tolerance band (Stage 1)
Description
Builds the CMRF Stage 1 band in copula space. This step is fast, so it can be run several times to choose 'kappa' and 'bw' (see [plot_band()]) before running the slower Stage 2 gate in [cepip_cmrf()]. It is optional: [cepip_cmrf()] builds the band itself with 'kappa = 0.5' and 'bw = 0.1' when no band is given.
Usage
cepip_band(
skcs,
ip,
kappa = 0.5,
bw = 0.1,
coverage = 0.95,
spine_k = "prss",
min_spine_cells = 50L,
min_donors = 3L,
y_spine = NULL,
workers = 1L,
verbose = FALSE
)
Arguments
skcs |
A 'cepip_skcs' object (with its SKCS-core fit). |
ip |
A 'cepip_ip' object of the same cells. |
kappa |
Height of the widening at IP_u, as a fraction of delta_base. |
bw |
Bandwidth (standard deviation) of the Gaussian kernel, in u units. |
coverage |
Quantile of |v - g(u)| that sets delta_base. |
spine_k |
'"prss"' or a fixed k for the spine. |
min_spine_cells |
Minimum number of SKCS-core TREP cells for the spine. |
min_donors |
Minimum number of donors for the donor random intercept. |
y_spine |
Optional composite values (one per cell, same order) whose ranks give the v coordinates of the spine fit. By default the composite of 'skcs' is used. The published SEA-AD analysis fitted the spine on the composite column of its data export, which is stored in single precision; passing that column reproduces the published spine exactly. |
workers |
Number of parallel workers for the spine k search. |
verbose |
Print progress messages. |
Details
* Copula coordinates over all cells: u = rank(x) / (n + 1) and v = rank(y) / (n + 1), with average ranks for ties. * IP_u = (number of cells with x <= IPx) / (n + 1) and IP_v = (number of cells with y <= IPy) / (n + 1). * The spine g(u) is a GAM, v ~ s(u, bs = "tp", k) (with the donor random intercept and inverse-donor weights when donors are known), fitted by REML ('select = TRUE', 'gamma = 1.5') on the SKCS-core TREP cells, with k chosen by PRSS as in [cepip_fit()]. The donor term is excluded when predicting. * delta_base is the 95 |v - g(u)| over all SKCS-core cells. * delta_adapt(u) = delta_base x (1 + kappa x K(u)), with the Gaussian kernel K(u) = exp(-(u - IP_u)^2 / (2 bw^2)), which is 1 at IP_u. The band is therefore widest at the IP, and at IP_u +- bw the widening term has decayed to exp(-0.5) = 60.65 * A cell passes Stage 1 when |v - g(u)| <= delta_adapt(u). SKCS-core cells that pass are the reference set (Original + Refined, OriRef); non-SKCS cells that pass are the Stage 2 candidates; SKCS-core cells that fail are "Original only" and leave the Module.
Value
An object of class 'cepip_band'.
Examples
d <- cepip_toy
s <- cepip_skcs(d$x, d$y, donor = d$donor, id = d$id)
ip <- cepip_ip(d$x, d$y, donor = d$donor, id = d$id)
b <- cepip_band(s, ip)
b$delta_base
table(b$cells$is_candidate)
Copula monotone recovery filter (CMRF)
Description
Runs CMRF Stage 2, the repeated stratified k-fold cross-validation (RSk-fCV) gate, on the Stage 1 band of [cepip_band()].
Usage
cepip_cmrf(
skcs = NULL,
ip = NULL,
band = NULL,
kappa = 0.5,
bw = 0.1,
folds = 5L,
repeats = 40L,
alpha = 0.05,
seed = 42L,
workers = 1L,
tests = TRUE,
nperm_jt = 10000L,
nperm_hd = 9999L,
R_dcor = 999L,
verbose = FALSE
)
Arguments
skcs, ip |
'cepip_skcs' and 'cepip_ip' objects of the same cells. Not needed when 'band' is given. |
band |
Optional 'cepip_band' object. When it is missing, the band is built with 'kappa' and 'bw'. |
kappa, bw |
Band settings used when 'band' is missing. |
folds, repeats |
Number of folds and of repeats of RSk-fCV. |
alpha |
Significance level of the gate. |
seed |
Base seed of the fold assignment. |
workers |
Number of parallel workers for the repeats. |
tests |
Run the class tests (needs the Suggests packages clinfun, Hmisc and energy). |
nperm_jt, nperm_hd, R_dcor |
Permutations of the Jonckheere-Terpstra, Hoeffding D and distance correlation tests. |
verbose |
Print progress messages. |
Details
* Curve A is an isotonic (PAVA) fit of v on u for the reference cells (OriRef) alone; curve B is the same fit for OriRef plus the active candidates. Each curve is cross-validated with 'folds' folds stratified into deciles of u, repeated 'repeats' times, which gives 'folds' x 'repeats' test RMSE values per curve. * The gate compares B with A by a one-sided paired Wilcoxon test (B greater than A) and by dRMSE = median(B - A) / median(A) against the background variability BV = IQR(A) / median(A). * While p < 'alpha' and dRMSE > BV, the active candidate with the largest |v - PAVA(u)| (PAVA fitted once on OriRef plus all candidates) is removed and curve B is refitted. The candidates left when the gate closes are the Recovered cells.
Classes: Original + Refined (OriRef; SKCS-core, inside the band), Recovered (non-SKCS, inside the band, passed the gate), Original only (SKCS-core, outside the band) and Excluded (all other cells, including candidates removed by the gate). The Module is OriRef plus Recovered.
The fold assignment of repeat r uses the seed 'seed' + r x 1000 + stratum, so results do not depend on 'workers'.
Value
An object of class 'cepip_cmrf'.
Examples
d <- cepip_toy
s <- cepip_skcs(d$x, d$y, donor = d$donor, id = d$id)
ip <- cepip_ip(d$x, d$y, donor = d$donor, id = d$id)
b <- cepip_band(s, ip)
# a short gate of 5 repeats (the default is 40)
m <- cepip_cmrf(band = b, repeats = 5, tests = FALSE)
m$counts
Composite expression of partner genes
Description
The composite is the arithmetic mean of the partner genes per cell, on the scale of the input values.
Usage
cepip_composite(expr, partners)
Arguments
expr |
Matrix or data frame of expression values (cells x genes). |
partners |
Names (or column indices) of the partner genes. |
Value
A numeric vector, one value per cell.
Examples
d <- cepip_toy
y <- cepip_composite(d, paste0("P", 1:8))
head(y)
Differential expression between cell groups, paired by donor
Description
Compares two of the four groups of [cepip_groups()] (foreground against background). These four contrasts are built in: 'PostIP_vs_PreIP', 'NMPost_vs_NMPre', 'PostIP_vs_NMPost' and 'PreIP_vs_NMPre'. Other pairs of the four groups can be given as '"A_vs_B"'.
Usage
cepip_deg(
x,
expr,
contrast = c("PostIP_vs_PreIP", "NMPost_vs_NMPre", "PostIP_vs_NMPost",
"PreIP_vs_NMPre"),
unit = c("donor", "cell"),
min_pct = 0.05,
min_cells_donor = 3L,
min_donors = 5L,
fdr = 0.05,
lfc = 0.25,
min_cells = 10L,
verbose = FALSE
)
Arguments
x |
A 'cepip', 'cepip_cmrf' or 'cepip_ip' object. |
expr |
Log1p-normalized expression, genes x cells, with cell ids as column names (the ids used in the analysis) and gene names as row names. |
contrast |
One or more contrasts. |
unit |
'"donor"' or '"cell"'. |
min_pct, min_cells_donor, min_donors, fdr, lfc |
See Details. |
min_cells |
Minimum number of cells per group. |
verbose |
Print progress messages. |
Details
Donor unit (default): 1. Genes expressed (value > 0) in at least 'min_pct' of the foreground or of the background cells are kept. 2. Donors with at least 'min_cells_donor' cells in both groups are paired; at least 'min_donors' paired donors are needed. 3. For each donor and gene, log2FC = log2(mean(expm1(x_fg)) + 1) - log2(mean(expm1(x_bg)) + 1), on log1p-normalized values. 4. Genes that change in fewer than 5 gene is tested by a Wilcoxon signed-rank test of its non-zero donor log2FCs against 0 (at least 3 donors); BH is applied across genes. avg_log2FC is the mean over all paired donors and r_rb the matched-pairs rank-biserial correlation.
The cell unit compares the cells directly with a Wilcoxon rank-sum test. A gene is significant when its BH-adjusted p < 'fdr' and |avg_log2FC| >= 'lfc'. Which groups to compare, and the thresholds, depend on the question and the data; the defaults are those of the SEA-AD analysis.
Value
A named list with one data frame per contrast ('NULL' when a contrast cannot be tested), with the attribute 'info' giving the cell and donor counts.
Examples
d <- cepip_toy
res <- cepip(d$x, d$y, donor = d$donor, id = d$id, verbose = FALSE)
expr <- t(as.matrix(d[, c("x", paste0("P", 1:8))]))
colnames(expr) <- d$id
deg <- cepip_deg(res, expr, contrast = "PostIP_vs_PreIP")
deg$PostIP_vs_PreIP
Per-donor share of cells in each group
Description
Counts each donor's cells per group and divides by the donor's total cell count. The total should include every cell of the donor in the population, also cells left out of the analysis (for example cells with GOI = 0), so 'n_total' can be given separately. Donors in 'n_total' without any cell in a group get 0.
Usage
cepip_donor(x, n_total = NULL)
Arguments
x |
A 'cepip', 'cepip_cmrf' or 'cepip_ip' object. |
n_total |
Optional named vector of total cells per donor. By default the number of analysed cells per donor. |
Value
A data frame with one row per donor: total cells, counts and percentages of PostIP, PreIP, NMPost and NMPre, and the Module percentage.
Examples
d <- cepip_toy
res <- cepip(d$x, d$y, donor = d$donor, id = d$id, verbose = FALSE)
head(cepip_donor(res))
Cell explanatory power and TREP cells
Description
Computes the explanatory power of each cell, EP_i = 1 - (y_i - f(x_i))^2 / (y_i - mean(y))^2, where f is the population-level curve (donor term excluded) and the null is the plain mean of y. The top round(n x DE) cells by EP are the top-ranked explanatory power (TREP) cells. DE here is the weighted population-level DE, 1 - RSS_w / TSS_w, from the same curve; without donors it equals the usual deviance explained.
Usage
cepip_ep(fit)
Arguments
fit |
A 'cepip_fit' object. |
Value
A data frame with one row per cell: id, x, y, donor, fitted value, EP, EP rank and TREP flag. The attribute 'de_pop' holds the DE used for the TREP count.
Examples
d <- cepip_toy
f <- cepip_fit(d$x, d$y, donor = d$donor, id = d$id)
ep <- cepip_ep(f)
head(ep)
sum(ep$is_TREP)
Download the SEA-AD example data
Description
Downloads the example data used in the paper from Hugging Face once and keeps them in a local cache. The data come from the Seattle Alzheimer's Disease Brain Cell Atlas (SEA-AD) and are subject to the Allen Institute Terms of Use (noncommercial research use, with citation); see <https://alleninstitute.org/legal/terms-of-use> and <https://alleninstitute.org/legal/citation-policy>.
Usage
cepip_example(
which = c("pvm3", "hkg8", "donors", "deg_NeoC", "deg_NonNeoC"),
cache = tools::R_user_dir("cepip", "cache"),
repo = "kahkengwong/cepip_example",
quiet = FALSE
)
Arguments
which |
Which file to load. |
cache |
Cache folder. |
repo |
Hugging Face dataset repository. |
quiet |
Suppress the download progress. |
Details
Files: * '"pvm3"': the 24,463 Micro-PVM_3-SEAAD cells with MT-CO3 > 0 (NeoC 13,331, NonNeoC 11,132): barcode, region set, region, donor, MT-CO3, the eight mitochondrial partner genes, the composite of the data export ('composite_part8'), ACTB and the eight other housekeeping genes, all as ln(CP10K + 1). * '"hkg8"': the 2,511 cells of the housekeeping-gene control arm of the paper (ACTB > 0 cells of the M4 arm): ACTB, the eight housekeeping genes and their export composite. * '"donors"': donor metadata (covariates and outcomes) and the total number of Micro-PVM_3-SEAAD cells per donor and region set. * '"deg_NeoC"', '"deg_NonNeoC"': sparse ln(CP10K + 1) matrices (genes x cells) of the PVM_3 cells for [cepip_deg()].
The files are kept in 'cache' (by default the user cache folder of R, see [tools::R_user_dir()]), so they are downloaded only once. [cepip_example_clear()] deletes them. When the download fails (for example without an internet connection), the function gives a message and returns 'NULL'.
Value
A data frame (or a sparse matrix for the DEG files), or 'NULL' when the download fails.
Examples
hk <- cepip_example("hkg8", cache = file.path(tempdir(), "cepip_cache"))
if (!is.null(hk)) dim(hk)
Delete the cached example data
Description
Deletes the files that [cepip_example()] saved in the cache folder, and the folder itself when it is then empty.
Usage
cepip_example_clear(cache = tools::R_user_dir("cepip", "cache"))
Arguments
cache |
Cache folder (the same as in [cepip_example()]). |
Value
The paths of the deleted files, invisibly.
Examples
cepip_example_clear(cache = file.path(tempdir(), "cepip_cache"))
Fit the CEP-IP GAM
Description
Fits 'y ~ x + s(x, bs = "tp", k)' by REML with 'select = TRUE' and 'gamma = 1.5', as in the published CEP-IP. When 'donor' is given, a donor random intercept 's(donor, bs = "re")' is added (if there are at least 'min_donors' donors) and each cell is weighted by 1 / (cells of its donor), so that every donor carries the same total weight.
Usage
cepip_fit(
x,
y,
donor = NULL,
id = NULL,
k = "prss",
k_max = 10L,
gamma = 1.5,
select = TRUE,
min_donors = 3L,
workers = 1L,
verbose = FALSE
)
Arguments
x |
Expression of the gene of interest (GOI), one value per cell. |
y |
Composite expression of the partner genes (see [cepip_composite()]). |
donor |
Optional donor label per cell. |
id |
Optional cell identifiers (for example barcodes). |
k |
'"prss"' for the PRSS search, or a fixed integer k. |
k_max |
Largest k in the PRSS search. |
gamma, select |
Passed to [mgcv::gam()]. |
min_donors |
Minimum number of donors for the donor random intercept. |
workers |
Number of parallel workers for the k search. |
verbose |
Print progress messages. |
Details
The basis dimension k is chosen by the lowest penalized residual sum of squares (PRSS) over k = 3 to min('k_max', unique x - 1), with ties going to the smaller k. The published pipeline ran a second phase of up to five iterations as a sanity check of this search; it cannot change the selected k because the first phase already fits every k in the range and each fit is deterministic, so it is not repeated here.
The package does not transform the data. The published CEP-IP analyses used log2(value + 1) of Seurat log-normalized data; the SEA-AD analysis used ln(CP10K + 1) directly.
Value
An object of class 'cepip_fit': the model, the selected k, the PRSS table, the deviance explained (DE) and the fitting frame.
Examples
d <- cepip_toy
f <- cepip_fit(d$x, d$y, donor = d$donor, id = d$id)
f$k
f$de
Cell groups for downstream analyses
Description
Returns the four groups used by [cepip_donor()] and [cepip_deg()]: Module Post-IP ('PostIP'), Module Pre-IP ('PreIP'), non-Module Post-IP ('NMPost') and non-Module Pre-IP ('NMPre'). With CMRF, the Module is OriRef plus Recovered; without CMRF (the published CEP-IP), the TREP cells take the place of the Module. The IP side follows the IP axis of the analysis (by default y < IPy is Pre-IP).
Usage
cepip_groups(x)
Arguments
x |
A 'cepip' object, a 'cepip_cmrf' object, or a 'cepip_ip' object. |
Value
A data frame with the cell id, donor, module flag, IP side and group.
Examples
d <- cepip_toy
res <- cepip(d$x, d$y, donor = d$donor, id = d$id, verbose = FALSE)
table(cepip_groups(res)$group)
Inflection point (IP) and IP reliability score (IPRS)
Description
Finds the CEP-IP inflection point.
Usage
cepip_ip(
x,
y,
donor = NULL,
id = NULL,
ip_k = "prss",
resid_k = 20L,
grid = 2000L,
z = 1.96,
ip_axis = c("y", "x"),
min_cells = 10L,
gamma = 1.5,
select = TRUE,
min_donors = 3L,
fit = NULL,
verbose = FALSE
)
Arguments
x |
Expression of the gene of interest (GOI), one value per cell. |
y |
Composite expression of the partner genes (see [cepip_composite()]). |
donor |
Optional donor label per cell. |
id |
Optional cell identifiers (for example barcodes). |
ip_k |
k of the GAM f: '"prss"' (default) or a fixed integer. The published SEA-AD analysis used 10. |
resid_k |
Largest k of the residual GAM R(x). |
grid |
Number of grid points for R(x). |
z |
Multiplier of the standard error for the R(x) bands. |
ip_axis |
'"y"' or '"x"', the axis that splits Pre-IP from Post-IP cells. |
min_cells |
Minimum number of top-EP cells for the residual GAM. |
gamma, select |
Passed to [mgcv::gam()]. |
min_donors |
Minimum number of donors for the donor random intercept. |
fit |
Optional 'cepip_fit' of the same cells, to reuse instead of refitting. |
verbose |
Print progress messages. |
Details
1. The CEP-IP GAM f is fitted on all cells ('ip_k'), and the top round(n x DE) cells by explanatory power are kept (see [cepip_ep()]). 2. The residuals y - f(x) of these cells are modelled by a residual GAM R(x) = s(x, bs = "tp", k = min('resid_k', unique x - 1)), fitted by REML (with the donor random intercept and inverse-donor weights when 'donor' is given; the donor term is excluded when predicting). 3. R(x) is evaluated on a grid of 'grid' points. IPx is the first negative-to-positive zero crossing, located by linear interpolation between the two grid points. 4. The 95 R(x) - 1.96 se(x) (upper), found by Brent's method ([stats::uniroot()]) from the band sign change nearest the IP. When a band has no sign change, the grid point where the band is closest to zero on that side is used. 5. IPy = f(IPx), the population-level curve at IPx. 6. IPRS = (C1 + C2 + C3 + C4 + C5) / 5: CI precision, a one-sided test of post-IP against pre-IP residuals, the slope of R(x) at the IP, the ribbon width of R(x), and the smoothness of R(x) combined with the share of post-IP residuals that are positive.
Cells are labelled Pre-IP or Post-IP on the axis given by 'ip_axis'. The default '"y"' splits on the response, y < IPy. The published CEP-IP analyses split on x (x < IPx), which 'ip_axis = "x"' reproduces.
Value
An object of class 'cepip_ip'. 'has_ip' is 'FALSE' when R(x) has no negative-to-positive crossing; the Pre-IP and Post-IP labels are then 'NA'.
Examples
d <- cepip_toy
ip <- cepip_ip(d$x, d$y, donor = d$donor, id = d$id)
ip$ip_x
ip$iprs
Save the tables and plots of one arm
Description
Writes all tables and plots of one arm into one folder, with no subfolders. Several region sets (or other subsets) of the same arm go into the same folder, each with its name in the file and sheet names.
Usage
cepip_save(
results,
dir,
arm = basename(dir),
extras = NULL,
plots = TRUE,
formats = c("png", "pdf"),
log = NULL,
max_cell_rows = 20000L
)
Arguments
results |
A 'cepip' object, or a named list of 'cepip' objects (one per region set, for example 'list(NeoC = a, NonNeoC = b)'). |
dir |
Output folder (created if needed). |
arm |
Arm name used in the workbook names. Defaults to the folder name. |
extras |
Optional named list (same names as 'results') of lists with any of 'screen' ([cepip_screen()]), 'donor' ([cepip_donor()]), 'deg' ([cepip_deg()]), 'associate' ([cepip_associate()]) and 'outcomes' (a donor table for [plot_donor()]) for that region set. |
plots |
Write the plots. |
formats |
Plot file formats. |
log |
Optional character vector written to '00_log.txt'. |
max_cell_rows |
Row limit for keeping cell-level sheets in the main workbook. |
Details
* '00_settings.txt': settings, package versions and the summary of each result. '00_log.txt': the text passed in 'log' (for example the console output of the run, captured with [utils::capture.output()] or [sink()]). * '<arm>_cepip_tables.xlsx': one sheet per table, numbered in pipeline order. When the cell-level tables are large (more than 'max_cell_rows' rows in total), they go into a second workbook, '<arm>_cepip_cells.xlsx'. * Plots as .png and .pdf, numbered in pipeline order.
Value
The paths of the files written, invisibly.
Examples
d <- cepip_toy
res <- cepip(d$x, d$y, donor = d$donor, id = d$id, verbose = FALSE)
out <- file.path(tempdir(), "toy")
files <- cepip_save(res, out, plots = FALSE)
basename(files)
unlink(out, recursive = TRUE)
Spearman-Kendall dual-filter screen for partner genes
Description
Screens candidate genes of interest (GOIs) for partner genes (dual-filter genes, DFGs) across the cells of one population, as in the published CEP-IP and the SEA-AD analysis.
Usage
cepip_screen(
expr,
goi = NULL,
donor = NULL,
strata = NULL,
min_pct = 0.25,
min_mean = 0.3,
min_partner_cells = 10L,
min_donor_cells = 5L,
rs = 0.6,
tau = 0.5,
negative = TRUE,
max_cells = NULL,
seed = 42L,
goi_block = 500L,
partner_block = 2000L,
verbose = FALSE
)
Arguments
expr |
Expression matrix, genes x cells (dense or sparse), with gene names as row names. |
goi |
Candidate GOIs. 'NULL' screens every gene. |
donor |
Optional donor label per cell. |
strata |
Optional stratum per cell (for example region) for 'max_cells'. |
min_pct, min_mean |
GOI expression gate. |
min_partner_cells |
Minimum number of expressing cells of a partner. |
min_donor_cells |
Minimum number of cells per donor. |
rs, tau |
Spearman and Kendall thresholds. |
negative |
Keep negative DFGs. |
max_cells |
Optional number of cells to screen. |
seed |
Seed of the subsample. |
goi_block, partner_block |
Block sizes for the Spearman step (memory). |
verbose |
Print progress messages. |
Details
1. Donors with fewer than 'min_donor_cells' cells are dropped. 2. A GOI is screened if it is expressed (value > 0) in at least 'min_pct' of the cells and its mean over the expressing cells is at least 'min_mean'. In the SEA-AD data (ln(CP10K + 1)) these defaults kept 2,500 to 4,000 genes per population. Adjust them for other data and normalizations. 3. Partner genes must be expressed in at least 'min_partner_cells' cells. 4. Spearman rs is computed for every GOI-partner pair (average ranks, zeros included). Kendall tau-b is computed for the pairs with |rs| >= 'rs'. 5. A partner is a DFG when rs >= 'rs' and tau-b >= 'tau' (positive), or, with 'negative = TRUE', rs <= -'rs' and tau-b <= -'tau' (negative).
GOIs are ranked by the number of DFGs, then by the largest rs. The choice of the GOI and of the DFG panel for [cepip_composite()] is left to the user.
Run time grows with cells x GOIs x partners. For large populations, 'max_cells' screens a random subsample (stratified by 'strata' when given). Kendall tau-b uses 'pcaPP::cor.fk()' when pcaPP is installed (fast, same result), else [stats::cor()].
Value
An object of class 'cepip_screen' with 'ranking' (one row per screened GOI with at least one DFG), 'partners' (one row per GOI-DFG pair), 'gate' (the GOI gate of every candidate) and 'cells' (the cells screened).
Examples
d <- cepip_toy
expr <- t(as.matrix(d[, c("x", paste0("P", 1:8))]))
colnames(expr) <- d$id
sc <- cepip_screen(expr, goi = "x", donor = d$donor)
sc$ranking
head(sc$partners)
Spearman-Kendall cell-level scores (SKCS)
Description
Splits the global Spearman and Kendall correlations between x and y into per-cell contributions.
Usage
cepip_skcs(
x,
y,
donor = NULL,
id = NULL,
s_min = 0.05,
k_min = 0.5,
fit = TRUE,
verbose = FALSE,
...
)
Arguments
x |
Expression of the gene of interest (GOI), one value per cell. |
y |
Composite expression of the partner genes (see [cepip_composite()]). |
donor |
Optional donor label per cell. |
id |
Optional cell identifiers (for example barcodes). |
s_min, k_min |
SKCS-core gates. |
fit |
Refit the GAM on the SKCS-core cells and define TREP cells there. |
verbose |
Print progress messages. |
... |
Passed to [cepip_fit()] for the SKCS-core fit (for example 'k'). |
Details
* S_i = (r_x,i - mean r_x)(r_y,i - mean r_y) / (n - 1), from average ranks; S_i,norm = S_i / ((n - 1) / 4). * K_i = (C_i - D_i) / sqrt((n - 1 - t_x,i)(n - 1 - t_y,i)), where C_i and D_i are the concordant and discordant pairs of cell i and t_x,i and t_y,i its ties in x and in y (tau-b form).
SKCS-core cells pass both gates, S_i,norm > 's_min' and K_i > 'k_min'. The CEP-IP GAM is then refitted on the SKCS-core cells (with the same settings as [cepip_fit()]) and the TREP cells are defined within the core.
The gates S_i,norm > 0.05 and K_i > 0.5 were introduced in the SEA-AD analysis; the published CEP-IP did not use cell-level scores.
Value
An object of class 'cepip_skcs': a per-cell table ('cells') with S_i, S_i,norm, K_i and the SKCS-core flag, the global rho and tau, and, when 'fit = TRUE', the SKCS-core fit and its EP and TREP table.
Examples
d <- cepip_toy
s <- cepip_skcs(d$x, d$y, donor = d$donor, id = d$id)
s$n_core
Simulated example data
Description
1,200 simulated cells from 12 donors with a gene of interest 'x', eight partner genes 'P1' to 'P8' and their composite 'y' (the row mean of the partners). The data are simulated (see 'data-raw/make_toy.R') and run through the whole method in under a minute.
Usage
cepip_toy
Format
A data frame with 1,200 rows and columns 'id', 'donor', 'region', 'x', 'P1' to 'P8' and 'y'.
Validate the CMRF Module by refitting the CEP-IP GAM (Validation1)
Description
Refits the CEP-IP GAM (same settings as [cepip_fit()], k by PRSS) on the Module cells (OriRef plus Recovered), recomputes EP and the TREP cells, and labels every Module cell Pre-IP or Post-IP on the IP axis. This gives the Recovered cells their own four-group labels. The agreement between the refitted and the original four-group labels of the OriRef cells is reported.
Usage
cepip_validate(cmrf, min_donors = 3L, workers = 1L, verbose = FALSE)
Arguments
cmrf |
A 'cepip_cmrf' object. |
min_donors |
Minimum number of donors for the donor random intercept. |
workers |
Number of parallel workers for the k search. |
verbose |
Print progress messages. |
Details
The refitted labels are kept in 'group_ref'. Plots of the CMRF classes use the original labels ('group_orig'), because the Original-only cells have no refitted label.
Value
An object of class 'cepip_validate': the refit, the per-cell table of the Module cells, the label agreement and the cross-tabulation.
Examples
d <- cepip_toy
# CMRF with a short gate of 5 repeats (the default is 40)
res <- cepip(d$x, d$y, donor = d$donor, id = d$id, cmrf = TRUE, repeats = 5,
tests = FALSE, verbose = FALSE)
v <- cepip_validate(res$cmrf)
v$de
v$agreement
Forest plot of donor-level associations
Description
The standardized beta of the predictor with its 95 [cepip_associate()], in the order given. Point size shows delta R2; filled points have FDR < 0.05. delta R2 and the FDR are written beside each row.
Usage
plot_associate(assoc, title = NULL)
Arguments
assoc |
Output of [cepip_associate()]. |
title |
Optional plot title. |
Value
A ggplot object.
Examples
d <- cepip_toy
res <- cepip(d$x, d$y, donor = d$donor, id = d$id, verbose = FALSE)
dt <- cepip_donor(res)
set.seed(1)
dt$score <- 0.05 * dt$pct_PostIP + rnorm(nrow(dt))
a <- cepip_associate(dt, "pct_PostIP", outcomes = "score", min_donors = 10)
plot_associate(a)
Plot the formation of the CMRF band (Figure 6)
Description
Panels: (1) the spine fitted on the SKCS-core TREP cells, (2) all SKCS-core cells in their four groups, (3) the flat delta_base band, (4) the adaptive delta_adapt band over delta_base, a QC panel zoomed to IP_u +- 0.25 (QC1 the width at IP_u, QC2 the widths at IP_u +- bw, where 60.65 percent of the widening is left, and QC3 the tangents of the band edges at IP_u +- bw), (5) the Stage 1 candidates, (6) the SKCS-core cells outside the band (Original only) and, when 'cmrf' is given, (7) the cells recovered by the gate and those it removed and (8) the final Module. This plot is useful to choose 'kappa' and 'bw' with [cepip_band()] before running [cepip_cmrf()].
Usage
plot_band(band, cmrf = NULL, title = NULL)
Arguments
band |
A 'cepip_band' object, or a 'cepip_cmrf' or 'cepip' object. |
cmrf |
Optional 'cepip_cmrf' object; adds panels 7 and 8. |
title |
Optional overall title. |
Value
A gtable of class 'cepip_figure' (printing it draws it; it can also be drawn with [grid::grid.draw()]). The attribute 'cepip_size' gives a suitable width and height in inches.
Examples
d <- cepip_toy
# CMRF with a short gate of 5 repeats (the default is 40)
res <- cepip(d$x, d$y, donor = d$donor, id = d$id, cmrf = TRUE, repeats = 5,
tests = FALSE, verbose = FALSE)
plot_band(res)
Plot the cells and the CMRF band in copula space
Description
(u, v) of all cells, the SKCS-core cells coloured by their four groups, the spine g(u), the flat delta_base band, the adaptive delta_adapt band and the IP_u and IP_v lines.
Usage
plot_copula(band, title = NULL)
Arguments
band |
A 'cepip_band' object, or a 'cepip_cmrf' or 'cepip' object. |
title |
Optional plot title. |
Value
A ggplot object.
Examples
d <- cepip_toy
# CMRF with a short gate of 5 repeats (the default is 40)
res <- cepip(d$x, d$y, donor = d$donor, id = d$id, cmrf = TRUE, repeats = 5,
tests = FALSE, verbose = FALSE)
plot_copula(res)
Volcano plot of one differential expression contrast
Description
avg_log2FC (x) against -log10 of the BH-adjusted p (y) for one contrast of [cepip_deg()]. Significant genes are coloured by direction and the top 10 are labelled; dashed lines mark the 'lfc' and 'fdr' thresholds.
Usage
plot_deg(deg, contrast = NULL, title = NULL)
Arguments
deg |
A 'cepip_deg' object. |
contrast |
Name of the contrast. By default the first one that was tested. |
title |
Optional plot title. |
Value
A ggplot object.
Examples
d <- cepip_toy
res <- cepip(d$x, d$y, donor = d$donor, id = d$id, verbose = FALSE)
expr <- t(as.matrix(d[, c("x", paste0("P", 1:8))]))
colnames(expr) <- d$id
deg <- cepip_deg(res, expr, contrast = "PostIP_vs_PreIP")
plot_deg(deg)
Plot the per-donor Module Post-IP percentage
Description
With only the table of [cepip_donor()], a sorted bar chart of 'pct_PostIP' per donor. With 'data' and 'outcome', a scatter of 'pct_PostIP' against the outcome with a linear fit.
Usage
plot_donor(donor_table, outcome = NULL, data = NULL, title = NULL)
Arguments
donor_table |
Output of [cepip_donor()], or a 'cepip', 'cepip_cmrf' or 'cepip_ip' object (then [cepip_donor()] is run on it). |
outcome |
Optional column name of 'data'. |
data |
Optional donor table with a 'donor' column and the outcome. |
title |
Optional plot title. |
Value
A ggplot object.
Examples
d <- cepip_toy
res <- cepip(d$x, d$y, donor = d$donor, id = d$id, verbose = FALSE)
plot_donor(res)
Plot the CEP-IP GAM with the inflection point
Description
GOI (x) against the composite (y). Cells are coloured by the four groups (Pre-IP or Post-IP, TREP or non-TREP) when an IP was found; otherwise TREP cells are coloured and non-TREP cells are gray. The population curve (donor term excluded) is drawn with its 95 as a shaded strip, and IPy as a horizontal line.
Usage
plot_fit(x, title = NULL)
Arguments
x |
A 'cepip' object (uses 'x$fit', 'x$ep' and 'x$ip'), a 'cepip_ip' object, or a 'cepip_fit' object (no IP lines). |
title |
Optional plot title. |
Value
A ggplot object.
Examples
d <- cepip_toy
res <- cepip(d$x, d$y, donor = d$donor, id = d$id, verbose = FALSE)
plot_fit(res)
Pitcher plot of the CMRF classes in copula space (Figure 8A)
Description
u_SK = F(S_i,norm) (x) against v_SK = F(K_i) (y). Module cells are the OriRef cells in their four groups and the Recovered cells (purple diamonds); non-Module cells are the Excluded cells (gray) and the Original only cells (tan triangles). The S_i,norm and K_i gates are drawn in copula coordinates.
Usage
plot_module(cmrf, title = NULL)
Arguments
cmrf |
A 'cepip_cmrf' object, or a 'cepip' object. |
title |
Optional plot title. |
Value
A ggplot object.
Examples
d <- cepip_toy
# CMRF with a short gate of 5 repeats (the default is 40)
res <- cepip(d$x, d$y, donor = d$donor, id = d$id, cmrf = TRUE, repeats = 5,
tests = FALSE, verbose = FALSE)
plot_module(res)
Pitcher plot of the SKCS scores
Description
K_i (y) against S_i,norm (x) for every cell, with the S_i,norm and K_i gates. SKCS-core cells are coloured by the four groups when 'ip' is given, otherwise in one colour.
Usage
plot_pitcher(skcs, ip = NULL, title = NULL)
Arguments
skcs |
A 'cepip_skcs' object, or a 'cepip' object (its 'skcs' is used, and its band or IP gives the groups). |
ip |
Optional 'cepip_band', 'cepip_cmrf' or 'cepip_ip' object of the same cells, for the four-group colours. |
title |
Optional plot title. |
Value
A ggplot object.
Examples
d <- cepip_toy
res <- cepip(d$x, d$y, donor = d$donor, id = d$id, skcs = TRUE, verbose = FALSE)
plot_pitcher(res)
Plot the PRSS search over the basis dimension k
Description
Plot the PRSS search over the basis dimension k
Usage
plot_prss(fit, title = NULL)
Arguments
fit |
A 'cepip_fit' object, or a 'cepip' or 'cepip_ip' object (its fit is used). |
title |
Optional plot title. |
Value
A ggplot object.
Examples
d <- cepip_toy
res <- cepip(d$x, d$y, donor = d$donor, id = d$id, verbose = FALSE)
plot_prss(res)
Plot the residual GAM R(x) and the IP
Description
R(x) with its 95 (faint points), the zero line and the IPx with its CI. The subtitle gives the IPRS, its tier and the components C1 to C5.
Usage
plot_residual(ip, title = NULL)
Arguments
ip |
A 'cepip_ip' object, or a 'cepip' object (its 'ip' is used). |
title |
Optional plot title. |
Value
A ggplot object.
Examples
d <- cepip_toy
res <- cepip(d$x, d$y, donor = d$donor, id = d$id, verbose = FALSE)
plot_residual(res)
Plot the RSk-fCV gate of CMRF Stage 2 (Figure 7)
Description
Panels for the first or the last gate iteration: (A) the fold layout of repeat 1, (B) the u strata of the pool, (C) curve A (isotonic fit on the OriRef cells) and curve B (on the pool) with the candidates' residual sticks, (D) the per-fold RMSE of repeat 1, (E) the RMSE_A and RMSE_B heatmaps (repeat x fold), (F) the RMSE distributions with the Wilcoxon p and the paired differences against BV, (G) the candidates ranked by |v - PAVA(u)|, (I) the gate convergence over all iterations and (Q) the class counts by IP side. When the class tests were run, (P) |v - g(u)| by class with the Jonckheere-Terpstra p and (O) the Recovered cells with the Hoeffding D and distance correlation p are added.
Usage
plot_rskfcv(cmrf, iteration = c("first", "last"), title = NULL)
Arguments
cmrf |
A 'cepip_cmrf' object, or a 'cepip' object. |
iteration |
'"first"' or '"last"'. |
title |
Optional overall title. |
Details
The pool of the first iteration is the OriRef cells plus all candidates; the pool of the last iteration is the OriRef cells plus the candidates still active (the Recovered cells).
Value
A gtable of class 'cepip_figure' (printing it draws it; it can also be drawn with [grid::grid.draw()]). The attribute 'cepip_size' gives a suitable width and height in inches.
Examples
d <- cepip_toy
# CMRF with a short gate of 5 repeats (the default is 40)
res <- cepip(d$x, d$y, donor = d$donor, id = d$id, cmrf = TRUE, repeats = 5,
tests = FALSE, verbose = FALSE)
plot_rskfcv(res, "first")
plot_rskfcv(res, "last")
Plot the partner-gene screen of one GOI
Description
Scatter of Spearman rs (x) against Kendall tau-b (y) for the partner genes of one GOI that passed the Spearman step of [cepip_screen()]. The dual-filter genes (DFGs) are highlighted and labelled; dashed lines mark the rs and tau thresholds.
Usage
plot_screen(screen, goi = NULL, title = NULL)
Arguments
screen |
A 'cepip_screen' object. |
goi |
The GOI to show. By default the top-ranked GOI. |
title |
Optional plot title. |
Value
A ggplot object.
Examples
d <- cepip_toy
expr <- t(as.matrix(d[, c("x", paste0("P", 1:8))]))
colnames(expr) <- d$id
sc <- cepip_screen(expr, goi = "x", donor = d$donor)
plot_screen(sc)
Plot the Validation1 refit of the CMRF Module
Description
Left: the Module cells (x against y) in their refitted four groups (Recovered cells as diamonds) with the refitted population curve and its 95 OriRef cells, with the agreement in the title.
Usage
plot_validate(validate, cmrf = NULL, title = NULL)
Arguments
validate |
A 'cepip_validate' object, or a 'cepip' object (its 'validate' and 'cmrf' are used). |
cmrf |
Optional 'cepip_cmrf' object, for the IPx and IPy lines. |
title |
Optional overall title. |
Value
A gtable of class 'cepip_figure' (printing it draws it; it can also be drawn with [grid::grid.draw()]). The attribute 'cepip_size' gives a suitable width and height in inches.
Examples
d <- cepip_toy
# CMRF with a short gate of 5 repeats (the default is 40)
res <- cepip(d$x, d$y, donor = d$donor, id = d$id, cmrf = TRUE, repeats = 5,
tests = FALSE, verbose = FALSE)
plot_validate(res)
Draw a multi-panel cepip figure
Description
The multi-panel plots ([plot_band()], [plot_rskfcv()], [plot_validate()]) return a gtable with the extra class 'cepip_figure', so that printing it (at the console or in a report) draws it on a new page.
Usage
## S3 method for class 'cepip_figure'
print(x, ...)
Arguments
x |
A 'cepip_figure' object. |
... |
Not used. |
Value
'x', invisibly.