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.
After release through Bioconductor, IntrinsicScan can be installed with:
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("IntrinsicScan")
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.
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 NA
The 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_medoid and high_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_p and exact_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.
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 3
Feature-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.
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.
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.267944
The 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 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.
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.
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.
sessionInfo()
#> 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] 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] Rcpp_1.1.2 tinytex_0.61
#> [7] SummarizedExperiment_1.43.0 magick_2.9.1
#> [9] GenomicRanges_1.65.4 jquerylib_0.1.4
#> [11] IRanges_2.47.5 Seqinfo_1.3.2
#> [13] yaml_2.3.12 fastmap_1.2.0
#> [15] lattice_0.23-1 R6_2.6.1
#> [17] XVector_0.53.0 S4Arrays_1.13.2
#> [19] knitr_1.52 DelayedArray_0.39.8
#> [21] bookdown_0.48 MatrixGenerics_1.25.0
#> [23] bslib_0.12.0 rlang_1.3.0
#> [25] cachem_1.1.0 xfun_0.61
#> [27] sass_0.4.10 otel_0.2.0
#> [29] SparseArray_1.13.4 cli_3.6.6
#> [31] magrittr_2.0.5 digest_0.6.39
#> [33] grid_4.6.1 lifecycle_1.0.5
#> [35] S4Vectors_0.51.10 evaluate_1.0.5
#> [37] abind_1.4-8 stats4_4.6.1
#> [39] rmarkdown_2.32 matrixStats_1.5.0
#> [41] tools_4.6.1 htmltools_0.5.9