方法一：R包DoubletFinder
目的：在实际情况中，单细胞RNA测序时有一些液滴会包含超过两个或两个以上细胞，被称之为doublet或multiplet。这种情况会影响后续聚类分亚群的过程中的准确性，所以需要使用一些方法去除doublets。R包DoubletFinder是各种高分文献里常用的方法，并且容易上手，所以今天尝试用它进行doublet检测。
内容： 1.学习doublet检测的原理。
2.网上扒DoubletFinder包的使用脚本，进行对比学习。
3. 下载符合要求的脚本并进行调试，修改为自己可用的脚本，并运行。
数据：GEO: GSE144024获取的文件，孕后四周卵黄囊单细胞测序结果（文献Decoding Human Megakaryocyte Development）。
脚本来源：参考公众号（TPO生物信息）
步骤：在Rstudio中运行下面脚本。
##单细胞分析实录(4): doublet检测
install.packages('devtools')
library(devtools)
devtools::install_github('chris-mcginnis-ucsf/DoubletFinder')
library("scde")
library("DEsingle")
library("Seurat")
library("dplyr")
library("magrittr")
library("SingleCellExperiment")
library("patchwork")
library(ggplot2)
library(DoubletFinder)
test=read.csv("E:/YS_Stage_gene_counts.csv",header = T,row.names = 1)
test.seu <- CreateSeuratObject(counts = test)
test.seu <- NormalizeData(test.seu, normalization.method = "LogNormalize", scale.factor = 10000)
test.seu <- FindVariableFeatures(test.seu, selection.method = "vst", nfeatures = 2000)
test.seu <- ScaleData(test.seu, features = rownames(test.seu))
test.seu <- RunPCA(test.seu, features = VariableFeatures(test.seu),npcs = 50)
DimHeatmap(test.seu, dims = 1:20, cells = 500, balanced = TRUE)
test.seu <- FindNeighbors(test.seu, dims = 1:15)
test.seu <- FindClusters(test.seu, resolution = 0.5)
test.seu <- RunUMAP(test.seu, dims = 1:15)
test.seu <- RunTSNE(test.seu, dims = 1:15)
sweep.res.list <- paramSweep_v3(test.seu, PCs = 1:10, sct = FALSE)
for(i in 1:length(sweep.res.list)){
if(length(sweep.res.list[[i]]$pANN[is.nan(sweep.res.list[[i]]$pANN)]) != 0){
if(i != 1){
sweep.res.list[[i]] <- sweep.res.list[[i - 1]]
}else{
sweep.res.list[[i]] <- sweep.res.list[[i + 1]]
}
}
}
sweep.stats <- summarizeSweep(sweep.res.list, GT = FALSE)
bcmvn <- find.pK(sweep.stats)
pk_v <- as.numeric(as.character(bcmvn$pK))
pk_good <- pk_v[bcmvn$BCmetric==max(bcmvn$BCmetric)]
nExp_poi <- round(0.1*length(colnames(test.seu)))
#指定期望的doublet数
test.seu <- doubletFinder_v3(test.seu, PCs = 1:10, pN = 0.25, pK = pk_good, nExp = nExp_poi, reuse.pANN = FALSE, sct = FALSE)
#这一行是主要代码，会在test.seu@meta.data数据框的基础上加上两列
colnames(test.seu@meta.data)[ncol(test.seu@meta.data)]="DoubletFinder"
#第二列换一个列名
head(test.seu@meta.data)
test.seu@meta.data$CB=rownames(test.seu@meta.data)
DF_df <- test.seu@meta.data[,c("CB","DoubletFinder")]
write.table(DF_df, file = "E:/YS/DoubletFinder_result.txt", quote = FALSE, sep = '\t', row.names = FALSE, col.names = TRUE)
DimPlot(test.seu,reduction = "umap",pt.size = 2,group.by = "DoubletFinder")
结果：CB为cellular barcode，doubletFinder为每个细胞的预测结果。Singlet表示该液滴中为单个细胞，doublet表示该液滴中为两个或两个以上细胞。
image.png
分析与讨论：DoubletFinder的使用应该是针对每一个sample进行单独分析，并在用seurat对数据进行任何处理之前进行。同时因为DoubletFinder非常占用内存，我们需要及时保留我们需要的信息，即每个细胞的预测结果以及CB，并把不需要的部分删除。得到这样的结果文件后，我们在后续seurat的过程中就可以选择只去除Doublet细胞，或者去除doublet的影响较大的细胞亚群。
方法二：python中scrublet包
目的：scrublet是python中的一个包，用于doublet检测，是文献（Decoding Human Megakaryocyte Development）中使用的方法。所以今天学习这个包来复现文章结果。最重要的一点是scrublet比DoubletFinder要准确得多。DoubletFinder最终的结果受设定的pANN阈值的影响，结果稳定性不高。
内容：1.学习回顾python的使用方法。
2. 在网上扒Scrublet 包的使用方法，找到合适的，修改调试为自己可用的脚本，并运行。
数据：GEO: GSE144024获取的文件，孕后四周卵黄囊单细胞测序结果（文献Decoding Human Megakaryocyte Development）。
脚本：参考公众号（TPO生物信息）
步骤：1.在python运行以下脚本
import scrublet as scr
import scipy.io
import matplotlib.pyplot as plt
import numpy as np
import os
import pandas as pd
infile = "E:/NCBI/YS_Stage_gene_counts.csv"
outfile = "E:/NCBI/Scrublet_result.txt"
finallist = []
with open(infile, 'r') as f:
header = next(f)
cell_barcodes = header.rstrip().split(',')
for line in f:
tmpline = line.rstrip().split(',')[1: ]
tmplist = [int(s) for s in tmpline]
finallist.append(tmplist)
finalarray = np.array(finallist)
print("finalarray的值为：")
print(finalarray)
count_matrix = np.transpose(finalarray)
print('Count matrix shape: {} rows, {} columns'.format(count_matrix.shape[0], count_matrix.shape[1]))
scrub = scr.Scrublet(count_matrix, expected_doublet_rate = 0.1)
doublet_scores, predicted_doublets = scrub.scrub_doublets()
predicted_doublets_final = scrub.call_doublets(threshold = 0.2)
with open(outfile, 'w') as f:
f.write('\t'.join(['CB', 'Scrublet', 'Scrublet_Score']) + '\n')
for i in range(len(doublet_scores)):
if predicted_doublets_final[i] == 0:
result = 'Singlet'
else:
result = 'Doublet'
f.write('\t'.join([cell_barcodes[i], result, str(doublet_scores[i])]) + '\n')
2. 得到Scrublet_result.txt文件后，在RStudio中运行以下脚本进行可视化。
##python包scrublet结果可视化
library("scde")
library("DEsingle")
library("Seurat")
library("dplyr")
library("magrittr")
library("SingleCellExperiment")
library("patchwork")
library(ggplot2)
library(tidyverse)
test=read.csv("E:/NCBI/YS_Stage_gene_counts.csv",header=TRUE,row.names = 1)
test.seu <- CreateSeuratObject(counts = test, project = "Yolk sac", min.cells = 3, min.features = 200)
head(test.seu@meta.data)
test.seu <- NormalizeData(test.seu, normalization.method = "LogNormalize", scale.factor = 10000)
test.seu <- FindVariableFeatures(test.seu, selection.method = "vst", nfeatures = 2000)
test.seu <- ScaleData(test.seu, features = rownames(test.seu))
test.seu <- RunPCA(test.seu, features = VariableFeatures(test.seu),npcs = 50)
test.seu <- FindNeighbors(test.seu, dims = 1:20)
test.seu <- FindClusters(test.seu, resolution = 0.5)
test.seu <- RunUMAP(test.seu, dims = 1:20)
test.seu <- RunTSNE(test.seu, dims = 1:20)
test.seu@meta.data$CB=rownames(test.seu@meta.data)
SR_df=read.table("E:/NCBI/Scrublet_result.txt",header = T,sep = "\t",stringsAsFactors = F)
test.seu@meta.data=inner_join(test.seu@meta.data,SR_df,by="CB")
rownames(test.seu@meta.data)=test.seu@meta.data$CB
DimPlot(test.seu,reduction = "umap",pt.size = 2,group.by = "Scrublet")
graph2ppt(file="E:/NCBI/Scrublet检测doublet.ppt", width=4, aspectr=5)
image.png
分析与讨论：
使用scrublet，我们最终能获得了一个分数和一个基于分数得到的判断结果。从图C中可以看出doublet分布均一，并且与之前聚类的结果比较来看，各个亚群中doulet分布比例相近，因此判断doublet的检测结果可靠，在后面的分析中可以移除被判定是doublet的细胞，就可以得到更准确的亚群分布。
网友评论