Predict OCR-driven Boolean rules from multi-omics datasets
ocrRBBR is an R package to infer Boolean rules linking chromatin accessibility (ATAC-seq peaks) to gene expression (RNA-seq) in both bulk and single-cell multiomic datasets. The package identifies combinations of OCRs (Open Chromatin Regions) that best predict the expression state of a target gene.
Features
- Supports bulk and single-cell multiome data with paired RNA-seq and ATAC-seq per cell or sample.
- Outputs interpretable Boolean rules with associated performance metrics.
- Enables parallel computing for faster processing of large datasets.
- RNA-seq (bulk): Data are quantile-normalized or TPM-normalized to adjust for sequencing depth and gene length.
- ATAC-seq (bulk): Signal intensities are quantile-normalized to adjust for sequencing depth differences across samples.
- Single-cell RNA-seq: LogNormalized with a scale factor of 10,000 (using Seurat).
- Single-cell ATAC-seq: Normalized using the ReadsInTSS method to adjust for cell-specific variations in sequencing depth.
- ocrRBBR is data-efficient and works well even with small sample sizes. Unlike neural networks, which require many parameters, ocrRBBR uses ridge regression with fewer parameters, making it suitable for datasets with limited samples. It has been tested on both bulk (85 cell types) and single-cell (9,834 cells) datasets and performs similarly on smaller single-cell datasets.
- In ocrRBBR, samples do not need to originate from the same tissue or cell type. When samples are from the same cell type or tissue, ocrRBBR partitions them into more homogeneous groups, each associated with a distinct Boolean rule within the inferred rule set.
install.packages("ocrRBBR")A toy dataset example is provided in Data/. Please see following examples for instructions.
ocrRBBR_single_cell() predicts Boolean rules for a given gene based on bulk-level multi-omics datasets (RNA-seq and ATAC-seq).
ocrRBBR_bulk(
rnaseq_data=rnaseq_data, # Matrix of RNA-seq gene expression (genes × samples)atacseq_data=atacseq_data, # Matrix of ATAC-seq peak accessibility (peaks × samples)gene_name="Rag2", # Gene for which Boolean rules will be inferredpeak_ids=peak_ids, # Candidate regulatory peaks for this genemax_feature=3, # Maximum number of OCRs allowed in a Boolean ruleslope=10, # Slope for sigmoid activation functionnum_cores=8, # Number of parallel workers for computationverbose=FALSE# Logical. If TRUE, progress messages and a progress bar are shown. Default is FALSE.
)Parameter Descriptions
# Required arguments# rnaseq_data A numeric matrix of RNA-seq expression values.# Rows correspond to genes, columns correspond to cell types or samples.# **Note:** ocrRBBR was tested using **quantile-normalized RNA-seq data**, but it should also work equally well on **TPM-normalized RNA-seq datasets**, provided the data is appropriately scaled across samples.## atacseq_data A numeric matrix of ATAC-seq signal intensities.# Rows correspond to peaks, columns correspond to cell types or samples.# Column names must match those of rnaseq_data.# **Note:** Similar to RNA-seq data, ocrRBBR is tested using **quantile-normalized ATAC-seq data** but is expected to work with other normalization methods, as long as the data distributions are comparable across samples.## gene_name A character string specifying the gene for which to infer Boolean rules.# peak_ids A vector of peak identifiers corresponding to rows in atacseq_data to be used as candidate regulatory regions for gene_name.## Optional arguments# max_feature An integer specifying the maximum number of input features allowed in a Boolean rule. Default is 3.# slope The slope parameter for the sigmoid activation function. Default is 10.# num_cores The number of parallel workers to use for computation. Adjust according to your system. Default is NA (automatic selection).# verbose Logical. If TRUE, progress messages and a progress bar are shown. Default is FALSE.link_peaks_to_tss() links ATAC-seq peaks to genes based on a user-defined window (±100kb by default) around the TSS (Transcription Start Site).
linked_peaks <- link_peaks_to_tss(
gtf_file = gtf_file, # Path to the GTF file with gene annotations
peaks_gr = peaks_gr, # GRanges object containing ATAC-seq peaks
gene_list = NA, # Optional: A list of specific genes to link peaks (default is NA, considering all genes)
tss_window = NA # Optional: A custom window size around the TSS (default is ±100kb)
)# Load the ATAC-seq data, RNA-seq data, and peak locations
data(multiome_human_mouse)
# List all objects in the current R environment
ls()
[1] "human_atacseq_data""human_cell_type""human_meta_data""human_peaks_gr""human_rnaseq_data""mouse_atacseq_data""mouse_peaks_gr""mouse_rnaseq_data"# Inspect the GRanges object peaks_gr, which contains peak consensus scores across mammalian genomes and associated peak p-values.
head(mouse_peaks_gr)
GRangesobjectwith6rangesand3metadatacolumns:seqnamesrangesstrand|peakIDphastCons_scoresmlog10_bestPvalue<Rle><IRanges><Rle>|<integer><numeric><numeric>
[1] chr13020761-3020811*|10.000.56
[2] chr13087201-3087251*|20.000.50
[3] chr13120084-3120134*|30.0710.80
[4] chr13121460-3121510*|40.153.02
[5] chr13372762-3372812*|50.031.31
[6] chr13399192-3399242*|60.062.39-------seqinfo:44sequencesfromanunspecifiedgenome; noseqlengths# Inspect the atacseq_data matrix, where rows correspond to peaks and columns correspond to cell types.# The values represent quantile-normalized ATAC-seq signal intensities.mouse_atacseq_data[1:5, 1:5]
LTHSC.34-.BMLTHSC.34+.BMSTHSC.150-.BMMPP4.135+.BMproB.CLP.BM2783454.372.570.904.843.592783460.410.100.900.113.212783520.410.101.850.831.152783532.360.711.850.112.522783540.411.640.900.830.90# Inspect the rnaseq_data matrix, where rows correspond to genes and columns correspond to cell types.# The values represent quantile-normalized RNA-seq signal intensities.mouse_rnaseq_data[ , 1:5]
LTHSC.34-.BMLTHSC.34+.BMSTHSC.150-.BMMPP4.135+.BMproB.CLP.BMRag21.0207951.021262.9853656.39795483.3949Spi1117.386316199.55519362.58834458.98175327.2625gene_name<-"Rag2"linked_peaks<- link_peaks_to_tss(
gtf_file= system.file("extdata", "gencode.vM25.annotation.sample.gtf", package="ocrRBBR"),
peaks_gr=mouse_peaks_gr,
gene_list=gene_name,
tss_window=100000# ±100kb
)
linked_peaks# A tibble: 83 × 7peakgene_idgene_namegene_typetranscript_idpeak_iddistance<chr><chr><chr><chr><chr><int><dbl>1chr2:101525902-101525952ENSMUSG00000032864Rag2protein_codingENSMUST0000004403127834598790.2chr2:101533903-101533953ENSMUSG00000032864Rag2protein_codingENSMUST0000004403127834690790.3chr2:101537021-101537071ENSMUSG00000032864Rag2protein_codingENSMUST0000004403127834787672.4chr2:101543845-101543895ENSMUSG00000032864Rag2protein_codingENSMUST0000004403127834880848.5chr2:101545763-101545813ENSMUSG00000032864Rag2protein_codingENSMUST0000004403127834978930.6chr2:101549272-101549322ENSMUSG00000032864Rag2protein_codingENSMUST0000004403127835075420.7chr2:101550344-101550394ENSMUSG00000032864Rag2protein_codingENSMUST0000004403127835174348.8chr2:101551169-101551219ENSMUSG00000032864Rag2protein_codingENSMUST0000004403127835273524.9chr2:101551591-101551641ENSMUSG00000032864Rag2protein_codingENSMUST0000004403127835373102.10chr2:101552534-101552584ENSMUSG00000032864Rag2protein_codingENSMUST0000004403127835472158.# ℹ 73 more rows# ℹ Use `print(n = ...)` to see more rows## -------------------------------# Get linked peaks for gene# -------------------------------linked_peaks_gene<-linked_peaks[linked_peaks$gene_name==gene_name, ]
if(nrow(linked_peaks_gene) ==0){
stop("No linked peaks found for gene: ", gene_name)
}
peak_ids<-linked_peaks_gene$peak_id# -------------------------------# This step removes ATAC-seq peaks with low signal intensities (based on p-values) or peaks not conserved across the mammalian genome.# This step helps reduce potential false positive predictions by ocrRBBR and can be omitted if desired.# -------------------------------peaks_gr_tmp<-mouse_peaks_gr[mouse_peaks_gr$peakID%in%peak_ids]
if(length(peaks_gr_tmp) ==0){
stop("No peak GRanges found for gene: ", gene_name)
}
phast_median<- median(peaks_gr_tmp$phastCons_scores, na.rm=TRUE)
pval_median<- median(peaks_gr_tmp$mlog10_bestPvalue, na.rm=TRUE)
peaks_gr_tmp<-peaks_gr_tmp[
(peaks_gr_tmp$phastCons_scores>phast_median) &
(peaks_gr_tmp$mlog10_bestPvalue>pval_median)
]
peak_ids<-peaks_gr_tmp$peakID
head(peak_ids)
[1] 278345278346278352278353278354278355boolean_rules<- ocrRBBR_bulk(mouse_rnaseq_data, mouse_atacseq_data, gene_name, peak_ids, max_feature=3, slope=10, num_cores=NA, verbose=TRUE)
Startingprocessingforgene:Rag2...Allinputcheckspassed.trainingprocessstartedwith8computingcores|====================|100%
head(boolean_rules)
Boolean_RuleR2BICRule_Coefficients1 [AND(278352,278381,278384),AND(¬278352,278381,278384),AND(278352,¬278381,278384)] 0.79-306.880.46:0.72:0.34:-2.45:-2.17:-0.4:-0.64:-0.052 [AND(278362,278381,278398),AND(¬278362,278381,278398),AND(278362,¬278381,278398),AND(¬278362,¬278381,¬278398)] 0.79-306.710.49:0.75:0.55:-1.6:-2.02:-0.48:-0.67:0.053 [AND(278381,278390,278398),AND(278381,¬278390,278398),AND(¬278381,¬278390,278398),AND(¬278381,278390,¬278398)] 0.79-306.690.18:-1.59:0.93:-0.35:0.58:0.13:-1.33:-0.854 [AND(278355,278381,278384),AND(¬278355,278381,278384),AND(278355,¬278381,278384)] 0.79-305.850.54:0.75:0.34:-1.99:-1.88:-0.49:-0.6:-0.085 [AND(278386,278398)] 0.75-304.720.71:-0.8:-1.04:-0.286 [278384] 0.73-304.670.54:-0.54# Top-ranked rule sets (high R², low BIC) highlight minimal yet robust combinations of OCRs that are most predictive of gene regulation.# Each row represents a Boolean rule set linking chromatin accessibility (OCRs; indexed by peak IDs) to gene expression.# Boolean_Rule shows the logical combinations of peaks (using AND and NOT, ¬) that best explain the gene’s expression pattern. Multiple AND clauses indicate alternative regulatory configurations that lead to similar transcriptional outcomes.# R² quantifies how much variance in gene expression is explained by the rule set; higher values indicate stronger explanatory power.# BIC (Bayesian Information Criterion) balances model fit and complexity; lower values indicate a better trade-off between accuracy and simplicity.# Rule_Coefficients provide the signed contribution of each Boolean rule in the fitted model. Positive coefficients indicate activating regulatory effects, while negative coefficients suggest repressive or inhibitory influences.A toy dataset example is provided in Data/. Please see following examples for instructions.
ocrRBBR_single_cell() infers OCR-driven Boolean regulatory rules for a given gene using single-cell paired RNA-seq and ATAC-seq multiome data, integrating chromatin accessibility states in gene-flanking regions with gene expression across cells.
ocrRBBR_single_cell(
rnaseq_data=rnaseq_data, # RNA-seq expression matrix (genes × cells)atacseq_data=atacseq_data, # ATAC-seq accessibility matrix (peaks × cells)gene_name="Rag2", # Gene for which Boolean rules are inferredpeak_ids=peak_ids, # Candidate regulatory peaks for the genemax_feature=3, # Maximum number of OCRs in a Boolean ruleslope=10, # Slope of the sigmoid activation functionnum_cores=8, # Number of parallel workersESS=500, # Effective sample sizemeta.data=meta.data, # Per-cell metadata (QC metrics)verbose=FALSE# Logical. If TRUE, progress messages and a progress bar are shown. Default is FALSE.
)Parameter Descriptions
# rnaseq_data A numeric matrix of RNA-seq expression values (genes × cells).# RNA-seq values are assumed to be normalized using Seurat’s LogNormalize method with a scale factor of 10,000:# NormalizeData(seurat_obj, normalization.method = "LogNormalize", scale.factor = 1e4)## atacseq_data A numeric matrix of ATAC-seq signal intensities (peaks × cells).# ATAC-seq counts are assumed to be normalized per cell using the ReadsInTSS method, where raw peak counts are divided by the number of Tn5 insertions within transcription start site (TSS) regions for each cell.# ReadsInTSS values are typically obtained from ArchR and applied as column-wise scaling factors.## gene_name A character string specifying the target gene for Boolean rule inference.## peak_ids A vector of peak identifiers corresponding to rows in atacseq_data, defining candidate regulatory regions for the gene.## meta_data A numeric matrix or data.frame containing additional per-cell covariates (rows = cells, columns = covariates), such as:# nCount_RNA# nFeature_RNA# Percentage of mitochondrial reads## Optional arguments# max_feature Maximum number of OCRs allowed in a Boolean rule. Default is 3.# slope Slope parameter of the sigmoid activation function used in the model. Default is 10.# num_cores Number of parallel workers for computation. Default is NA (automatic selection).# ESS Effective sample size of the single-cell data after accounting for noise and cell-to-cell correlation. Default is NA.# verbose Logical. If TRUE, progress messages and a progress bar are shown. Default is FALSE.link_peaks_to_tss() links ATAC-seq peaks to genes based on a user-defined window (±100kb by default) around the TSS (Transcription Start Site).
linked_peaks <- link_peaks_to_tss(
gtf_file = gtf_file, # Path to the GTF file with gene annotations
peaks_gr = peaks_gr, # GRanges object containing ATAC-seq peaks
gene_list = NA, # Optional: A list of specific genes to link peaks (default is NA, considering all genes)
tss_window = NA # Optional: A custom window size around the TSS (default is ±100kb)
)ESS() estimates the effective sample size (ESS) of single-cell RNA-seq data by accounting for correlation among cells within the same cell type.
ess <- ESS(
rnaseq_data = rnaseq_data, # Normalized RNA-seq expression matrix (genes × cells)
cell_type = cell_type # Cell-level metadata with cell type annotations
)Parameter Descriptions
# rnaseq_data must be a numeric matrix with genes as rows and cells as columns. RNA-seq data should be normalized (e.g., Seurat LogNormalize, scale factor = 10,000).# cell_type is a data frame with row names matching the columns of rnaseq_data and a required column named cell_type indicating cell identities.# Load the RData file containing the ATAC-seq data, RNA-seq data, and peak locations
data(multiome_human_mouse)
# Inspect the GRanges object peaks_gr, which contains peak consensus scores across mammalian genomes and associated peak p-values.
head(human_peaks_gr)
GRangesobjectwith6rangesand1metadatacolumn:seqnamesrangesstrand|peakID<Rle><IRanges><Rle>|<integer>
[1] chr110109-10357*|1
[2] chr1180730-181630*|2
[3] chr1191491-191736*|3
[4] chr1267816-268196*|4
[5] chr1586028-586373*|5
[6] chr1629721-630172*|6-------seqinfo:33sequencesfromanunspecifiedgenome; noseqlengths# Inspect the atacseq_data matrix, where rows correspond to peaks and columns correspond to cells.# The values represent normalized ATAC-seq counts per cell, using the ReadsInTSS method.human_atacseq_data[1:5, 1:5]
5x5sparseMatrixofclass"dgCMatrix"AAACAGCCAAGGAATC-1AAACAGCCAATCCCTT-1AAACAGCCAATGCGCT-1AAACAGCCACACTAAT-1AAACAGCCACCAACCG-183441.....83442.....83443.....83444.....83445.....# Inspect the rnaseq_data matrix, where rows correspond to genes and columns correspond to cells.# The values represent normalized RNA-seq signal intensities, calculated using Seurat's LogNormalize method with a scale factor of 10,000..human_rnaseq_data[1:5, 1:5]
5x5sparseMatrixofclass"dgCMatrix"AAACAGCCAAGGAATC-1AAACAGCCAATCCCTT-1AAACAGCCAATGCGCT-1AAACAGCCACACTAAT-1AAACAGCCACCAACCG-1MIR1302-2HG.....FAM138A.....OR4F5.....AL627309.1.....AL627309.3.....# Inspect meta data
head(human_meta_data)
nCount_RNAnFeature_RNApercent.mtAAACAGCCAAGGAATC-1838033087.470167AAACAGCCAATCCCTT-13771189610.527711AAACAGCCAATGCGCT-1687629046.457243AAACAGCCACACTAAT-1173384618.003462AAACAGCCACCAACCG-1541522826.500462AAACAGCCAGGATAAC-1275913536.922798# Inspect cell_type
head(human_cell_type)
cell_idcell_typeAAACAGCCAAGGAATC-1AAACAGCCAAGGAATC-1CD4NaiveAAACAGCCAATCCCTT-1AAACAGCCAATCCCTT-1CD4TCMAAACAGCCAATGCGCT-1AAACAGCCAATGCGCT-1CD4NaiveAAACAGCCACACTAAT-1AAACAGCCACACTAAT-1CD8NaiveAAACAGCCACCAACCG-1AAACAGCCACCAACCG-1CD8NaiveAAACAGCCAGGATAAC-1AAACAGCCAGGATAAC-1CD4Naive# Generate a list of peaks located within a ±100kb window around each gene.gene_name<-"CD74"linked_peaks<- link_peaks_to_tss(
gtf_file= system.file("extdata", "gencode.v48.basic.annotation.sample.gtf", package="ocrRBBR"),
peaks_gr=human_peaks_gr,
gene_list=gene_name, tss_window=250000# ±250kb
)
linked_peaks# A tibble: 43 × 7peakgene_idgene_namegene_typetranscript_idpeak_iddistance<chr><chr><chr><chr><chr><int><dbl>1chr5:150175660-150175697ENSG00000019582CD74protein_codingENST00000353334834412370712chr5:150180579-150181363ENSG00000019582CD74protein_codingENST0000035333483442231778.3chr5:150185269-150186122ENSG00000019582CD74protein_codingENST00000353334834432270544chr5:150208861-150209820ENSG00000019582CD74protein_codingENST00000353334834442034095chr5:150217881-150218594ENSG00000019582CD74protein_codingENST00000353334834451945126chr5:150234576-150235054ENSG00000019582CD74protein_codingENST0000035333483446177934.7chr5:150264137-150264569ENSG00000019582CD74protein_codingENST0000035333483447148396.8chr5:150350443-150351368ENSG00000019582CD74protein_codingENST0000035333483448618449chr5:150352283-150352763ENSG00000019582CD74protein_codingENST000003533348344960226.10chr5:150356700-150359453ENSG00000019582CD74protein_codingENST000003533348345054673# ℹ 33 more rows# ℹ Use `print(n = ...)` to see more rowslinked_peaks_gene<-linked_peaks[linked_peaks$gene_name==gene_name, ]
if(nrow(linked_peaks_gene) ==0){
stop("No linked peaks found for gene: ", gene_name)
}
peak_ids<-linked_peaks_gene$peak_id
head(peak_ids)
[1] 834418344283443834448344583446ess_value<- ESS(rnaseq_data=human_rnaseq_data, cell_type=human_cell_type)
Processing:BintermediateProcessing:BmemoryProcessing:BnaiveProcessing:CD14MonoProcessing:CD16MonoProcessing:CD4NaiveProcessing:CD4ProliferatingProcessing:CD4TCMProcessing:CD4TEMProcessing:CD8NaiveProcessing:CD8TCMProcessing:CD8TEMProcessing:cDC2Processing:dnTProcessing:gdTProcessing:HSPCProcessing:MAITProcessing:NKProcessing:NKProliferatingProcessing:NK_CD56brightProcessing:pDCProcessing:PlasmablastProcessing:Tregboolean_rules<- ocrRBBR_single_cell(human_rnaseq_data, human_atacseq_data, gene_name, peak_ids, slope=6, num_cores=NA, ESS=261, meta_data=human_meta_data, verbose=TRUE)
Startingprocessingforgene:CD74...Allinputcheckspassed.trainingprocessstartedwith8computingcores|====================|100%
head(boolean_rules)
Boolean_RuleR2BICRule_Coefficients1 [83456] 0.40220.130.39:-0.382 [AND(83456,83458),AND(83456,¬83458)] 0.42222.680.55:-0.09:0.17:-0.683 [AND(83456,83460),AND(83456,¬83460)] 0.42224.480.4:-0.08:0.3:-0.684 [AND(83456,83475),AND(83456,¬83475)] 0.41226.150.3:-0.08:0.49:-0.645 [AND(83456,83482),AND(83456,¬83482)] 0.41226.650.62:-0.42:0.17:-0.36# Top-ranked rule sets (high R², low BIC) highlight minimal yet robust combinations of OCRs that are most predictive of gene regulation.# Each row represents a Boolean rule set linking chromatin accessibility (OCRs; indexed by peak IDs) to gene expression.# Boolean_Rule shows the logical combinations of peaks (using AND and NOT, ¬) that best explain the gene’s expression pattern. Multiple AND clauses indicate alternative regulatory configurations that lead to similar transcriptional outcomes.# R² quantifies how much variance in gene expression is explained by the rule set; higher values indicate stronger explanatory power.# BIC (Bayesian Information Criterion) balances model fit and complexity; lower values indicate a better trade-off between accuracy and simplicity.# Rule_Coefficients provide the signed contribution of each Boolean rule in the fitted model. Positive coefficients indicate activating regulatory effects, while negative coefficients suggest repressive or inhibitory influences.