Chapter 10 Ambient contamination
10.1 Motivation
Ambient contamination is often present in scRNA-seq datasets generated by high-throughput protocols. During dissociation, cell lysis releases RNA molecules into the ambient solution. These molecules are captured and sequenced in each cell’s reaction chamber, contributing counts to genes that are not actually expressed in that cell. This is most pronounced when the ambient solution contains very high-abundance transcipts, e.g., hemoglobin from lysed erythrocytes.
In most situations, ambient contamination is just a mild nuisance that can be ignored. Most of our single-cell analysis focuses on the differences between cell types/states, so the addition of the same contaminating expression profile to each cell mostly cancels out in these comparisons. There is some potential for artificial differences when the same volume of contaminating solution is added to cells with different total RNA content, but these should be relatively minor compared to the log-fold changes from marker genes between cell types.
For DE analyses of multi-condition datasets (Section 9.2), ambient contamination becomes more problematic if the ambient expression profile differs across conditions. Such differences in the ambient profile are not uncommon in studies involving a strong experimental perturbation. If a gene is upregulated in one subpopulation in one condition, lysis of that subpopulation’s cells will increase the abundance of those transcripts in the ambient solution. Similarly, changes in the abundance of a subpopulation may introduce differences in the ambient expression of its marker genes. As a result, we might detect spurious DE genes that are driven by differences in the ambient profiles rather than any intrinsic change in gene regulation. This is obviously undesirable.
10.2 Spurious DE from contamination
To illustrate, let’s consider the Tal1-knockout (KO) mouse gastrulation data from Pijuan-Sala et al. (2019). Tal1 is a transcription factor that has known roles in erythroid differentiation. This study knocked out Tal1 in the early mouse embryo to determine whether blocking the erythroid lineage diverts cells to other developmental fates.
## class: SingleCellExperiment
## dim: 29453 56122
## metadata(0):
## assays(1): counts
## rownames(29453): ENSMUSG00000051951 ENSMUSG00000089699 ...
## ENSMUSG00000095742 tomato-td
## rowData names(2): ENSEMBL SYMBOL
## colnames(56122): cell_1 cell_2 ... cell_56121 cell_56122
## colData names(9): cell barcode ... pool sizeFactor
## reducedDimNames(1): pca.corrected
## mainExpName: NULL
## altExpNames(0):We compute pseudo-bulk profiles using the author-provided cell type labels and the sample identities. This experiment contains two technical replicates for each of its two conditions, i.e., wild-type (WT) and KO embryos. We’ll pretend that our replicates are biological so that we can proceed to the DE analysis, though readers should keep in mind that many of the \(p\)-values will be greatly understated in this section.
library(scrapper)
pseudo.tal1 <- aggregateAcrossCells.se(
sce.tal1,
factor=colData(sce.tal1)[,c("sample", "celltype.mapped")]
)We perform a voom-based DE analysis between the WT and KO pseudo-bulk samples for “spinal cord” cells.
We focus on the spinal cord as this should be mostly unrelated to the erythroid lineage.
Thus, we would not expect to see a strong change in expression upon loss of Tal1.
library(edgeR)
pseudo.spinal <- pseudo.tal1[,pseudo.tal1$celltype.mapped == "Spinal cord"]
y.spinal <- DGEList(
assay(pseudo.spinal, "sums"),
samples=as.data.frame(colData(pseudo.spinal)),
genes=as.data.frame(rowData(pseudo.spinal))
)
# KO samples are those labelled with tdTomato.
keep.spinal <- filterByExpr(y.spinal, group=y.spinal$samples$tomato)
y.spinal <- y.spinal[keep.spinal,]
y.spinal <- normLibSizes(y.spinal)
# Using a log-fold change threshold of 0.5 here, to avoid prioritizing genes
# with very small log-fold changes when the variance is also very low.
library(limma)
design.spinal <- model.matrix(~tomato, y.spinal$samples)
v.spinal <- voomLmFit(y.spinal, design.spinal, sample.weights=TRUE)
fit.spinal <- treat(v.spinal, robust=TRUE, lfc=0.5)
res.spinal <- topTreat(fit.spinal, sort.by="none", n=Inf, coef="tomatoTRUE")
head(res.spinal[order(res.spinal$P.Value),])## ENSEMBL SYMBOL logFC AveExpr t
## ENSMUSG00000086503 ENSMUSG00000086503 Xist -6.688946 5.8780449 -34.06156
## ENSMUSG00000052217 ENSMUSG00000052217 Hbb-bh1 -8.128804 5.9893873 -31.26849
## ENSMUSG00000055609 ENSMUSG00000055609 Hba-x -7.448039 5.6788928 -30.64385
## ENSMUSG00000052187 ENSMUSG00000052187 Hbb-y -9.115354 4.6302664 -19.51437
## ENSMUSG00000069919 ENSMUSG00000069919 Hba-a1 -8.255853 3.3076057 -17.47957
## tomato-td tomato-td tomato-td 7.263182 0.0972628 15.84265
## P.Value adj.P.Val
## ENSMUSG00000086503 1.651289e-28 1.736990e-24
## ENSMUSG00000052217 2.942424e-27 1.547568e-23
## ENSMUSG00000055609 5.749406e-27 2.015933e-23
## ENSMUSG00000052187 1.513873e-20 3.981107e-17
## ENSMUSG00000069919 4.826428e-19 1.015384e-15
## tomato-td 1.881580e-17 3.298723e-14Most of the top DE genes are the hemoglobins, which are downregulated in the KO cells. This is rather surprising as the spinal cord cells do not belong to the erythroid lineage and should not express hemoglobins at all. The most sober explanation is that the WT samples contain hemoglobin transcripts in the ambient solution due to leakage from erythrocytes or their precursors. In contrast, the erythroid lineage is blocked in the KO samples so no hemoglobin transcripts are present.
10.3 Estimating the ambient profile
To quantify the impact of ambient contamination on the interpretation of our DE results, we obtain an estimate of the ambient expression profile from the unfiltered count matrix for each sample. We consider all barcodes with total counts at or below 100 to represent empty droplets (A. Lun et al. 2019), and we sum the counts for each gene across these barcodes to obtain an expression vector representing the ambient profile for each sample.
# Extracting an ambient profile from the unfiltered count matrix,
# i.e., prior to any cell calling from Cellranger or emptyDrops.
ambient.tal1 <- list()
for (s in pseudo.spinal$sample) {
raw.tal1 <- Tal1ChimeraData(type="raw", samples=s)[[1]]
mat <- counts(raw.tal1)
ambient.tal1[[s]] <- SummarizedExperiment(
list(sums=cbind(rowSums(mat[,colSums(mat) < 100]))),
rowData=rowData(raw.tal1)
)
}
ambient.tal1 <- do.call(cbind, ambient.tal1)We repeat our DE analysis on the ambient count matrix, which yields a similar set of DE genes as in res.spinal.
This supports our hypothesis that ambient contamination drives the apparent expression of hemoglobins in the WT cells and their downregulation in KO cells.
y.ambient <- DGEList(
assay(ambient.tal1),
genes=rowData(ambient.tal1),
samples=colData(pseudo.spinal) # re-using the colData, as it's the same samples.
)
keep.ambient <- filterByExpr(y.ambient, group=y.ambient$samples$tomato)
y.ambient <- y.ambient[keep.ambient,]
y.ambient <- normLibSizes(y.ambient)
design.ambient <- model.matrix(~tomato, y.ambient$samples)
v.ambient <- voomLmFit(y.ambient, design.ambient)
fit.ambient <- treat(v.ambient, robust=TRUE, lfc=0.5)
res.ambient <- topTreat(fit.ambient, sort.by="none", n=Inf, coef="tomatoTRUE")
head(res.ambient[order(res.ambient$P.Value),])## ENSEMBL SYMBOL logFC AveExpr t
## ENSMUSG00000052187 ENSMUSG00000052187 Hbb-y -5.245640 11.142235 -142.20912
## ENSMUSG00000055609 ENSMUSG00000055609 Hba-x -4.805669 11.667507 -134.13325
## ENSMUSG00000069919 ENSMUSG00000069919 Hba-a1 -4.640627 9.369705 -100.34617
## ENSMUSG00000069917 ENSMUSG00000069917 Hba-a2 -4.498210 8.198450 -72.51491
## ENSMUSG00000027562 ENSMUSG00000027562 Car2 -3.495525 7.710538 -45.12594
## ENSMUSG00000086503 ENSMUSG00000086503 Xist -4.359582 6.320087 -40.10206
## P.Value adj.P.Val
## ENSMUSG00000052187 4.487863e-42 5.034933e-38
## ENSMUSG00000055609 2.318553e-41 1.300592e-37
## ENSMUSG00000069919 8.067189e-38 3.016860e-34
## ENSMUSG00000069917 7.310262e-34 2.050346e-30
## ENSMUSG00000027562 4.114727e-28 9.232625e-25
## ENSMUSG00000086503 1.088653e-26 2.035600e-23Alternatively, we could use other definitions for empty droplets to obtain slightly different estimates of the ambient profile.
For example, more recent versions of CellRanger use barcodes with the 45,000-th to 90,000-th largest total counts to construct their ambient profiles.
If an ambient profile is not available, we can also approximate it with the sum of counts across all cells in each sample.
This assumes that all subpopulations contribute equally to the ambient solution, e.g., all cell types have the same propensity for damage/lysis.
10.4 Dealing with ambient DE
10.4.1 Removing ambient DE genes
A simple fix is to just remove the ambient DE genes from the spinal cord results.
Any ambient DE is considered to compromise the interpretation of the affected gene’s results, so we just filter it out to be on the safe side.
One might think to filter out genes in res.spinal that are significant in res.ambient,
but a better approach is to combine the ambient and spinal cord counts into a single object and repeat the relevant comparisons.
This avoids differences in detection power between comparisons due to fluctuations in the variance estimates and/or differences in abundance-based filtering.
Otherwise, a gene could have the same log-fold change in both comparisons but only be detected as significant in one of the comparisons.
y.combined <- DGEList(
counts = cbind(assay(pseudo.spinal, "sums"), assay(ambient.tal1, "sums")),
genes = as.data.frame(rowData(pseudo.spinal))
)
y.combined$samples$group <- sprintf(
"%s.%s",
ifelse(rep(pseudo.spinal$tomato, 2), "KO", "WT"),
rep(c("spinal", "ambient"), each=ncol(pseudo.spinal))
)
keep.combined <- filterByExpr(y.combined, group=y.combined$samples$mode)
y.combined <- y.combined[keep.combined,]
y.combined <- normLibSizes(y.combined)As we’ve combined all of the counts into a single analysis, the design matrix is slightly more complicated.
We also need to set contrasts= in each call to voomLmFit() to specify which comparison we want to perform.
design.combined <- model.matrix(~0 + group, y.combined$samples)
con.combined.spinal <- makeContrasts(groupKO.spinal - groupWT.spinal, levels=design.combined)
v.combined.spinal <- voomLmFit(y.combined, design.combined, sample.weights=TRUE, contrasts=con.combined.spinal)
fit.combined.spinal <- treat(v.combined.spinal, robust=TRUE, lfc=0.5)
res.combined.spinal <- topTreat(fit.combined.spinal, sort.by="none", n=Inf)
con.combined.ambient <- makeContrasts(groupKO.ambient - groupWT.ambient, levels=design.combined)
v.combined.ambient <- voomLmFit(y.combined, design.combined, sample.weights=TRUE, contrasts=con.combined.ambient)
fit.combined.ambient <- treat(v.combined.ambient, robust=TRUE, lfc=0.5)
res.combined.ambient <- topTreat(fit.combined.ambient, sort.by="none", n=Inf)We filter out genes with any suggestion of DE in the ambient comparison using a \(p\)-value threshold of 0.05, i.e., the false negative rate of incorrectly removing non-DE genes is 5%. The exception is that of genes where the signs of the log-fold changes are not the same between comparisons, as a change in the spinal cord expression could not be caused by an ambient difference in the other direction.
# Putting relevant statistics in the same table for easier inspection.
res.combined <- res.combined.spinal
res.combined$ambient.P.Value <- res.combined.ambient$P.Value
res.combined$ambient.logFC <- res.combined.ambient$logFC
keep.combined <- res.combined$ambient.P.Value > 0.05 |
sign(res.combined$ambient.logFC) != sign(res.combined$logFC)
summary(keep.combined)## Mode FALSE TRUE
## logical 586 11049res.combined <- res.combined[keep.combined,]
res.combined$adj.P.Val <- p.adjust(res.combined$P.Value, method="BH")
head(res.combined[order(res.combined$P.Value),])## ENSEMBL SYMBOL logFC AveExpr
## ENSMUSG00000042814 ENSMUSG00000042814 Mcts2 1.0901073 6.255298
## ENSMUSG00000039936 ENSMUSG00000039936 Pik3cd -4.0077289 -0.411373
## ENSMUSG00000086429 ENSMUSG00000086429 Gt(ROSA)26Sor 0.8738376 5.859638
## ENSMUSG00000027406 ENSMUSG00000027406 Idh3b 0.7278997 6.859530
## ENSMUSG00000048327 ENSMUSG00000048327 Ckap2l 0.8750045 5.056830
## ENSMUSG00000029499 ENSMUSG00000029499 Pxmp2 -1.0996698 3.358019
## t P.Value adj.P.Val ambient.P.Value
## ENSMUSG00000042814 7.864491 3.697419e-10 4.085278e-06 0.3458162
## ENSMUSG00000039936 -4.332072 4.584555e-05 2.532738e-01 0.4353765
## ENSMUSG00000086429 4.118059 8.511292e-05 3.134709e-01 0.9995373
## ENSMUSG00000027406 3.877086 1.784586e-04 4.929472e-01 0.7599793
## ENSMUSG00000048327 3.755445 2.576580e-04 5.693727e-01 0.1722221
## ENSMUSG00000029499 -3.257494 1.099132e-03 1.000000e+00 0.9386596
## ambient.logFC
## ENSMUSG00000042814 0.5292301
## ENSMUSG00000039936 -0.5929610
## ENSMUSG00000086429 0.2312355
## ENSMUSG00000027406 0.4580299
## ENSMUSG00000048327 0.6036546
## ENSMUSG00000029499 -0.1981913Now, the hemoglobins and its friends are absent from the set of top DE genes. This removal approach is robust as it does not make quantitative assumptions about the contributions of ambient contamination to the spinal cord DE. However, it is also conservative and may remove genes that are clearly driven by genuine changes in expression within the cell population. For example, consider a DE gene that was detected in both the ambient and spinal cord comparisons but has a much larger log-fold change in the latter. Our approach would remove this gene despite the fact that the spinal cord DE would not be explained by the ambient DE.
10.4.2 More extreme log-fold changes
For the sake of completeness, we’ll also demonstrate a more sophisticated approach based on an interaction model. Consider a gene that is upregulated upon KO within the spinal cord subpopulation. We expect that the log-fold change for this gene in the spinal cord comparison should be more extreme than that in the ambient comparison. This is because any DE that leaks into the ambient solution from this subpopulation is effectively “diluted” by the non-DE contributions from all of the other subpopulations. Conversely, we ignore any DE gene where the spinal cord log-fold change is equal to or weaker than the ambient log-fold change as this could be driven by ambient contamination.
The interaction approach is more powerful than our previous removal approach as the former takes advantage of the magnitudes of the differences in each comparison. By comparing the log-fold changes, we can recover DE genes in the spinal cord comparison even if they are also significant in the ambient comparison. This “difference of differences” is also known as an interaction term, which we can formulate into a contrast with our combined count matrix.
con.combined.diff <- makeContrasts(
(groupKO.spinal - groupWT.spinal) - (groupKO.ambient - groupWT.ambient),
levels=design.combined
)
v.combined.diff <- voomLmFit(y.combined, design.combined, sample.weights=TRUE, contrasts=con.combined.diff)
fit.combined.diff <- treat(v.combined.diff, robust=TRUE, lfc=0.5)
res.combined.diff <- topTreat(fit.combined.diff, sort.by="none", n=Inf)
head(res.combined.diff)## ENSEMBL SYMBOL logFC AveExpr t
## ENSMUSG00000025902 ENSMUSG00000025902 Sox17 1.07509732 0.9683625 0.6326995
## ENSMUSG00000033845 ENSMUSG00000033845 Mrpl15 0.16394271 7.8057481 0.0000000
## ENSMUSG00000025903 ENSMUSG00000025903 Lypla1 -0.17385802 5.7232183 0.0000000
## ENSMUSG00000033813 ENSMUSG00000033813 Tcea1 -0.04511493 6.2962425 0.0000000
## ENSMUSG00000033793 ENSMUSG00000033793 Atp6v1h -0.15443631 4.9033011 0.0000000
## ENSMUSG00000025907 ENSMUSG00000025907 Rb1cc1 -0.29690973 3.4989041 0.0000000
## P.Value adj.P.Val
## ENSMUSG00000025902 0.3102854 1
## ENSMUSG00000033845 0.9999983 1
## ENSMUSG00000025903 0.9954541 1
## ENSMUSG00000033813 0.9999886 1
## ENSMUSG00000033793 0.9851374 1
## ENSMUSG00000025907 0.7822975 1A significant interaction means that the the spinal cord and ambient log-fold changes are significantly different. However, we want genes where the spinal log-fold change is more extreme than the ambient log-fold change. The interaction term itself represents the difference in the spinal cord and ambient log-fold changes, so we prioritize genes where the sign of the interaction is the same as the spinal cord log-fold change. (This also retains genes where the spinal cord and ambient log-fold changes are of different sign, which is fine.) We also refine the \(p\)-value to ensure that the gene is significant in the standalone spinal cord comparison, to avoid detecting irrelevant genes where there is a non-significant change in the spinal cord but a significant change in the opposite direction in the ambient solution.
# Putting relevant statistics in the same table for easier inspection.
res.combined2 <- res.combined.diff
res.combined2$spinal.logFC <- res.combined.spinal$logFC
res.combined2$spinal.P.Value <- res.combined.spinal$P.Value
res.combined2$ambient.logFC <- res.combined.ambient$logFC
res.combined2$ambient.P.Value <- res.combined.ambient$P.Value
# Recall that 'res.combined.diff$logFC' is the value of the interaction term.
# Setting uncooperative genes to have a p-value of 1 forces them to the bottom of the ranking.
res.combined2$P.Value[sign(res.combined2$logFC) != sign(res.combined2$spinal.logFC)] <- 1
# A gene must be significant in both the difference-of-differences contrast and
# in the standalone spinal cord comparison to achieve a low p-value here.
res.combined2$P.Value <- pmax(res.combined2$P.Value, res.combined2$spinal.P.Value)
res.combined2$adj.P.Val <- p.adjust(res.combined2$P.Value, method="BH")
head(res.combined2[order(res.combined2$P.Value),])## ENSEMBL SYMBOL logFC AveExpr t
## ENSMUSG00000086503 ENSMUSG00000086503 Xist -2.359372 6.099066 -9.118323
## ENSMUSG00000052187 ENSMUSG00000052187 Hbb-y -3.899549 7.886251 -6.272608
## ENSMUSG00000069919 ENSMUSG00000069919 Hba-a1 -3.645244 6.338655 -5.504155
## ENSMUSG00000055609 ENSMUSG00000055609 Hba-x -2.672150 8.673200 -6.428622
## ENSMUSG00000069917 ENSMUSG00000069917 Hba-a2 -3.316159 5.417905 -4.125483
## ENSMUSG00000037664 ENSMUSG00000037664 Cdkn1c -1.875127 4.483854 -3.313040
## P.Value adj.P.Val spinal.logFC spinal.P.Value
## ENSMUSG00000086503 6.596811e-12 7.675389e-08 -6.688018 2.401462e-33
## ENSMUSG00000052187 2.273993e-07 1.322895e-03 -9.114164 3.039058e-17
## ENSMUSG00000069919 9.608587e-07 3.726530e-03 -8.254954 1.794992e-17
## ENSMUSG00000055609 2.803563e-06 8.154865e-03 -7.446747 5.924820e-14
## ENSMUSG00000069917 8.389170e-05 1.952160e-01 -7.783435 5.309619e-14
## ENSMUSG00000037664 9.397119e-04 1.000000e+00 -5.506987 4.095503e-16
## ambient.logFC ambient.P.Value
## ENSMUSG00000086503 -4.328646 1.916941e-35
## ENSMUSG00000052187 -5.214615 6.588314e-47
## ENSMUSG00000069919 -4.609710 1.420022e-53
## ENSMUSG00000055609 -4.774597 9.754411e-27
## ENSMUSG00000069917 -4.467276 3.602415e-46
## ENSMUSG00000037664 -3.631860 9.075532e-32And after all that… we end up with our hemoglobins dominating the list again39, albeit with much weaker \(p\)-values. This highlights the main weakness of the interaction model - namely, its stronger assumptions regarding the behavior of the ambient contamination. Indeed, the spinal cord log-fold changes for the hemoglobins are more extreme than their ambient values, which violates our expectations for spurious ambient-derived DE. Our guess is that the actual contamination in each droplet is not accurately quantified by our ambient counts, e.g., due to enrichment of schizocytes in low-count barcodes. The interaction approach is much more sensitive to such inaccuracies than the removal approach, which is why we prefer the latter.
Even if the ambient profiles were estimated accurately, we would still recommend inspecting DE genes from both approaches. This is due to a particular quirk of using interaction models on count data whereby significant interactions will not be detected for genes with low expression in any comparison. For example, consider a gene with strong DE in the spinal cord comparison but no expression in either of the WT or KO ambient profiles. This gene should be prioritized but it will not have a significant interaction term as an ambient log-fold change cannot be reliably estimated. In contrast, the removal approach has no such problems detecting this gene due to the lack of ambient DE. By using both methods, we can use each method’s strengths to cover the other’s false negatives.
10.5 Further comments
Neither of the two approaches described above will detect a gene that is consistently DE in all subpopulations. In this case, the DE within each subpopulation is indistinguishable from the effect of ambient contamination.
Alternatively, we could attempt to subtract the ambient contamination from the count matrix itself, e.g., using tools like cellbender or SoupX.
I don’t have much to say on this topic other than to point out that
(i) estimating the contaminating proportion in each cell/pseudo-bulk sample is hard, and
(ii) subtracting the contamination without interfering with downstream statistical analyses (i.e., preserving the mean-variance relationship, null distribution of \(p\)-values) is hard40.
So, your mileage may vary.
In our opinion, the best solution to ambient contamination lies in careful experimental design. Specifically, multiple samples from different conditions should be multiplexed together in each run of a high-throughput scRNA-seq protocol. This ensures that samples from all conditions are subjected to the same ambient contamination, which should cancel out in the subsequent DE analysis.
Session information
## 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] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] edgeR_4.99.12 limma_3.99.4
## [3] scrapper_1.7.11 MouseGastrulationData_1.27.0
## [5] SpatialExperiment_1.23.0 SingleCellExperiment_1.35.2
## [7] SummarizedExperiment_1.43.0 Biobase_2.73.2
## [9] GenomicRanges_1.65.4 Seqinfo_1.3.2
## [11] IRanges_2.47.5 S4Vectors_0.51.10
## [13] BiocGenerics_0.59.12 generics_0.1.4
## [15] MatrixGenerics_1.25.0 matrixStats_1.5.0
## [17] BiocStyle_2.41.0
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.1 dplyr_1.2.1 blob_1.3.0
## [4] filelock_1.0.3 Biostrings_2.81.9 fastmap_1.2.0
## [7] BumpyMatrix_1.21.0 BiocFileCache_3.3.1 digest_0.6.39
## [10] lifecycle_1.0.5 statmod_1.5.2 KEGGREST_1.53.6
## [13] RSQLite_3.53.3 magrittr_2.0.5 compiler_4.6.1
## [16] rlang_1.3.0 sass_0.4.10 tools_4.6.1
## [19] yaml_2.3.12 knitr_1.52 S4Arrays_1.13.2
## [22] bit_4.6.0 curl_8.1.0 DelayedArray_0.39.8
## [25] abind_1.4-8 withr_3.0.3 purrr_1.2.2
## [28] grid_4.6.1 ExperimentHub_3.3.2 beachmat_2.29.3
## [31] cli_3.6.6 rmarkdown_2.32 crayon_1.5.3
## [34] otel_0.2.0 httr_1.4.9 rjson_0.2.23
## [37] BiocBaseUtils_1.15.1 DBI_1.3.0 cachem_1.1.0
## [40] parallel_4.6.1 AnnotationDbi_1.75.2 BiocManager_1.30.27
## [43] XVector_0.53.0 vctrs_0.7.3 Matrix_1.7-6
## [46] jsonlite_2.0.0 bookdown_0.48 BiocNeighbors_2.7.3
## [49] bit64_4.8.6 magick_2.9.1 jquerylib_0.1.4
## [52] glue_1.8.1 BiocVersion_3.24.0 tibble_3.3.1
## [55] pillar_1.11.1 rappdirs_0.3.4 htmltools_0.5.9
## [58] R6_2.6.1 dbplyr_2.6.0 httr2_1.3.0
## [61] evaluate_1.0.5 lattice_0.23-1 AnnotationHub_4.3.3
## [64] png_0.1-9 memoise_2.0.1 bslib_0.12.0
## [67] Rcpp_1.1.2 SparseArray_1.13.4 xfun_0.62
## [70] pkgconfig_2.0.3References
Sigh.↩︎
This reminds me of the old problem of input subtraction in ChIP-seq data, which I tried to solve during my PhD. In the end, we just decided to filter out regions with high input coverage rather than attempting subtraction. It was a lose-lose situation; either the input counts were small and subtraction had no meaningful effect, or the input counts were large and subtraction broke the expected mean-variance relationship in our count data.↩︎