单细胞代码解析-妇科癌症单细胞转录组及染色质可及性分析3

###########################################################
# Part 2: Doublet detection and removal
###########################################################

#' # Doublet detection: 1) DoubletDecon 2) DoubletFinder 3) Take intersection of calls
######################################################################################

# # Doublet Decon
# # Add if else statement to regress out nCount RNA if needed before running DoubletDecon
##DoubletDecon：该方法结合反卷积分析（deconvolution）和鉴定细胞独有的基因表达来检测Doublet
#软件链接：[https://github.com/EDePasquale/DoubletDecon](https://github.com/EDePasquale/DoubletDecon)
library(DoubletDecon)
#length() 函数用于获取或设置向量(列表)或其他对象的长度；levels专门是处理factor变量的level属性用的；as.numeric：将因子变量（factor）转化为数值变量
idents.length <- length(levels(Idents(rna)))
for (i in levels(Idents(rna))){
  i <- as.numeric(i)
  levels(Idents(rna))[i] <- i -1
}
#Improved_Seurat_Pre_Process()
#as.factor () R语言中的函数用于将传递的对象 (通常是Vector)转换为Factor。
#Idents(rna) <- as.factor(Idents(rna))
seuratObject=rna
#Seurat创建对象和细胞过滤
newFiles=Improved_Seurat_Pre_Process(seuratObject, num_genes=50, write_files=F)
#dev.off()
#rawDataFile：表达矩阵（gene by cell）；groupsFile：包含组分配的文件名（3列：cell，组，组（数字或字符））；filename：唯一的文件名，输入文件的名字；location：应在其中存储输出的目录
#fullDataFile：包含完整表达式数据的文件名（gene by cell）。默认值为NULL。removeCC：通过KEGG富集去除细胞周期基因。默认值为FALSE。
#species：科学物种名称，KEGG ID，三个字母的物种缩写或NCBI ID。默认值为“ mmu”。rhop：平均值x中的x * SD以确定黑名单中相关性的上限。默认值为1。
#write：将输出文件写为.txt文件。默认值为TRUE。recluster：recluster反卷积使用Hopach或反卷积分类分别对doublet和非doublet进行分类。
#PMF：在双重确定标准中使用步骤3（独特的基因表达）。默认值为TRUE。useFull：使用完整的基因列表进行PMF分析。需要fullDataFile。默认值为FALSE。
#heatmap：是否生成热图的布尔值。默认值为TRUE。大于约3000个像元的数据集可能比较慢。重心：在解卷积中，将重心用作参考，而不是默认重心。
#num_doubs：用户定义的每对集群要生成的双峰数目。默认值为100。only50：仅使用由50％/ 50％的父单元格混合创建的合成对偶，而不是30％/ 70％和70％/ 30％的扩展选项，默认为FALSE。
#min_uniq：挽救群集所需的最小独特基因数，默认值为4。
results=Main_Doublet_Decon(rawDataFile=newFiles$newExpressionFile,
                           groupsFile=newFiles$newGroupsFile,
                           filename=SAMPLE.ID,
                           location="./DoubletDecon",
                           fullDataFile=NULL,
                           removeCC=FALSE,
                           species="hsa",
                           rhop=1.1,
                           write=F,
                           PMF=TRUE,
                           useFull=FALSE,
                           heatmap=FALSE,
                           centroids=TRUE,
                           num_doubs=100,
                           only50=FALSE,
                           min_uniq=4,
                           nCores=4)
#rownames获取和设置数据框架的行名；gsub ()函数是2R语言中处理正则表达式中的一种
decon.doublets <- rownames(results$Final_doublets_groups)
decon.doublets <- gsub("\\.","-",decon.doublets)
# DoubletFinder
# Add modfieid Parameter sweep function to regress out nCOunt RNA if needed
# Add if else statement to regress out nCount RNA if needed before running DoubletFinder
#########################################################################################
#DoubletFinder是一种R包，可预测单细胞RNA 测序 数据中的doublet，具体解析[https://www.jianshu.com/p/b1947c4156ad](https://www.jianshu.com/p/b1947c4156ad)
#软件链接：[https://github.com/chris-mcginnis-ucsf/DoubletFinder](https://github.com/chris-mcginnis-ucsf/DoubletFinder)
library(DoubletFinder)
# Estimate Doublets
#################################################################################################
## pK Identification (no ground-truth) ---------------------------------------------------------------------------------------
#DoubletFinder 去除双细胞 
sweep.res.list <- paramSweep_v3(rna, PCs = 1:50, sct = FALSE)

sweep.stats <- summarizeSweep(sweep.res.list , GT = FALSE)
#找最佳 PK
bcmvn<- find.pK(sweep.stats)

pK.1 <- as.numeric(unfactor(dplyr::arrange(bcmvn,desc(BCmetric))$pK[1]))


nExp_poi <- round(doublet.rate*length(colnames(counts)))  ## Assuming 4.6% doublet formation rate - tailor for your dataset
# Run doubletFinder with pk.1
rna <- doubletFinder_v3(rna, PCs = 1:50, pN = 0.25, pK = pK.1, nExp = nExp_poi, reuse.pANN = FALSE, sct = FALSE)
meta <- rna@meta.data

#R 语言中的 paste0 () 函数用于连接所有没有分隔符的元素
doublet.column <- paste0("DF.classifications_0.25_",pK.1,"_",nExp_poi)


doublet.calls <- rna[[doublet.column]]
colnames(doublet.calls) <- "Call"
##dplyr::filter：数据过滤
rna.dub <- dplyr::filter(doublet.calls, Call == "Doublet")
rna.singlet <- dplyr::filter(doublet.calls, Call == "Singlet")

DF.doublets <- rownames(rna.dub)
# Intersect doublet calls
###########################################################################

head(DF.doublets)
head(decon.doublets)
#intersect () R语言中的函数用于查找两个对象的交集
doublets <- intersect(decon.doublets,DF.doublets)
rna@meta.data$Doublet.Call <- ifelse(rownames(rna@meta.data) %in% doublets,"TRUE","FALSE")
#FeatureScatte：在一组单个单元格中创建两个特征（通常是特征表达式）的散点图。细胞的颜色是由它们的身份类别决定的。皮尔逊两个特征之间的相关性显示在绘图上方。 
FeatureScatter(rna,feature1 = "nCount_RNA",feature2 = "nFeature_RNA",group.by = "Doublet.Call")+ggsave("Doublet_calls_FeatureScatter.png")
DimPlot(rna,group.by = "Doublet.Call")+ggsave("Doublet_calls_UMAP.png")
saveRDS(doublets,"Intersection_DoubletDecon_DoubletFinder_doublets.rds")
cells <- colnames(rna)