IntrinsicScan: label-blind feature-wise binary partitioning
Introduction
IntrinsicScan asks a deliberately narrow question for
each continuous feature: can the feature define its own binary partition
without using phenotype labels, and how strongly does that fixed
partition subsequently concord with a binary phenotype?
The partition is the exact optimum of a one-dimensional two-medoid L1 objective, where L1 loss is the sum of absolute distances between observations and the medoid of their assigned state. Phenotype labels are not available to the partitioning step. Labels are used only afterward to orient the low/high states and calculate concordance and inference. A high concordance score is not evidence by itself for biological bimodality.
The package is designed for high-throughput continuous molecular
measurements such as gene-expression data and interoperates with
SummarizedExperiment and ExpressionSet.
Related label-blind discretization approaches include StepMiner, which
fits step functions using a squared-error objective.
IntrinsicScan instead uses an L1/medoid objective and
explicitly separates intrinsic boundary discovery from phenotype
concordance, inferential evidence, and resampling stability.
Installation
After release through Bioconductor, IntrinsicScan can be
installed with:
How the intrinsic partition is defined
For each molecular feature, IntrinsicScan identifies two
intrinsic states without using phenotype labels to determine the
partition.
Let the observed values for one feature, after sorting, be
\[ x_{(1)} \leq x_{(2)} \leq \cdots \leq x_{(n)}. \]
IntrinsicScan considers every admissible split between
two distinct adjacent values. For a candidate split at position \(k\), the observations are divided into a
lower state and a higher state, and the within-state L1 loss is
\[ L(k) = \sum_{i \leq k} |x_{(i)} - m_L(k)| + \sum_{i > k} |x_{(i)} - m_H(k)|, \]
where \(m_L(k)\) and \(m_H(k)\) are the medoids of the lower and higher states, respectively.
The L1 loss is the sum of absolute distances between observations and their state medoid. Unlike a squared-error objective, deviations are therefore not squared.
The selected partition is the split that minimizes this total L1 loss:
\[ k^* = \arg\min_k L(k). \]
The threshold is placed halfway between the two adjacent observed values that define the selected split. If several candidate splits have the same optimum within numerical tolerance, ties are resolved deterministically by first preferring the largest gap between the two states and, if necessary, the split closest to the middle of the ordered observations.
Importantly, phenotype labels are not used during this optimization. The low and high states are therefore determined entirely from the observed values of the feature. Only after the partition has been fixed are phenotype labels used to assess concordance.
For a fixed partition, IntrinsicScan evaluates phenotype
concordance using the orientation-free Youden index
\[ J = |\mathrm{TPR}_{high} - \mathrm{FPR}_{high}|. \]
After orienting the two states so that the state more associated with the phenotype-positive class is treated as positive, this is equivalent to
\[ J = \mathrm{sensitivity} + \mathrm{specificity} - 1. \]
Thus, \(J = 0\) indicates no concordance between the intrinsic partition and the phenotype, whereas \(J = 1\) indicates perfect concordance in the analyzed samples.
Statistical significance is assessed only after the partition is
fixed. Conditional on the numbers of phenotype-positive samples,
phenotype-negative samples and samples assigned to the high state, the
number of phenotype-positive samples in the high state follows a
hypergeometric distribution under exchangeable phenotype labels.
exact_J_p is the exact probability of obtaining an
orientation-free \(J\) at least as
large as the observed value under this conditional null model.
exact_J_fdr is the Benjamini-Hochberg-adjusted value across
the analyzed features.
Features are ranked primarily by decreasing \(J\), with increasing
exact_J_fdr and decreasing relative_gap used
as successive tie breakers. These quantities remain separate result axes
rather than being combined into a composite score.
Minimal example
library(IntrinsicScan)
x <- rbind(
feature_A = c(1, 1.1, 1.2, 1.3, 1.4, 5, 5.1, 5.2, 5.3, 5.4),
feature_B = c(1, 2, 3, 2, 3, 2, 3, 4, 3, 4),
feature_C = rep(7, 10)
)
colnames(x) <- paste0("S", seq_len(ncol(x)))
y <- rep(c("control", "case"), each = 5)
fit <- intrinsicScan(
x,
phenotype = y,
positive = "case",
verbose = FALSE
)
res <- intrinsicResults(fit)
res[, c(
"feature_id",
"valid",
"reason",
"threshold",
"low_medoid",
"high_medoid",
"J",
"positive_state",
"exact_J_p",
"exact_J_fdr",
"relative_gap"
)]
#> DataFrame with 3 rows and 11 columns
#> feature_id valid reason threshold low_medoid
#> <character> <logical> <character> <numeric> <numeric>
#> feature_A feature_A TRUE NA 3.2 1.2
#> feature_B feature_B TRUE NA 2.5 2.0
#> feature_C feature_C FALSE zero-variance feature NA NA
#> high_medoid J positive_state exact_J_p exact_J_fdr
#> <numeric> <numeric> <character> <numeric> <numeric>
#> feature_A 5.2 1.0 high 0.00793651 0.0238095
#> feature_B 3.0 0.4 high 0.52380952 0.7857143
#> feature_C NA 0.0 NA 1.00000000 1.0000000
#> relative_gap
#> <numeric>
#> feature_A 0.911392
#> feature_B 1.000000
#> feature_C NAThe returned table contains one row per analyzed feature. The columns shown above summarize the fitted intrinsic partition and its subsequent relationship with the phenotype.
feature_A illustrates a feature with a clear intrinsic
low/high partition that also aligns perfectly with the phenotype in this
small example. Its threshold is determined from the feature values
alone; the phenotype labels are used only afterward to quantify
concordance.
feature_B also admits an intrinsic two-state partition,
but the resulting states are less concordant with the phenotype. This
distinction is central to the method: IntrinsicScan does
not ask whether the phenotype can be used to create a good split, but
whether a split discovered independently of the phenotype subsequently
aligns with it.
feature_C is constant and therefore cannot define a
valid two-state partition. It is retained in the result with
valid = FALSE, and reason describes why the
partition could not be fitted.
For an initial interpretation, the most useful columns are:
threshold: the boundary separating the fitted low and high states;low_medoidandhigh_medoid: representative observed values of the two fitted states;J: the orientation-free Youden index describing phenotype concordance, ranging from 0 for no concordance to 1 for perfect concordance;positive_state: whether the phenotype-positive class is associated with the fitted high or low state;exact_J_pandexact_J_fdr: exact conditional statistical evidence for the observed concordance and its multiple-testing-adjusted value;relative_gap: the distance between the two observations surrounding the fitted threshold, scaled by the feature interquartile range.
The complete set of returned variables, including optional
resampling-stability summaries, is documented in
?intrinsicResults.
Output data structure
intrinsicScan() returns a
SummarizedExperiment. The analyzed measurements are stored
in assay "expression", with rows reordered according to the
IntrinsicScan ranking. By default, fitted high/low state
assignments are also retained as the logical assay
"state".
SummarizedExperiment::assayNames(fit)
#> [1] "expression" "state"
rd <- SummarizedExperiment::rowData(fit)
rd[, c(
"scan_feature_id",
"scan_threshold",
"scan_J",
"scan_exact_J_fdr",
"scan_rank_J"
)]
#> DataFrame with 3 rows and 5 columns
#> scan_feature_id scan_threshold scan_J scan_exact_J_fdr scan_rank_J
#> <character> <numeric> <numeric> <numeric> <integer>
#> feature_A feature_A 3.2 1.0 0.0238095 1
#> feature_B feature_B 2.5 0.4 0.7857143 2
#> feature_C feature_C NA 0.0 1.0000000 3Feature-level results are appended to rowData() with the
prefix scan_. The same results can be obtained without this
prefix using intrinsicResults(). The original sample
metadata are retained in colData(), together with
IntrinsicScanLabel, which records the phenotype labels used
for the scan.
Additional analysis metadata, including the positive and negative
phenotype labels, analysis parameters, result table, and optional
resampling details, are stored in the IntrinsicScan entry
of metadata().
names(S4Vectors::metadata(fit)$IntrinsicScan)
#> [1] "engine" "method" "positive_label"
#> [4] "negative_label" "parameters" "analyzed_assay"
#> [7] "results" "resampling_details" "source_class"This organization keeps the fitted partitions, input measurements, sample annotations and feature-level statistics aligned in a standard Bioconductor container.
Bioconductor containers
For a SummarizedExperiment, phenotype may
either be supplied directly or name a column in colData().
The assay argument selects the assay to analyze.
For an ExpressionSet, phenotype may
similarly name a column in pData(), and the expression
matrix is obtained from exprs().
IntrinsicScan currently converts the selected
measurements to an in-memory numeric matrix before analysis. No
normalization, transformation, imputation, or feature filtering is
performed by the package, so the supplied measurements must already be
appropriately preprocessed and finite.
Biological example: AML versus non-leukemic bone marrow
The minimal example above is useful for understanding the mechanics
of the method but is intentionally artificial. A more realistic
application can be illustrated using the Bioconductor experiment-data
package leukemiasEset.
leukemiasEset contains RMA-normalized microarray
expression measurements from 60 bone-marrow samples representing acute
lymphoblastic leukemia (ALL), acute myeloid leukemia (AML), chronic
lymphocytic leukemia (CLL), chronic myeloid leukemia (CML), and
non-leukemic bone marrow (NoL), with 12 samples in each group. The data
are represented as an ExpressionSet with 20,172
Ensembl-level features.
Because IntrinsicScan operates on a binary phenotype,
here we restrict the dataset to AML and NoL samples. For vignette speed,
the example uses the first 1,000 features in their original dataset
order; this deterministic subset is not selected using the phenotype
labels. A research analysis would ordinarily scan all features that pass
the analyst’s prespecified preprocessing and quality-control
criteria.
The biological question in this example is therefore: which gene-expression features in this real dataset contain an intrinsic two-state expression boundary, discovered without using the diagnostic labels, that subsequently aligns with AML versus non-leukemic bone marrow?
library(leukemiasEset)
data("leukemiasEset")
keep <- leukemiasEset$LeukemiaType %in% c("AML", "NoL")
aml_nol <- leukemiasEset[, keep]
aml_nol <- aml_nol[seq_len(min(1000L, nrow(aml_nol))), ]
table(as.character(aml_nol$LeukemiaType))
#>
#> AML NoL
#> 12 12
bio_fit <- intrinsicScan(
aml_nol,
phenotype = "LeukemiaType",
positive = "AML",
nResamples = 0L,
keepStates = FALSE,
verbose = FALSE
)
bio_res <- intrinsicResults(bio_fit)
utils::head(
bio_res[, c(
"feature_id",
"threshold",
"J",
"positive_state",
"exact_J_p",
"exact_J_fdr",
"relative_gap"
)],
10
)
#> DataFrame with 10 rows and 7 columns
#> feature_id threshold J positive_state exact_J_p
#> <character> <numeric> <numeric> <character> <numeric>
#> ENSG00000065618 ENSG00000065618 5.61341 1.000000 low 7.39602e-07
#> ENSG00000008438 ENSG00000008438 9.42937 0.916667 low 9.61483e-06
#> ENSG00000012779 ENSG00000012779 7.12663 0.916667 low 9.61483e-06
#> ENSG00000047597 ENSG00000047597 6.19368 0.916667 low 9.61483e-06
#> ENSG00000055118 ENSG00000055118 6.82603 0.916667 low 9.61483e-06
#> ENSG00000024526 ENSG00000024526 4.66474 0.833333 low 6.73038e-05
#> ENSG00000010292 ENSG00000010292 6.66544 0.833333 low 6.73038e-05
#> ENSG00000064763 ENSG00000064763 8.01576 0.833333 low 6.73038e-05
#> ENSG00000002586 ENSG00000002586 9.80279 0.833333 high 6.73038e-05
#> ENSG00000057704 ENSG00000057704 4.61738 0.833333 low 6.73038e-05
#> exact_J_fdr relative_gap
#> <numeric> <numeric>
#> ENSG00000065618 0.000739602 0.432746
#> ENSG00000008438 0.001922966 0.323340
#> ENSG00000012779 0.001922966 0.254630
#> ENSG00000047597 0.001922966 0.192030
#> ENSG00000055118 0.001922966 0.137847
#> ENSG00000024526 0.004206488 0.466115
#> ENSG00000010292 0.004206488 0.314841
#> ENSG00000064763 0.004206488 0.296103
#> ENSG00000002586 0.004206488 0.290987
#> ENSG00000057704 0.004206488 0.267944The first rows correspond to features whose independently fitted intrinsic partitions show the strongest concordance with AML versus NoL according to the package ranking. A high-ranking feature therefore has two distinct properties: its low/high boundary was obtained without using the diagnostic labels, and the resulting fixed partition subsequently showed strong phenotype concordance.
The fitted threshold of the top-ranked valid feature can also be visualized against the two phenotype groups:
top_id <- as.character(bio_res$feature_id[which(bio_res$valid)[1]])
top_threshold <- bio_res$threshold[match(top_id, bio_res$feature_id)]
top_values <- Biobase::exprs(aml_nol)[top_id, ]
group <- factor(
as.character(aml_nol$LeukemiaType),
levels = c("NoL", "AML")
)
boxplot(
top_values ~ group,
xlab = "Phenotype",
ylab = "Expression",
main = paste("Top IntrinsicScan feature:", top_id)
)
abline(h = top_threshold, lty = 2)
The horizontal line is the threshold discovered from the feature values before the AML/NoL labels were used for concordance assessment.
This example demonstrates the workflow and interpretation of an
IntrinsicScan result; it is not intended as independent
validation of a leukemia biomarker. The same samples are used to
discover the intrinsic partition and to assess its phenotype
concordance. Independent datasets or appropriately designed validation
analyses are required for biological or clinical claims.
Stability
Stability assessment is optional and disabled by default in the public interface because it can dominate runtime for genome-wide scans.
When enabled, partitions are repeatedly fitted in resampled training
data and the resulting thresholds and phenotype orientations are
evaluated both in the training observations and, when available, in
held-out observations. IntrinsicScan reports threshold
variability, training and held-out \(J\), orientation stability, partition
stability, and the fractions of resamples showing exact concordance.
Each feature receives a deterministic random-number stream derived from the master seed and its feature identifier. Consequently, reranking or reordering other features does not change that feature’s resampling history. Random-number seeding is scoped locally so that package calls do not leave the user’s global RNG state altered.
Interpretation limits
IntrinsicScan is not a general clustering algorithm, a
supervised cutoff optimizer, or a formal test of bimodality. The
intrinsic partition and phenotype-concordance assessment answer distinct
questions and should be interpreted separately.
A strong \(J\) indicates that the fixed intrinsic partition aligns with the supplied binary phenotype in the analyzed samples. It does not establish that the feature has two biologically discrete states, that the association is causal, or that it will generalize to independent data.
The L1/medoid objective is designed to provide a conservative, gross-outlier-resistant intrinsic binary partition. This design can reduce sensitivity to genuinely rare expression states.
References
Sahoo D, Dill DL, Tibshirani R, Plevritis SK (2007). Extracting binary signals from microarray time-course data. Nucleic Acids Research 35:3705-3712. doi:10.1093/nar/gkm284.
Aibar S, Fontanillo C, De Las Rivas J. leukemiasEset:
Leukemia’s microarray gene expression data (ExpressionSet).
Bioconductor experiment-data package. doi:10.18129/B9.bioc.leukemiasEset.
Haferlach T, Kohlmann A, Wieczorek L, et al. (2010). Clinical utility of microarray-based gene expression profiling in the diagnosis and subclassification of leukemia: report from the International Microarray Innovations in Leukemia Study Group. Journal of Clinical Oncology 28:2529-2537. doi:10.1200/JCO.2009.23.4732.
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] leukemiasEset_1.49.0 Biobase_2.73.2 BiocGenerics_0.59.12
#> [4] generics_0.1.4 IntrinsicScan_0.99.4 BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] Matrix_1.7-6 jsonlite_2.0.0
#> [3] compiler_4.6.1 BiocManager_1.30.27
#> [5] SummarizedExperiment_1.43.0 GenomicRanges_1.65.4
#> [7] jquerylib_0.1.4 IRanges_2.47.5
#> [9] Seqinfo_1.3.2 yaml_2.3.12
#> [11] fastmap_1.2.0 lattice_0.23-1
#> [13] XVector_0.53.0 R6_2.6.1
#> [15] S4Arrays_1.13.2 knitr_1.52
#> [17] DelayedArray_0.39.8 MatrixGenerics_1.25.0
#> [19] maketools_1.3.2 bslib_0.12.0
#> [21] rlang_1.3.0 cachem_1.1.0
#> [23] xfun_0.61 sass_0.4.10
#> [25] sys_3.4.3 otel_0.2.0
#> [27] SparseArray_1.13.4 cli_3.6.6
#> [29] digest_0.6.39 grid_4.6.1
#> [31] lifecycle_1.0.5 S4Vectors_0.51.10
#> [33] evaluate_1.0.5 buildtools_1.0.0
#> [35] abind_1.4-8 stats4_4.6.1
#> [37] rmarkdown_2.32 matrixStats_1.5.0
#> [39] tools_4.6.1 htmltools_0.5.9