Hello.

I am reading through the "best practices for scRNA-seq without multiplexing" and would just like to confirm that for my experiment which used the Chromium i7 Multiplex chip, we had pairs of sample multiplexed and the pairs were then demultiplexed using CellRanger. So this section referring to distinct sample types across lanes does not apply to my multiplexed experiment correct?

For those using multiplexing we can refer to the "real-world applications" using the standard Seurat normalization correct?

Thank you very much,
AC

Hi,
I was running the function "modelHomotypic". however, I got a error as follow:

homotypic.prop <- modelHomotypic(ctrl_doublet@ident)
Error in modelHomotypic(ctrl_doublet@ident) : object 's' not found

I wonder that what is the mean of object 's', and how to debug it.

Hi there,

thanks for the great package and the really nice preprint. I really enjoyed reading the preprint and I think it is amazing that the implementation directly runs upon a Seurat object, since many people (myself included) work with Seurat directly.

I ran into an issue when using the doubletFinder() function on a Seurat object for which more than 20 PCs (the default value) were used to perform the tSNE.

I think I found the culprit for this, which is in line 66 and of doubletFinder.R:

Which will select the number of PCs that were used in the original tSNE clustering.
To get the PCA coordinates though, we try to subset the PCs calculated in seu_wdoublets using the number of PCs used in the original Seurat object. If this number is >20 PCs, the following line:

will return a subscript out of bounds error.

I created a pull request in which I added pcs.compute to the RunPCS command so that the number of PCs that will be calculated corresponds to the number of PCs initially used in the original Seurat object.

Let me know if there are any problems with this or if there is a better way to implement this.
Hope this helps,

Best,

Florian

Hi All,
I work on immune cells, as you know macrophages are phagocytic, they can eat other cells.
When I analyse macropahges, a cluster of them express markers of other immune cells, like T cells.
Do you think using, your tool can differentiate between the two? I mean between a doublet and between a macrophages that ate a T cell.
Thank you

Appreciate this is probably not all that helpful without me uploading the data, but I'm getting this error:

>     sweep.res.list_kidney <- DoubletFinder::paramSweep(seu_kidney)
[1] "Creating artificial doublets for pN = 5%"
[1] "Creating Seurat object..."
[1] "Normalizing Seurat object..."
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
[1] "Finding variable genes..."
Calculating gene means
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating gene variance to mean ratios
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
[1] "Scaling data..."
Scaling data matrix
  |=====================================================================================================================================| 100%
[1] "Running PCA..."
Error in irlba(A = t(x = data.use), nv = pcs.compute, ...) : 
  max(nu, nv) must be positive

Ran the standard steps prior to this:

    seu_kidney <- CreateSeuratObject(raw.data = bb,min.cells = min.cells, min.genes = min.genes)
    seu_kidney <- NormalizeData(seu_kidney)
    seu_kidney <- ScaleData(seu_kidney, vars.to.regress = "nUMI")
    seu_kidney <- FindVariableGenes(seu_kidney, x.low.cutoff = 0.0125, y.cutoff = 0.25, do.plot=FALSE)
    seu_kidney <- RunPCA(seu_kidney, pc.genes = [email protected], pcs.print = 0)
    seu_kidney <- RunTSNE(seu_kidney, dims.use = 1:10, verbose=TRUE)

Hi,

