scMOSAIC is an R package for processing targeted PacBio long-read sequencing generated from 10x Genomics Multiome ATAC libraries and linking locus-specific DNA variation to matched single-nucleus transcriptomes.
The package was developed for single-cell measurement of HTT CAG repeat length in Huntington's disease, but the general strategy can be adapted to other repeat expansions or genomic loci that can be recovered from ATAC-derived DNA.
Important distinction: this repository contains both (1) the reusable scMOSAIC R package and (2) an
examples/directory containing analysis scripts used for specific studies. The package provides general-purpose functions; the study-specific scripts document how those functions were applied to reproduce a particular analysis.
scMOSAIC combines three pieces of information from the same nucleus:
- Gene expression from the 10x Genomics Multiome RNA library, used for cell-type identification.
- ATAC-derived genomic DNA, which retains the 10x cell barcode.
- Targeted PacBio long-read sequencing, used to resolve the genomic locus of interest across repetitive or otherwise difficult-to-sequence sequence.
For the HTT application, the computational workflow extracts the 10x barcode and DNA insert from each PacBio read, identifies reads spanning HTT exon 1, measures the uninterrupted CAG repeat tract, assigns reads to samples/cells, and summarizes repeat length at the nucleus, donor, cell-type, or anatomical-region level.
The experimental scMOSAIC workflow branches from 10x Multiome pre-amplified material into:
- a gene-expression track for cell identity;
- an optional standard ATAC track for chromatin accessibility; and
- a long-read ATAC track for targeted repeat/locus measurement.
The long-read track uses targeted hybridization capture followed by PacBio sequencing.
scMOSAIC/
├── R/ # Reusable R-package functions
│ ├── demultiplex_samples.R
│ ├── parse_reads.R
│ ├── count_repeats.R
│ ├── demultiplex_reads.R
│ └── calculate_alleles.R
├── data/ # Small package datasets / lookup tables
├── data-raw/ # Scripts or source material used to construct package data
├── man/ # R documentation
├── tests/ # Basic package tests
├── visuals/ # Read-structure diagrams and other documentation figures
scMOSAIC is currently distributed through GitHub.
# install.packages("remotes")
remotes::install_github("KellisLab/scMOSAIC")
library(scMOSAIC)You can also install with devtools:
# install.packages("devtools")
devtools::install_github("KellisLab/scMOSAIC")
library(scMOSAIC)First confirm that the repository is reachable from the machine running R:
url("https://api.github.com/repos/KellisLab/scMOSAIC")Then update the installation helper and retry:
install.packages("remotes")
remotes::install_github(
"KellisLab/scMOSAIC",
upgrade = "never"
)On institutional clusters, GitHub API access, HTTPS proxy settings, firewalls, or outdated SSL/certificate libraries can also prevent install_github() from downloading the repository.
A small example dataset is included with the package so that users can test the core parsing workflow without downloading the full study dataset.
library(scMOSAIC)
# Parse long reads into DNA inserts, 10x barcodes,
# sample-index information, and read orientation.
read_info <- parse_reads(
input = scMOSAIC::Fitzwalter2025,
out.file = "none"
)
head(read_info)
# Count CAG repeats, collapse reads to nuclei,
# and map ATAC barcodes onto the corresponding RNA barcode.
cell_info <- demultiplex_reads(
reads = read_info,
out.file = "None"
)
head(cell_info)demultiplex_reads() returns cell-level repeat summaries including fields such as:
atac
cag_repeats_mean
cag_repeats_median
cag_repeats_mode
cag_repeats_var
cag_repeats_sd
number_of_reads
sample
barcode_sample
The exact fields may evolve as the package is developed.
PacBio processing can be thought of as seven major stages:
- Orient reads from P5 to P7.
- Demultiplex sequencing runs by sample index.
- Identify the expected 10x library structural elements.
- Retain reads spanning the target genomic locus.
- Measure the repeat or other locus-specific feature.
- Assign reads to individual nuclei using 10x barcodes.
- Assign alleles and calculate repeat-expansion metrics.
For the published HTT application, PacBio reads were first oriented using the PacBio Lima workflow before downstream processing in scMOSAIC.
For multiplexed PacBio runs, demultiplex_samples() can split an oriented FASTQ file into sample-specific files using the sample-index sequence at the end of the read.
library(scMOSAIC)
dir.create("results/sample_fastqs", recursive = TRUE, showWarnings = FALSE)
demultiplex_samples(
input = "pacbio_oriented.fastq",
out.folder = "results/sample_fastqs/",
run.name = "run_01"
)Important arguments include:
| Argument | Description |
|---|---|
input |
Input FASTQ/FASTA file containing oriented PacBio reads |
out.folder |
Existing directory for sample-level output files |
run.name |
Name used to identify the sequencing run |
SampleIndexes |
Optional subset of expected sample indexes |
Read2N.max.mm |
Maximum mismatch allowed for the Read 2N sequence |
sampleid.length |
Expected length of the sample-index barcode |
sampleid.max.mm |
Maximum mismatch allowed during sample-index matching |
min.reads |
Minimum reads required for a sample to be retained/written |
file.format |
Sequence-file format, e.g. "fastq" |
n.cores |
Number of CPU cores used for parallel processing |
---
## 2. Parse the 10x/PacBio read structure
`parse_reads()` is the core read-parsing function. It filters reads for expected structural components, identifies reads spanning the target locus, and extracts the insert, cell barcode, sample barcode, and strand.
A minimal example is:
```r
read_info <- parse_reads(
input = "sample.fastq",
out.file = "results/sample_parsed.csv",
plot = FALSE
)
The function can accept:
- a path to one sequence file;
- a character vector containing multiple file paths; or
- a
DNAStringSet.
For the HTT implementation, default sequences are provided for:
- P5;
- spacer;
- Read 1N;
- Read 2N;
- P7;
- the sequence upstream of the HTT CAG tract; and
- the sequence downstream of the HTT CAG tract.
The defaults reproduce the HTT-oriented implementation. Users adapting scMOSAIC to another locus should replace the target-flanking sequences and empirically optimize mismatch tolerances.
read_info <- parse_reads(
input = "sample.fastq",
out.file = "results/custom_locus_parsed.csv",
HTTexon1_p1 = "YOUR_UPSTREAM_FLANK_SEQUENCE",
HTTexon1_p2 = "YOUR_DOWNSTREAM_FLANK_SEQUENCE",
HTTexon1_p1.max.mm = 1,
HTTexon1_p2.max.mm = 1,
plot = FALSE
)Although the current argument names contain HTTexon1, these arguments simply define the target-flanking sequences used to identify locus-spanning reads.
Do not assume the HTT mismatch thresholds are optimal for another target. Sequence length, local composition, repeat structure, sequencing chemistry, and capture design can all affect the recovery/fidelity tradeoff.
count_cag_repeats() operates on parsed reads and expects at least two columns:
insertstrand
read_info$cag_repeats <- count_cag_repeats(
read_info,
min.repeats = 3
)
summary(read_info$cag_repeats)For forward-strand reads the function searches for contiguous CAG runs; for reverse-strand reads it searches for the reverse-complement CTG sequence.
The current implementation sums repeat runs exceeding the min.repeats threshold.
Example visualization:
hist(
read_info$cag_repeats,
breaks = 100,
xlab = "CAG repeat count",
main = "Read-level HTT CAG repeat lengths"
)Multiple PacBio reads may map to the same 10x nucleus. demultiplex_reads() combines read-level calls into cell-level summaries.
cell_info <- demultiplex_reads(
reads = read_info,
out.file = "None"
)
head(cell_info)To save the table directly:
demultiplex_reads(
reads = read_info,
out.file = "results/cell_level_repeat_calls.csv"
)The function calculates per-cell summary statistics including mean, median, mode, variance, standard deviation, and number of long reads.
The output also uses package barcode lookup information to connect the ATAC barcode to the corresponding RNA barcode, allowing repeat measurements to be joined to the matched 10x gene-expression profile.
For heterozygous repeat-expansion disorders such as Huntington's disease, a donor commonly contributes reads from both a normal inherited allele and an expanded inherited allele.
The scMOSAIC allele-calling strategy is:
- aggregate repeat calls at the donor level;
- remove extreme outlier repeat lengths;
- estimate a two-mode repeat-length distribution;
- identify the antimode separating the two distributions;
- assign values below the antimode to the normal allele and values at/above the antimode to the expanded allele; and
- summarize repeat lengths and expansion above the inherited expanded allele.
The package contains calculate_alleles() for this purpose.
library(scMOSAIC)
library(ggplot2)
library(locmode)
#for example lets make some toy wt and hd distributions
#donor 1: 10,000 reads collected, and has a inherited normal allele of 15 and a inheritied expanded allele of 40
d1 <- data.frame(repeats = c(rnorm(10000, mean = 15, sd = 2),
rnbinom(10000, 40, .5)),
donor = "donor_1")
#donor 2: 10,000 reads collected, and has a inherited normal allele of 16 and a inheritied expanded allele of 45
d2 <- data.frame(repeats = c(rnorm(10000, mean = 16, sd = 2),
rnbinom(10000, 45, .5)),
donor = "donor_2")
#donor 2: 10,000 reads collected, and has a inherited normal allele of 15 and a inheritied normal allele of 18
d3 <- data.frame(repeats = c(rnorm(10000, mean = 15, sd = 2),
rnorm(10000, mean = 20, sd = 2)),
donor = "donor_3")
#combine all three donors together
d <- do.call(rbind, list(d1, d2, d3))
#now run calculate allele on all donors to approximate their hd and wt alleles
allele_info <- lapply(levels(factor(d$donor)), function(donor){
res <- calculate_alleles(d[d$donor == donor,]$repeats, plot = TRUE)
})
names(allele_info) <- levels(factor(d$donor))
#view results
allele_info
Expected output is a small table containing the estimated normal-allele mode, antimode, and expanded-allele mode.
In the scMOSAIC analysis framework, expansion can be summarized relative to the major inherited expanded allele. Conceptually:
expanded_reads <- donor_reads$cag_repeats[
donor_reads$cag_repeats >= antimode
]
expansion_index <- mean(expanded_reads - expanded_allele_mode)For published analyses, calculate allele boundaries at the donor level before comparing expansion metrics across cell types or anatomical regions.
Development note:
calculate_alleles()is currently under active development. If you are using the current development branch, inspect/test donor-level allele calls before applying them to a large dataset.
After cell-level repeat calls have been produced, they can be joined to RNA-derived metadata using the shared 10x barcode/sample identifier.
For example, if rna_metadata contains one row per RNA nucleus:
# Example only; adjust column names to match your Cell Ranger / Seurat / Scanpy object.
merged <- merge(
rna_metadata,
cell_info,
by = "barcode_sample",
all.x = TRUE
)
table(is.na(merged$cag_repeats_mode))The resulting table can be used to compare somatic repeat expansion across:
- cell classes;
- transcriptional subtypes;
- anatomical regions;
- donors;
- disease states; or
- continuous transcriptional phenotypes.
scMOSAIC itself focuses on processing the long-read measurements and barcode linkage; downstream single-cell normalization, clustering, cell-type annotation, differential expression, and visualization can be performed using standard tools such as Seurat or Scanpy.
| Function | Purpose |
|---|---|
demultiplex_samples() |
Split multiplexed PacBio reads into samples using sample-index sequences |
parse_reads() |
Filter reads on expected structural elements and extract locus-spanning inserts, barcodes, sample IDs, and strand |
count_cag_repeats() |
Count CAG/CTG repeat tracts in parsed DNA inserts |
demultiplex_reads() |
Collapse read-level repeat calls into nucleus-level summaries and map ATAC to RNA barcodes |
calculate_alleles() |
Estimate normal/expanded allele modes and the antimode separating them |
R help pages can be opened with:
?demultiplex_samples
?parse_reads
?count_cag_repeats
?demultiplex_reads
?calculate_allelesA schematic of the expected 10x-derived long-read structure:
The parser uses these structural elements to determine whether a read is suitable for downstream locus and barcode extraction.
The overall framework is not intrinsically specific to HTT. A new target should satisfy two broad requirements:
- the locus must be recoverable in the ATAC-derived DNA library; and
- the targeted long-read library must contain enough sequence on either side of the feature of interest to identify locus-spanning reads reliably.
For a new locus, users will generally need to modify:
- the hybridization-capture probe design;
- upstream/downstream target-flanking sequences;
- mismatch tolerances;
- repeat-counting or variant-calling logic; and
- filtering thresholds.
The HTT flanking-sequence lengths and mismatch thresholds were empirically optimized for that locus and should not be treated as universal defaults.
10x Multiome
|
+-- RNA library --------------------------> cell type / cell state
|
+-- ATAC-derived DNA
|
+-- standard ATAC (optional)
|
+-- targeted hybridization capture
|
+-- PacBio HiFi sequencing
|
+-- orient reads
+-- sample demultiplexing
+-- read-structure parsing
+-- target-locus filtering
+-- repeat / variant calling
+-- 10x barcode assignment
+-- cell-level summarization
+-- allele / expansion metrics
|
+-- join to RNA metadata
The R package covers the computational portion of scMOSAIC. Wet-lab generation of scMOSAIC libraries requires the accompanying experimental protocol, including:
- compatible nuclei isolation;
- 10x Genomics Multiome processing;
- modified long-read ATAC library preparation;
- targeted hybridization capture;
- capture QC;
- SMRTbell library preparation; and
- PacBio long-read sequencing.
Users attempting to generate new scMOSAIC data should follow the experimental protocol rather than treating this README as a wet-lab SOP.
A few practical points are especially important:
- scMOSAIC long-read calls are sparse relative to the total number of profiled nuclei;
- multiple reads from the same nucleus should be summarized rather than treated as independent cells;
- repeat distributions should be inspected separately for each donor before allele assignment;
- low-read cells can have unstable estimates of mean, median, or mode repeat length;
- sample-index and 10x-barcode whitelist matching should be checked carefully;
- target-flank filtering trades sensitivity against false-positive locus recovery; and
- sequencing/capture batches should be evaluated for consistent enrichment and read-length distributions.
For the HTT workflow, read-level QC, donor-level allele distributions, and cell-level read depth should all be inspected before biological comparisons are made.
scMOSAIC is research software under active development.
The package has been tested on the datasets used during method development and now includes basic automated tests, but users may encounter edge cases in new sequencing runs, library chemistries, targets, or compute environments.
If you encounter a reproducible problem, please open a GitHub issue and include:
- the function call;
- the error message;
sessionInfo();- a minimal example if possible; and
- enough information about the input format to reproduce the issue without sharing protected human data.
If you use scMOSAIC, please cite the experimental protocol and the corresponding application paper.
Fitzwalter BE, Fass SB, Cameron J, Linville RM, Pineda SS, Kellis M, Heiman M. (2026).
Single-Cell Multi-Omic Somatic Analysis with Integrated Capture (scMOSAIC), Version 1.
Zenodo. https://doi.org/10.5281/zenodo.22083276
Linville RM, James BT, Galani K, Ho L-L, Shin JH, Oliver E, Fass SB, Cameron JC, Fitzwalter BE, Bock R, et al. (2026).
Cross-species single-cell atlas of the striatum defines cell type and subregion disease vulnerabilities.
Cell 189, 1–32. https://doi.org/10.1016/j.cell.2026.08.006
Please cite the protocol version corresponding to the version used in your experiment.
Bug reports, documentation improvements, and pull requests are welcome through the GitHub repository:
https://github.com/KellisLab/scMOSAIC
When proposing changes to the locus-specific parsing logic, please indicate whether the change is intended to:
- preserve compatibility with the published HTT workflow; or
- generalize scMOSAIC to another locus or sequencing design.
This distinction helps keep the validated HTT implementation reproducible while allowing the software to expand to new applications.