R/Find_double_jarning.R at main · wzwxh2023/R · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
PreprocessSeurat <- function(seu, PCs=1:10) {
  message("  Normalize ...")
  seu <- NormalizeData(seu, verbose = F)
  message("  Find variable genes ...")
  seu <- FindVariableFeatures(seu, selection.method = "vst", nfeatures = 2000, verbose = F)
  message("  Scale data ...")
  seu <- ScaleData(seu, verbose = F)
  message("  Run PCA ...")
  seu <- RunPCA(seu, verbose = F)
  message("  Find neighbors ...")
  k <- round(ncol(seu)*0.02)
  k <- ifelse(k < 20, 20, k)
  seu <- FindNeighbors(seu, dims = 1:10, reduction = "pca", k.param = k, verbose = F)
  message("  Find clusters ...")
  seu <- FindClusters(seu, resolution = 0.4, verbose = F)
  seu
}

FindOptimalpK <- function(seu, PCs=1:10, num.cores = 1) {
  sweep.res.list <- paramSweep_v3(seu, PCs = PCs, sct = F, num.cores = 1)
  sweep.stats <- summarizeSweep(sweep.res.list, GT = F)
  bcmvn <- find.pK(sweep.stats)
  pK <- bcmvn[which.max(bcmvn$BCmetric), ]$pK
  as.numeric(as.character(pK))
}


DF <- function(seu, PCs=1:10, auto.pK=TRUE, auto.cluster=TRUE) {
  message("1. Preprocessing ...")
  # this step is for getting clusters
  if (auto.cluster) {
    seu <- PreprocessSeurat(seu, PCs = PCs)
  }

  message("2. Find optimal pK ...")
  if (auto.pK) {
    optimal.pK <- FindOptimalpK(seu, PCs = PCs, num.cores = 1)
    message(paste0("   Optimal pK = ", optimal.pK))
  } else {
    optimal.pK <- 0.09 # default
  }

  message("3. Mark doublets ...")
  homotypic.prop <- modelHomotypic(seu$seurat_clusters)
  nExp_poi <- round(0.075*ncol(seu))
  nExp_poi.adj <- round(nExp_poi*(1-homotypic.prop))
  pN = 0.25 # default 0.25
  pANN_name <- paste("pANN", pN, optimal.pK, nExp_poi, sep = "_")
  DF_name <- paste("DF.classifications", pN, optimal.pK, nExp_poi, sep = "_")
  DF_name.adj <- paste("DF.classifications", pN, optimal.pK, nExp_poi.adj, sep = "_")

  seu <- doubletFinder_v3(seu, PCs = PCs, pN = pN, pK = optimal.pK, nExp = nExp_poi, reuse.pANN = F, sct = F)
  seu <- doubletFinder_v3(seu, PCs = PCs, pN = pN, pK = optimal.pK, nExp = nExp_poi.adj, reuse.pANN = pANN_name, sct = F)
  # summarize results
  seu$DF.classifications <- ifelse(seu@meta.data[[DF_name.adj]] == "Doublet", "Doublet.hc",
                                   ifelse(seu@meta.data[[DF_name]] == "Doublet", "Doublet.lc", "Singlet"))
  seu@meta.data[[pANN_name]] <- NULL
  seu@meta.data[[DF_name]] <- NULL
  seu@meta.data[[DF_name.adj]] <- NULL
  seu
}

MarkDoublets <- function(seu, PCs=1:10, split.by=NULL) {
  if (is.null(split.by)) {
    seu.list <- list(seu)
  } else {
    seu.list <- SplitObject(seu, split.by = split.by)
  }
  names(seu.list) <- NULL
  metadata.new <- pbapply::pblapply(seu.list, function(xx) {
    seu.tmp <- DF(xx, PCs=PCs, auto.pK=TRUE, auto.cluster=TRUE)
    seu.tmp@meta.data
  })
  metadata.new <- do.call(rbind, metadata.new)
  seu$DF.classifications <- metadata.new[rownames(seu@meta.data), ]$DF.classifications
  return(seu)
}