Thanks so much for developing this amazing, user-friendly tool. I was wondering if you could help me understand why each time I run pK Identification (no ground-truth), I get very different results. I saw from rdrr.io (https://rdrr.io/github/chris-mcginnis-ucsf/DoubletFinder/man/paramSweep_v3.html) that each time it is run, it will select 10,000 different cells (from my total of 35,000 cells), but I'm a bit surprised that the results vary so much. Attached here are the two most disparate pK plots so you can get a sense of the differences that I am seeing. Do you have any recommendation on how to select an appropriate pK value given this observation?

Thank you so much for any help; I greatly appreciate it.

Best,
Maria

bcmvn.AllIslets_9.pdf
bcmvn.AllIslets_2.pdf

Hi

I am currently following your "Example code for 'real-world' applications".
Does it always find a fixed number of doublets?

The example code first find the fixed number of doublets based on the "expected doublet rate" in

Exp_poi <- round(doubletRatio.expected * length(Cells(seurat)))
seurat <- doubletFinder(seurat, PCs = 1:10, pN = 0.25, pK = 0.09, nExp = nExp_poi, reuse.pANN = FALSE)
The numbero doublets seems to be always the same with the expected number from my own manual parameter, "doubletRatio.expected".

Then, it identifies high confident doublets among them by

nExp_poi.adj <- round(nExp_poi*(1-homotypic.prop))
seurat <- doubletFinder(seurat, PCs = 1:10, pN = 0.25, pK = 0.09, nExp = nExp_poi.adj, reuse.pANN = "pANN_0.25_0.09_913")

It seems like the expected doublet rate is a critical parameter. Can I get some advice about how to set it? or please let me know if my understanding is wrong.

Thank you.

Hi,

  I wanted to ask how one can calculate DE genes between tSNE clusters in a Seurat object upon removal of doublet cells,

Best,
Theo

Hi, DoubletFinder seems to be a very nice tool that I would like to apply on my own aligned Seurat dataset. Here, the Seurat object is pre-processed slightly different as different commands are used e.g. RunCCA instead of RunPCA etc (see also here: https://satijalab.org/seurat/immune_alignment.html).
Already with the first command in pK Identification I get an error:

[1] "Running PCA..."
Error in irlba(A = t(x = data.use), nv = pcs.compute, ...) :
max(nu, nv) must be positive

I assume it is due to differences in the commands. Is it possible to adapt the underlying functions of DoubletFinder so that it can be used for aligned Seurat objects? That would be really great!

I can't see a single-peak after running paramSweep_v3
the BCmetric is increasing as pK increases
which pK shall I choose?

The bcmvn i got are:
ParamID | pK | MeanBC | VarBC | BCmetric
1 | 0.005 | 0.6322526 | 0.0062336671 | 101.42546
2 | 0.01 | 0.6314945 | 0.0043372477 | 145.59797
3 | 0.02 | 0.6406362 | 0.0072198519 | 88.73260
4 | 0.03 | 0.6560862 | 0.0083522676 | 78.55187
5 | 0.04 | 0.6724027 | 0.0067474174 | 99.65335
6 | 0.05 | 0.6764357 | 0.0066137246 | 102.27757
7 | 0.06 | 0.6784079 | 0.0074017078 | 91.65560
8 | 0.07 | 0.6799124 | 0.0072759030 | 93.44714
9 | 0.08 | 0.6813135 | 0.0072258698 | 94.28809
10 | 0.09 | 0.6811553 | 0.0079572110 | 85.60226
11 | 0.1 | 0.6798985 | 0.0090639511 | 75.01127
12 | 0.11 | 0.6894440 | 0.0084931059 | 81.17690
13 | 0.12 | 0.6904313 | 0.0102618693 | 67.28124
14 | 0.13 | 0.6963331 | 0.0093759757 | 74.26780
15 | 0.14 | 0.6995890 | 0.0097656529 | 71.63771
16 | 0.15 | 0.7064218 | 0.0089377553 | 79.03794
17 | 0.16 | 0.7099229 | 0.0089569153 | 79.25975
18 | 0.17 | 0.7130803 | 0.0084972736 | 83.91872
19 | 0.18 | 0.7191425 | 0.0078496397 | 91.61471
20 | 0.19 | 0.7274438 | 0.0065485259 | 111.08512
21 | 0.2 | 0.7332291 | 0.0052203877 | 140.45491
22 | 0.21 | 0.7420419 | 0.0039294137 | 188.84290
23 | 0.22 | 0.7478681 | 0.0031535904 | 237.14815
24 | 0.23 | 0.7499830 | 0.0028707247 | 261.25215
25 | 0.24 | 0.7517944 | 0.0028600436 | 262.86118
26 | 0.25 | 0.7542440 | 0.0020662547 | 365.02952
27 | 0.26 | 0.7596140 | 0.0014461647 | 525.26108
28 | 0.27 | 0.7642677 | 0.0009825212 | 777.86381
29 | 0.28 | 0.7678856 | 0.0006758615 | 1136.15822
30 | 0.29 | 0.7678427 | 0.0005725801 | 1341.02242
31 | 0.3 | 0.7684884 | 0.0004597986 | 1671.35856

Thanks for your help!

Hello, I have the same question as this post(#44 (comment)) .I got the final picture ,but I want to know why find.pk alway show NULL when i do pK identification with and without ground-truth? This is my code and I provide the information you need.

# Do not apply DoubletFinder to aggregated scRNA-seq data representing multiple distinct samples
SeuratOb <- base::subset(Merge.integrated,time.sum %in% "P7")
DefaultAssay(SeuratOb) <- "RNA"
SeuratOb@reductions
SeuratOb <- NormalizeData(SeuratOb)
SeuratOb <- ScaleData(SeuratOb)
SeuratOb <- FindVariableFeatures(SeuratOb, selection.method = "vst", nfeatures = 2000)
SeuratOb <- RunPCA(SeuratOb)
ElbowPlot(SeuratOb,ndims = 50)
SeuratOb <- RunUMAP(SeuratOb, dims = 1:20)
## pK Identification (no ground-truth) ---------------------------------------------------------------------------------------
sweep.res.list_kidney <- paramSweep_v3(SeuratOb, PCs = 1:20, sct = FALSE)
sweep.stats_kidney <- summarizeSweep(sweep.res.list_kidney, GT = FALSE)
bcmvn_kidney <- find.pK(sweep.stats_kidney)
## pK Identification (ground-truth) ------------------------------------------------------------------------------------------
sweep.res.list_kidney <- paramSweep_v3(SeuratOb, PCs = 1:20, sct = FALSE)
gt.calls <- [email protected][rownames(sweep.res.list_kidney[[1]]), "GT"]
sweep.stats_kidney <- summarizeSweep(sweep.res.list_kidney, GT = TRUE, GT.calls = gt.calls)
bcmvn_kidney <- find.pK(sweep.stats_kidney)

## Homotypic Doublet Proportion Estimate -------------------------------------------------------------------------------------
annotations <- [email protected]$CellType.sum1
homotypic.prop <- modelHomotypic(annotations)           ## ex: annotations <- [email protected]$ClusteringResults
nExp_poi <- round(0.075*(nrow(SeuratOb)))  ## Assuming 7.5% doublet formation rate - tailor for your dataset
nExp_poi.adj <- round(nExp_poi*(1-homotypic.prop))

## Run DoubletFinder with varying classification stringencies ----------------------------------------------------------------
SeuratOb <- doubletFinder_v3(SeuratOb, PCs = 1:20, pN = 0.20, pK = 0.01, nExp = nExp_poi, reuse.pANN = FALSE , sct = FALSE)
head(SeuratOb[[]])
SeuratOb <- doubletFinder_v3(SeuratOb, PCs = 1:20, pN = 0.20, pK = 0.01, nExp = nExp_poi.adj, reuse.pANN = "pANN_0.2_0.01_2166", sct = FALSE)
a <- DimPlot(SeuratOb,pt.size = 1,label=TRUE, label.size = 5,reduction = "umap",group.by = "DF.classifications_0.2_0.01_2166" )+theme(aspect.ratio = 1)
a

> sweep.stats_kidney <- summarizeSweep(sweep.res.list_kidney, GT = FALSE)
> bcmvn_kidney <- find.pK(sweep.stats_kidney)
NULL
> sweep.stats_kidney <- summarizeSweep(sweep.res.list_kidney, GT = TRUE, GT.calls = gt.calls)
Error in model.matrix.default(mt, mf, contrasts) : 
  variable 1 has no levels
> bcmvn_kidney <- find.pK(sweep.stats_kidney)
NULL
> head(sweep.res.list_kidney[[1]])
                         pANN
P7_AAAGTGACAAGACCGA 0.0000000
P7_AACCCAACAGATCCTA 0.0000000
P7_AACCTTTGTGAATTGA 0.0212766
P7_AACGTCAGTGAGATCG 0.0000000
P7_AACTTCTGTGCATACT 0.0000000
P7_AAGCCATCACCGTGCA 0.0212766
> head(sweep.stats_kidney)
    pN    pK    BCreal
1 0.05 0.005 0.8992749
2 0.05  0.01 0.9082518
3 0.05  0.02 0.9088909
4 0.05  0.03 0.9113203
5 0.05  0.04 0.9013240
6 0.05  0.05 0.8913376
> head(bcmvn_kidney)
  ParamID    pK    MeanBC        VarBC BCmetric
1       1 0.005 0.8815419 0.0002655315 3319.915
2       2  0.01 0.8719980 0.0006510888 1339.292
3       3  0.02 0.8608932 0.0006136421 1402.924
4       4  0.03 0.8845970 0.0002277028 3884.875
5       5  0.04 0.8488235 0.0007731505 1097.876
6       6  0.05 0.8483289 0.0006561837 1292.822

Dear all,

Do we select pk corresponding to the highest BCmetric in the variable bcmvn_kidney?

Thank you

It's interesting that 14% of cells in each of 5 sample are decided as doublets. Is it due to some coincidence or due to settings in DoubletFinder? Or may be some of my code is wrong. Could you please tell me?

The cell numbers before and after subset (by "percent.mt <= Percent.mt", "DF_hi.lo == "Singlet"", or both) is listed below:

the original code is:

library(Seurat)
library(DoubletFinder)
Seurat.list <- c('GTE001', 'GTE002', 'GTE003', 'GTE009', 'GTE012')
for (i in 1:length(Seurat.list)) {
  Filename.brief <- paste('20191122_', Seurat.list[i], sep = '')
  data_dir <- paste('D:/Cell Ranger Processed Data/', Seurat.list[i], '/filtered_feature_bc_matrix', sep = '')
  Seurat.data <- Read10X(data.dir = data_dir)
  dense.size <- object.size(x = as.matrix(x = Seurat.data))
  Seurat.object <- CreateSeuratObject(counts = Seurat.data, project = Seurat.list[i], min.cells = 5)
  Seurat.object[["percent.mt"]] <- PercentageFeatureSet(Seurat.object, pattern = "^MT-")
  Seurat.object <- NormalizeData(Seurat.object, verbose = FALSE)
  Seurat.object <- FindVariableFeatures(Seurat.object, selection.method = "vst", nfeatures = 2000)
  Seurat.object <- ScaleData(Seurat.object, verbose = FALSE)
  Seurat.object <- RunPCA(Seurat.object, npcs = 30, verbose = FALSE, resolution = 0.1)
  Seurat.object <- RunTSNE(Seurat.object, dims = 1:10, resolution = 0.1)
  Seurat.object <- FindNeighbors(Seurat.object, reduction = "tsne", dims = 1:2)
  Seurat.object <- FindClusters(Seurat.object, resolution = 0.1)
  #Doublet Finder
  sweep.res.list_Seurat <- paramSweep_v3(Seurat.object, PCs = 1:10, sct=FALSE)
  sweep.stats_Seurat <- summarizeSweep(sweep.res.list_Seurat, GT = FALSE)
  bcmvn_Seurat <- find.pK(sweep.stats_Seurat)
  mpK<-as.numeric(as.vector(bcmvn_Seurat$pK[which.max(bcmvn_Seurat$BCmetric)]))
  annotations <- [email protected]$seurat_clusters
  homotypic.prop <- modelHomotypic(annotations)           ## ex: annotations <- [email protected]$ClusteringResults
  nExp_poi <- round(0.075*length(Seurat.object$seurat_clusters))  ## Assuming 7.5% doublet formation rate - tailor for your dataset
  nExp_poi.adj <- round(nExp_poi*(1-homotypic.prop))
  Seurat.object <- doubletFinder_v3(Seurat.object, PCs = 1:20, pN = 0.25, pK = mpK, nExp = nExp_poi, reuse.pANN = FALSE, sct = FALSE)
  Seurat.object <- doubletFinder_v3(Seurat.object, PCs = 1:10, pN = 0.25, pK = 0.09, nExp = nExp_poi.adj, reuse.pANN =colnames([email protected][1,])[13], sct = FALSE)
  [email protected][,'DF_hi.lo'] <- "Singlet"
  Code_of_Temp <- paste("[email protected]$DF_hi.lo[which([email protected]$", colnames(as.matrix([email protected][1,]))[13], " == 'Doublet')] <- 'Doublet_lo'", sep = '')
  eval(parse(text = Code_of_Temp))
  Code_of_Temp <- paste("[email protected]$DF_hi.lo[which([email protected]$", colnames(as.matrix([email protected][1,]))[14], " == 'Doublet')] <- 'Doublet_hi'", sep = '')
  eval(parse(text = Code_of_Temp))
  #Subset: MT%<=12% & ==singlet
  Percent.mt <- 12
  Seurat.object_temp <- subset(Seurat.object, subset = percent.mt <= Percent.mt & DF_hi.lo == "Singlet" )
}

hi, when useparamSweep_v3, it reports error:
Error in NormalizeData.default(object = GetAssayData(object = object, : trying to get slot "params" from an object of a basic class ("NULL") with no slots
but when I set sct argument to TRUE, it just worked fine, the question is I didn't use SCTransform to normalize my data, is this a normal thing to see?

Hi,
I'm actually testing multiple tools for scRNAseq analysis.
My actual pipeline is working with sctransform instead of the standard Seurat workflow. I've found a way to make it work with your tool.
But do you plan to allow the user to choose the normalization method?
Thanks,
JonathanCruard

Hey,

thanks for developing doubletFinder. I find it really useful.

The structure of Seurat objects changed with Seurat 3.0 and your function is currently incompatible with it. It returns "no slot of name "data" for this object of class "Seurat"".

Best Max

Hi Chris,

I'm excited to try DoubletFinder on our data!
Going through the paper, I realized I have a (perhaps too naive) question about the doublet rate correction, which I can't wrap my head around. I'm hoping you can help clarify the logic behind it:

Estimating Homotypic Doublet Proportions
In this study, homotypic doublet proportions were modeled as the sum of squared cell state frequencies. For example, consider a scRNA-seq dataset with five unique cell states present at the following proportions:

pCi = {0.40, 0.25, 0.15, 0.1, 0.1}

Where pCi is the proportion of cell state i.

The proportion of homotypic doublets present in this data is then estimated as:

The final number of detectable (i.e., heterotypic) doublets is then defined by adjusting the total number of doublets (i.e., as determined by the Poisson doublet formation rate) by the homotypic doublet proportion:

DDR = pHomo*TDR
Where DDR and TDR are the detectable and total doublet rates, respectively.

If pHomo is the proportion of homotypic (and undetectable) doublets, then why isn't the detectable doublet rate the complement of that proportion? (i.e. 1 - pHomo ) ?

Thanks a lot!

Hi all,
I'm excited to try DoubletFinder. I'm running into a silly error, which has me stalled in the very first step.
I have a Seurat object that has been through all the necessary steps (NormalizeData, FindVariableGenes, ScaleData, RunPCA, and RunTSNE). I also successfully installed DoubletFinder and all the dependencies. However, when I run doubletFinder_ParamSweep, I get an error saying:
could not find function "doubletFinder_ParamSweep.

Here is my session info:

> sessionInfo()
R version 3.5.2 (2018-12-20)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Sierra 10.12.6

Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods   base    

other attached packages:
 [1] ROCR_1.0-7          gplots_3.0.1        modes_0.7.0         KernSmooth_2.23-15  fields_9.6          maps_3.3.0        
 [7] spam_2.2-1          dotCall64_1.0-0     DoubletFinder_2.0.0 Seurat_2.3.4        Matrix_1.2-15       cowplot_0.9.4     
[13] ggplot2_3.1.0     

loaded via a namespace (and not attached):
  [1] snow_0.4-3             backports_1.1.3        Hmisc_4.1-1            VGAM_1.0-6             plyr_1.8.4           
  [6] igraph_1.2.2           lazyeval_0.2.1         splines_3.5.2          densityClust_0.3       fastICA_1.2-1        
 [11] digest_0.6.18          foreach_1.4.4          htmltools_0.3.6        viridis_0.5.1          lars_1.2             
 [16] gdata_2.18.0           magrittr_1.5           checkmate_1.9.1        cluster_2.0.7-1        mixtools_1.1.0       
 [21] limma_3.36.5           matrixStats_0.54.0     R.utils_2.7.0          docopt_0.6.1           colorspace_1.4-0     
 [26] ggrepel_0.8.0          xfun_0.4               dplyr_0.7.8            sparsesvd_0.1-4        crayon_1.3.4         
 [31] jsonlite_1.6           bindr_0.1.1            survival_2.43-3        zoo_1.8-4              iterators_1.0.10     
 [36] ape_5.2                glue_1.3.0             gtable_0.2.0           kernlab_0.9-27         prabclus_2.2-6       
 [41] BiocGenerics_0.26.0    DEoptimR_1.0-8         scales_1.0.0           pheatmap_1.0.12        mvtnorm_1.0-8        
 [46] bibtex_0.4.2           Rcpp_1.0.0             metap_1.0              dtw_1.20-1             viridisLite_0.3.0    
 [51] htmlTable_1.13.1       reticulate_1.10        foreign_0.8-71         bit_1.1-14             proxy_0.4-22         
 [56] mclust_5.4.2           SDMTools_1.1-221       Formula_1.2-3          tsne_0.1-3             stats4_3.5.2         
 [61] htmlwidgets_1.3        httr_1.4.0             FNN_1.1.2.2            RColorBrewer_1.1-2     fpc_2.1-11.1         
 [66] acepack_1.4.1          modeltools_0.2-22      ica_1.0-2              pkgconfig_2.0.2        R.methodsS3_1.7.1    
 [71] flexmix_2.3-14         nnet_7.3-12            tidyselect_0.2.5       rlang_0.3.1            reshape2_1.4.3       
 [76] munsell_0.5.0          tools_3.5.2            ggridges_0.5.1         stringr_1.3.1          npsurv_0.4-0         
 [81] knitr_1.21             bit64_0.9-7            fitdistrplus_1.0-11    robustbase_0.93-3      caTools_1.17.1.1     
 [86] DDRTree_0.1.5          purrr_0.2.5            RANN_2.6.1             bindrcpp_0.2.2         pbapply_1.3-4        
 [91] nlme_3.1-137           monocle_2.8.0          slam_0.1-44            R.oo_1.22.0            hdf5r_1.0.1          
 [96] compiler_3.5.2         rstudioapi_0.9.0       png_0.1-7              lsei_1.2-0             tibble_2.0.1         
[101] stringi_1.2.4          lattice_0.20-38        trimcluster_0.1-2.1    HSMMSingleCell_0.114.0 pillar_1.3.1         
[106] combinat_0.0-8         Rdpack_0.10-1          lmtest_0.9-36          data.table_1.12.0      bitops_1.0-6         
[111] irlba_2.3.2            gbRd_0.4-11            R6_2.3.0               latticeExtra_0.6-28    gridExtra_2.3        
[116] codetools_0.2-16       MASS_7.3-51.1          gtools_3.8.1           assertthat_0.2.0       withr_2.1.2          
[121] qlcMatrix_0.9.7        diptest_0.75-7         parallel_3.5.2         doSNOW_1.0.16          rpart_4.1-13         
[126] tidyr_0.8.2            class_7.3-15           segmented_0.5-3.0      Rtsne_0.15             Biobase_2.40.0       
[131] base64enc_0.1-3      
> ## pK Identification ---------------------------------------------------------------------------------------------------------
> sweep.res.list_Mixing <- doubletFinder_ParamSweep(seu_Mixing)
Error in doubletFinder_ParamSweep(seu_Mixing) :
  could not find function "doubletFinder_ParamSweep"

Again, apologies if there is something minor I'm missing. But does anyone know why this is happening?

Hi,
Thank you for your work, I'm testing it right now ! I'm working on a huge 10X single cell dataset and it seems doublets will be an issue to caracterize cell populations..
But I was wondering about one thing with your approach : in my dataset, whatever clustering approach I use, it obtain what may be a doublet cluster, certainly containing cells from the two known major populations in my samples. Don't you think it might interfere with doubletFinder as it will create synthetic doublets with doublets from the cluster and other valid cells ?
Thanks !

Hi,

Thanks for the tool - I am looking forward to trying it. I am running into an issue early on. I have a Seurat2 object. When I run param_sweep(seuratObj, PCs=10), I get the following error:

Traceback:

1. paramSweep(novaS2, PCs = 10)
2. RunPCA(seu_wdoublets, pc.genes = [email protected], pcs.print = 0, 
 .     pcs.compute = length(PCs))

Any idea where the problem could lie?

Thanks in advance!

Hi,
Thanks for developing such a useful tool! I'm using DoubletFinder to identify doublets from my 10X data, but it throws an error when I run summarizeSweep. Below is my code:

infile <- args[1]
outfile <- args[2]
indata <- read.table(infile, sep='\t', header = T, row.names = 1)
seu <- CreateSeuratObject(indata)
seu <- NormalizeData(seu, normalization.method = 'LogNormalize', scale.factor = 100000)
seu <- ScaleData(seu, vars.to.regress = 'nUMI')
seu <- FindVariableGenes(seu, x.low.cutoff = 0.0125, y.cutoff = 0.5, do.plot = FALSE)
seu <- RunPCA(seu, pc.genes = [email protected], pcs.print = 0, pcs.compute = 30)
seu <- FindClusters(seu, dims.use = 1:30, resolution = 0.5, print.output = FALSE)
seu <- RunTSNE(seu, dims.use = 1:30, do.fast = TRUE)

Above is code for some regular Seurat preprocessing, there is no error and I got a Seurat object "seu". Then I run doubletfinder:

sweep_res_list <- paramSweep(seu)
sweep_stats <- summarizeSweep(sweep_res_list, GT = FALSE)
bcmvn <- find.pK(sweep_stats)
pk <- as.numeric(as.vector(bcmvn$pK)[which.max(bcmvn$BCmetric)])
nExp_poi <- round(0.075 * length([email protected]))
seu <- doubletFinder(seu, pN = 0.25, pK = pk, nExp = nExp_poi, reuse.pANN = FALSE)

The function paramSweep work well, but when I run summarizeSweep, there is an error:
Error in seq.default(a, b, length = M) : 'from' must be a finite number
I upload the sweep_res_list (an rds file) here:
sweep_res_list.zip
I found that this script can run successfully on some samples, but failed on the samples with fewer cells (1000-2000). Could you figure out this problem? Thanks a lot!

Yang

Hi,

I was running DoubletFinder with a Seurat object. I just got an error message as follows;

[1] "Creating artificial doublets..."
[1] "Creating Seurat object..."
[1] "Normalizing Seurat object..."
[1] "Finding variable genes..."
Calculating gene means
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating gene variance to mean ratios
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Error in seq.int(rx[1L], rx[2L], length.out = nb) :
'to' must be a finite number

Just a question - I'm working with a highly heterogeneous population and haven't quite decided the right level of resolution. How dependent is the % doublet detection on the number of annotated clusters (i.e. if you chose a low resolution and the dataset contains 20 clusters, will you find a drastically different doublet percentage if you increase resolution and find 40)?

Hi, I just started using your package and it's working great. There's a branch in the Seurat repo that has support for loom objects, so we don't have to load the entire dataset into memory. Right now, Seurat supports normalization, scaling, finding variable genes, PCA, and tSNE for loom objects, so the Seurat pipeline for loom is very similar to that for sparse matrices in memory. Perhaps loom will become more and more common as we're getting more large datasets. However, your package only works for Seurat objects, because AddMetaData in Seurat doesn't have a method for loom objects yet; I don't know if a method for loom will be implemented in the near future. The way to add metadata into loom objects is obj$add.col.attribute, see this vignette. Can you write a method of your function for loom objects? Or may I do a pull request?

Hi Chris,
I tested DoubletFinder_v3 on 1 sample and got promising result. Now I am trying to use DoubletFinder_v3 for my combined object (40 samples, size ~12GB). I ran it on our cluster with 200GB memory and my job was killed by R all the time due to too much memory.
Here is the code:
seu <- doubletFinder_v3(seu, PCs = 1:10, pN = 0.25, pK = pk, nExp = nExp_poi, reuse.pANN = FALSE, sct = T)
Do you have any suggestions for a better approach?
Thanks,
Linh

Hi there,

Thank you for sharing this. I tried to use the code provided to analysis my own data and got just one doublet. I am not sure if I have done something wrong but I am not confident about the result.
Let me describe my scRNAseq data before hand:
About 9000 cells from a type of gland were processed in 10x. Since the tissue is not easy to dissociate, it is expected to have doublets in single-cell samples.
I process the scRNAseq data in Seurat v3 (no filtering for potential doublets), and then use the codes below to find doublets
sweep.res.list_Gland <- paramSweep_v3(Gland, PCs = 1:10) sweep.stats_Gland <- summarizeSweep(sweep.res.list_Gland, GT = FALSE) bcmvn_Gland <- find.pK(sweep.stats_Gland) annotations <- [email protected]$RNA_snn_res.0.6 homotypic.prop <- modelHomotypic(annotations) nExp_poi <- round(0.075*length(Gland)) nExp_poi.adj <- round(nExp_poi*(1-homotypic.prop)) Gland <- doubletFinder_v3(Gland, PCs = 1:10, pN = 0.25, pK = 0.09, nExp = nExp_poi, reuse.pANN = FALSE) Gland <- doubletFinder_v3(Gland, PCs = 1:10, pN = 0.25, pK = 0.09, nExp = nExp_poi.adj, reuse.pANN = "pANN_0.25_0.09_0") [email protected][,"DF_hi.lo"] <- [email protected]$DF.classifications_0.25_0.09_0 [email protected]$DF_hi.lo[which([email protected]$DF_hi.lo == "Doublet" & [email protected]$DF.classifications_0.25_0.09_0 == "Singlet")] <- "Doublet_lo" [email protected]$DF_hi.lo[which([email protected]$DF_hi.lo == "Doublet")] <- "Doublet_hi" TSNEPlot(Gland, group.by="DF_hi.lo", plot.order=c("Doublet_hi","Doublet_lo","Singlet"), colors.use=c("black","gold","red"))
Here is some output from the code above:

bcmvn_NonCF <- find.pK(sweep.stats_NonCF)
NULL
nExp_poi value: 0
nExp_poi.adj value: 0
Doublet found: 1

Thanks a lot!

Is it possible to run the method without Seurat? For example starting from a raw count or normalized gene expression matrix?

When running package functions R raised an error that packages fields and modes are missing. You can add them as dependencies in DESCRIPTION to make sure users install them before running the package.

Hi!

This looks like a great tool and I look forward to using it.

I was wondering, if I have multiple 10x runs that I am analyzing together, would you recommend identifying doublets from each 10x run separately, removing these cells and then merging all runs together to perform downstream analyses? Or do you think it would be okay to analyze all samples after a merge?

thanks!

I may solve the problem by myself... hence I left the "correct" code below, which was editted from the websitd. Thanks!

Hi, when I have the tSNE plot, there is no difference between "Doublet_hi" and "Doublet_lo", as the figure attached.

I am using the deveping version of Seurat3. Hence, I have to change the tutorial code from chris-mcginnis-ucsf/DoubletFinder to the followings. So the problem of "Doublet_hi" and "Doublet_lo" may be due to the change I induced. Coule you help me?

  sweep.res.list_Seurat <- paramSweep_v3(Seurat.object, PCs = 1:10, sct=FALSE)
  sweep.stats_Seurat <- summarizeSweep(sweep.res.list_Seurat, GT = FALSE)
  bcmvn_Seurat <- find.pK(sweep.stats_Seurat)
  mpK<-as.numeric(as.vector(bcmvn_Seurat$pK[which.max(bcmvn_Seurat$BCmetric)]))
  annotations <- [email protected]$seurat_clusters
  homotypic.prop <- modelHomotypic(annotations)           ## ex: annotations <- [email protected]$ClusteringResults
  #nExp_poi <- round(0.075*length([email protected]))  ## Assuming 7.5% doublet formation rate - tailor for your dataset
  nExp_poi <- round(0.075*length(Seurat.object$seurat_clusters))  ## Assuming 7.5% doublet formation rate - tailor for your dataset
  nExp_poi.adj <- round(nExp_poi*(1-homotypic.prop))
  Seurat.object <- doubletFinder_v3(Seurat.object, PCs = 1:20, pN = 0.25, pK = mpK, nExp = nExp_poi, reuse.pANN = FALSE, sct = FALSE)
  Seurat.object <- doubletFinder_v3(Seurat.object, PCs = 1:10, pN = 0.25, pK = 0.09, nExp = nExp_poi.adj, reuse.pANN =colnames(as.matrix([email protected][1,]))[9], sct = FALSE)
  png(paste(Filename.brief, '_Doublets_UMAP1.png', sep = ''), width = 500, height = 500)
  plot(DimPlot(Seurat.object, reduction = "umap", group.by = colnames(as.matrix([email protected][1,]))[9], label = TRUE))
  dev.off()
  png(paste(Filename.brief, '_Doublets_tSNE1.png', sep = ''), width = 500, height = 500)
  plot(DimPlot(Seurat.object, reduction = "tsne", group.by = colnames(as.matrix([email protected][1,]))[9], label = TRUE))
  dev.off()
  png(paste(Filename.brief, '_Doublets_UMAP2.png', sep = ''), width = 500, height = 500)
  plot(DimPlot(Seurat.object, reduction = "umap", group.by = colnames(as.matrix([email protected][1,]))[10], label = TRUE))
  dev.off()
  png(paste(Filename.brief, '_Doublets_tSNE2.png', sep = ''), width = 500, height = 500)
  plot(DimPlot(Seurat.object, reduction = "tsne", group.by = colnames(as.matrix([email protected][1,]))[10], label = TRUE))
  dev.off()
  [email protected][,'DF_hi.lo'] <- "Singlet"
  Code_of_Temp <- paste("[email protected]$DF_hi.lo[which([email protected]$", colnames(as.matrix([email protected][1,]))[10], " == 'Singlet' & [email protected]$", colnames(as.matrix([email protected][1,]))[9], " == 'Doublet')] <- 'Doublet_lo'", sep = '')
  eval(parse(text = Code_of_Temp))
  Code_of_Temp <- paste("[email protected]$DF_hi.lo[which([email protected]$", colnames(as.matrix([email protected][1,]))[10], " == 'Doublet')] <- 'Doublet_hi'", sep = '')
  eval(parse(text = Code_of_Temp))
  png(paste(Filename.brief, '_Doublets_UMAP3.png', sep = ''), width = 500, height = 500)
  plot(UMAPPlot(Patient.object, group.by="DF_hi.lo", plot.order=c("Doublet_hi","Doublet_lo","Singlet"), cols=c("red","orange","blue")))
  dev.off()
  png(paste(Filename.brief, '_Doublets_tSNE3.png', sep = ''), width = 500, height = 500)
  plot(TSNEPlot(Patient.object, group.by="DF_hi.lo", plot.order=c("Doublet_hi","Doublet_lo","Singlet"), cols=c("red","orange","blue")))
  dev.off()

Hi,
Thank you for developing this powerful tool to identify doublets. I used it for two single-cell RNA sequencing datasets with 9000 cells and the estimated percentage of doublets was 7.6%. A small part of identified doubles changes every time when I tried to re-run the script using the same pk value. Is that because the pseudodoublets are created randomly so that it won't provide a constant result for cells with a lower possibility of doublet? And does that means the results of each run is unreplicable? Is there any recommendation for this issue? Many thanks.

Given that you have found pN optimization is largely unnecessary, would it be possible to update the paramSweep function to allow for a static pN?

Hello! I am interested in using this package however i was wondering if there is an example workflow?

I am little confused in the supplementary information from the cell systems paper that states the doubledetection starts with pre-processed data, but the steps in the first section take the data through to PCA and tsne cluster detection in the seurat pipeline. please can you clarify?

Also, can you show an example of how to split the seurat object in to their different batches (i.e individual 10x lanes) and running doubletdetection on all individual batches simultaneously?

does the seurat pipeline have to be re-run after this?

many thanks!

Hi,
Thank you for your very nice tools.
Could you give me codes to visualize doublet cells on tsne plot when working on seu_kidney object?
Thank you!

what annotation file should be fill here?
homotypic.prop <- modelHomotypic(annotations)
I added '[email protected]' , is that correct?
And I noticed that the number of finally detected doublet is set to 0.075*cell number. Is 0.075 based on your experience?

Hi there!

I wanted to report 2 minor points/issues:

1. If you start from a Seurat object which was created using the SubsetData function, the raw.data slot contains the raw data from the full object - this means that the new Seurat object will contain everything which was present in the initial dataset.
2. The RunPCA function by default only computes 20 PCs. If the original Seurat object used more, this causes an error since the ones >20 are out of bounds in the new object containing doublets. This could also easily be adjusted by either automatically computing a larger number of PCs (pcs.compute) or having this parameter as an option.

I hope this makes sense. Thanks!!!

Gabriela

Dear All,
Thank you so much for this interesting paper and package.
I was wondering if there was any tutorial easy to follow for those that aren't very good at coding.
I'm using Seurat V3
Thank you.

I installed DoubletFinder today, but am not able to use the paramSweep_V3 and doubletFinder_V3 functions. Thanks for any help you can provide.

library(DoubletFinder)
sweep.res.list_MDS <- paramSweep_V3(MDS.list[[1]])
Error in paramSweep_V3(MDS.list[[1]]) :
could not find function "paramSweep_V3"
session_info()
─ Session info ──────────────────────────────────────────────────────────────────────────────────────────
setting value
version R version 3.5.2 (2018-12-20)
os macOS Mojave 10.14.3
system x86_64, darwin15.6.0
ui RStudio
language (EN)
collate en_US.UTF-8
ctype en_US.UTF-8
tz America/Chicago
date 2019-04-05
...
DoubletFinder * 2.0.0 2019-04-05 [1] Github (6643b82)

Hi,

Thanks for developing doubletFinder. It's really useful.

When I use "doubletFinder_ParamSweep" to identification pK, it shows errors.

An example of the bottom of the output file when errors are encountered is shown below.

sweep.res.list_pbmc4k <- doubletFinder_ParamSweep(pbmc4k_0.6)
[1] "Creating artificial doublets for pN = 5%"
[1] "Creating Seurat object..."
[1] "Normalizing Seurat object..."
Performing log-normalization
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
[1] "Finding variable genes..."
Calculating gene means
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating gene variance to mean ratios
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
[1] "Scaling data..."
Scaling data matrix
|======================================================================| 100%
[1] "Running PCA..."
[1] "Calculating PC distance matrix..."
[1] "Defining neighborhoods..."
[1] "Computing pANN across all pK..."
[1] "pK = 5e-04..."
Error in row.names<-.data.frame(*tmp*, value = value) :
invalid 'row.names' length

Best

I just updated to Seurat's latest version and latest version of DoubletFinder.

Getting the following error when I try to run paramSweep_v3:

[1] "Creating artificial doublets for pN = 5%" [1] "Creating Seurat object..." [1] "Normalizing Seurat object..." Error in NormalizeData.default(object = GetAssayData(object = object, : trying to get slot "params" from an object of a basic class ("NULL") with no slots

Hi there,

I have tested your package with Seurat V3 and it works great. I do have one naive question about the algorithm: seems like the artificial doublets are generated by averaging two random cells UMIs, whereas I would think of doublets as double/multiple UMIs, could you explain the assumption behind it (or does that limit the algorithm to a certain type of cell population)? Thanks.

Hi, I was using DoubletFinder with Seurat3, but when I estimate homotypic rate, it turns out that nExp_poi.adj value is negative, and I went to see the R code of modelHomotypic:
modelHomotypic <- function(annotations) { anno.freq <- table(annotations)/length(annotations) x <- sum(anno.freq^2) return(x) }
I wonder if length(annotations) also means the cell numbers? cause in Seurat3, length(object[["cluster"]]) is 1, if I change this to dim(object[["cluster"]])[1], nExp_poi.adj is a normal value between 0 and 1.

Thanks for the convenient tools to detect the doublets!
I met the problem when runing:
sweep.res.list_seurat <- paramSweep_v3(seurat, PCs = 1:75, sct = T)
Should this step take a long time to run? My desktop computer did not finish running it in 24h with the R console filled with the mesages of SCTransform.
Like this:
[1] "Creating artificial doublets for pN = 5%"
Loading required package: sctransform
[1] "Creating Seurat object..."
[1] "Running SCTransform..."
Calculating cell attributes for input UMI matrix
Variance stabilizing transformation of count matrix of size 15867 by 6415
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 2000 genes, 6415 cells
|==================================================================================================================================| 100%non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = 0.500000non-integer x = ...............

I am new to the scRNA-seq. I want to folllow the example and know more about DoubleFinder. How can I get the example data?
Thanks!

Hi,
I am trying to install this package using the following code:
devtools::install_github('chris-mcginnis-ucsf/DoubletFinder')
however, I get an error message like this:

Downloading GitHub repo chris-mcginnis-ucsf/DoubletFinder@master
Error: HTTP error 422.
No commit found for SHA: release/3.0

Rate limit remaining: 57/60
Rate limit reset at: 2019-06-21 11:00:20 UTC

Do you have any suggestions about how can I solve this problem? Thank you very much!

The function summarizeSweep cannot handle cases when all pANN values for a certain parameter set are 0. It gives the following error message:

Error in if (L == 0) warning("Binning grid too coarse for current (small) bandwidth: consider increasing 'gridsize'") :
missing value where TRUE/FALSE needed

This is because bkde function from KernSmooth cannot deal with singular input. You need to treat this as a special case.

Also, is it a sign of some kind of problem when all pANN are 0 for some parameters? In my case, all pN=0.3 give all-zero pANN, as well as some other combinations. This seems suspicious to me. If this is the case, it is worth notifying the user about this problem and possible causes.

Hi,
Thank you for developing / maintaining an excellent tool for single-cell biology.
I found this tool works well for scRNA-seq datasets.
It should theoretically also work for scATAC-seq datasets (peak-bc matrix) - What do you think?
Could you think it is possible to make the vignette when you were available or in the next version?
Thanks,
Yoshi

Hi,

Where is the list of doublets stored?

Thanks

Hi,I've read your latest paper in Cell System.And I have a question:

When using DoubletFinder,the user should input the doublets number which could be estimated by Poisson distribution according to the paper.But if I don't know the loading number of cells in the process of library building,how could I do? Is it possible to judge by the density plot of pANN ?
Thanks~