[QC] DoubletFinder: 2024 업데이트 버전

이전포스트 및 관련 페이지

Mechanism & Processing, option관련 설명 이전 포스트 참조
👉DoubletFinder: doublet 제거법 01

🔗github page
🔗publication link / Published: April 03, 2019 / DOI: https://doi.org/10.1016/j.cels.2019.03.003

반년전쯤 DoubletFinder관련된 사용법을 정리할때까지만해도 Seurat v.5를 지원하지 않아서 따로 seurat v.3로 변환후 사용하는 방법을 사용했었는데, 얼마전에 (23년 11월) DoubletFinder 업데이트에서 Seurat v.5 호환이 가능하게 변경됐다! (Thanks for your update!)

이 업뎃에서 ‘_v3’ function들이 다 제거 돼서 나도 코드 업뎃하는겸 다시 포스팅하게 됨. + 이전에는 log-normalization을 사용했었는데 이번에는 single-cell 전용 툴인 SCTransform을 사용해서 진행해보려함.

Installation

remotes::install_github('chris-mcginnis-ucsf/DoubletFinder')

## for SCTransforom
BiocManager::install("glmGamPoi")

Usage

• Seurat object

  있어야겠지, 난 QC이전에 CellRanger 결과 읽어들인 raw data(paste0('00.',i,'_raw.rds'))를 준비했다. Seurat v.5

• sample 별 예상 dublet rate

  이전 포스트 참고해서 시퀀싱이후 얻은 cell 갯수 대비 예상되는 doublet rate을 계상해서 준비 line 55

• paramSweep, doubletFinder

  이전에는 뒤에 다 _v3붙었었는데 빼줌

• paramSweep(sct = TRUE)

  sct normalization 사용하고 싶으면 parameter설정해주기

library(Seurat)
options(Seurat.object.assay.version = "v5")
library(ggplot2)
library(patchwork)
library(tidyverse)

setwd('/your/working/dir')

SampleList <- c('S1','S2')

##############################################################################
## Doublet Finder
##############################################################################
library(DoubletFinder)

# load raw dataset
for(i in SampleList){
  assign(i,readRDS(paste0('00.',i,'_raw.rds')))
}

for(i in SampleList){
  file_name = paste0("DoubletFinder_",i,"_Doublet_cells.rds")
  print(paste("Checking for file:", file_name))  # Debugging line
  print(paste("File exists:", file.exists(file_name)))  # Debugging line
  

  if(file.exists(file_name)){
    cat(paste("The files ", file_name, "exists. Skipping to next loop. \n"))
    next
  }
  ## Pre-process Seurat object (sctransform) -----------------------------------------------------------------------------------
  seu <- get(i)
  
  ## Pre-process Seurat object (sctransform)
  seu <- SCTransform(seu)
  seu <- RunPCA(seu)
  seu <- RunUMAP(seu, dims = 1:10)
  seu <- FindNeighbors(seu, dims=1:10)
  seu <- FindClusters(seu)
  
  ## pK Identification (no ground-truth) ---------------------------------------------------------------------------------------
  sweep.res.list <- paramSweep(seu, PCs = 1:10, sct = TRUE)
  sweep.stats <- summarizeSweep(sweep.res.list, GT = FALSE)
  bcmvn <- find.pK(sweep.stats)

  ggplot(bcmvn, aes(pK, BCmetric, group=1)) + geom_point() + geom_line()
  ggsave(paste0("DoubletFinder_",i,"_peak.png"))

  pK <- bcmvn %>% filter(BCmetric==max(BCmetric)) %>% select(pK)
  pK <- as.numeric(as.character(pK[[1]])) 
  print(paste0(i,":",pK))

  ## Homotypic Doublet Proportion Estimate -------------------------------------------------------------------------------------
  annotations <- seu$seurat_clusters    
  multiRate= c('S1'=0.039,'S2'=0.054)

  homotypic.prop <- modelHomotypic(annotations)           ## ex: annotations <- seu_kidney@meta.data$ClusteringResults
  nExp_poi <- round(multiRate[i]*ncol(seu))  ## 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 ----------------------------------------------------------------
  seu <- doubletFinder(seu, PCs = 1:10, pN = 0.25, pK = pK, nExp = nExp_poi.adj, reuse.pANN = FALSE, sct = TRUE)

  DF.name = colnames(seu@meta.data)[grepl("DF.classification", colnames(seu@meta.data))]

  DimPlot(seu, group.by = DF.name) + NoAxes()
  ggsave(paste0("DoubletFinder_",i,".png"))

  VlnPlot(seu, features = "nFeature_RNA", group.by = DF.name, pt.size = 0.1)
  ggsave(paste0("DoubletFinder_",i,"_vln.png"),width=14)

  saveRDS(colnames(seu[, seu@meta.data[, DF.name] == "Singlet"]),paste0("DoubletFinder_",i,"_Singlet_cells.rds"))
  saveRDS(colnames(seu[, seu@meta.data[, DF.name] == "Doublet"]),paste0("DoubletFinder_",i,"_Doublet_cells.rds"))
}

끝