DoubletFinder 去除双细胞

双细胞去除与Seurat的交互:在seurat标准流程进行到TSNE和UMAP降维之后,FindAllMarkers之前进行DoubletFinder操作。

确定是否剔除预测双细胞的一般步骤:

  1. 查看降维图上预测出来的双细胞位置。在降维图上,真实的双细胞更倾向于处于主群边缘游离出来
  2. 如果预测的双细胞单独聚成群,查看其marker基因是否具有明显的多种细胞marker。
  3. 不是我们主要研究对象的细胞群中的双细胞可以不去除。

注意的问题:https://www.jianshu.com/p/8de4153d6250

DoubletFinder是一个基于R的doublet检测工具,在它的GitHub主页上,作者给出了如何将它与seurat结合进行分析的代码。在实践过程中,小L发现两个问题。

一是作者给出的范例仅针对我们有一个样本数据的情况下。如果我们有多个样本时,在分析的哪一步进行doublet检测就变得复杂起来了。我们是在对数据进行基础的筛选以及去除批次效应(batch effect)之后呢,还是之前呢?是使用原始的reads数呢?还是使用标准化(normalization)之后的reads数呢?

二是使用DoubletFinder最终会产生一个非常大的seurat对象(object),但其中可能有很多信息是我们不需要的,如果不进行合理的信息提取,会占用大量的计算资源,在样本数量多的情况下,甚至会导致更严重的后果。小L习惯使用Rstudio server进行数据分析,在使用DoubletFinder就曾经因为要分析多个sample,数据太大,导致Rstudio server崩溃。几经调试才找到合适的处理的方法。

小L进行了不同的尝试,翻阅了一些资料,还参考了运用类似思路进行doublet检测的不同工具的使用方法,最后得出结论:应该使用DoubletFinder根据它的tutorial直接对多个样本进行依次分析,先判断出每个细胞是否为doublet,再对它们一起进行过滤(filtering)、整合(integration)等操作。因此,我们在DoubletFinder进行分析后,只需要保留它的结论以及pANN值添加到我们原始数据的metadate里,就能进行后续分析了,不需要保留整个DoubletFinder分析产生的seurat对象(object)。

由于这一块的内容比较多,关于scrublet算法的分析、DoubletFinder和scrublet的比较以及相关的代码范例就放到下一篇内容里写。

作者:小L的读博日常
链接:https://www.jianshu.com/p/8de4153d6250
来源:简书
著作权归作者所有。商业转载请联系作者获得授权,非商业转载请注明出处。

下面内容来自 小L的读博日常 同名公众号

小L用过DoubletFinder和scrublet,它俩各有特点:

1. scrublet是基于python的。尽管它有非常详尽的tutorial,对于没有python基础的同学,还是有些许挑战。而DoubletFinder则是基于R的,使用R对于很多做scRNA-seq生信分析的同学来说,还是比较自然。

2. 因为scrublt是基于python的,相较于DoubletFinder,它的速度也要快很多。一个10x的样本只需要几分钟到十几分钟就能跑完,而DoubletFinder则需要几十分钟到一个多小时。

3. 最重要的一点是scrublet比DoubletFinder要准确得多。DoubletFinder最终的结果受设定的pANN阈值的影响,结果稳定性不高。小L曾经有一次得到超过90%数据都是doublet的结果。而scrublet的结果则要稳定得多,也符合我们的预期。

总结完毕,下面就是具体的实操了。正如小L在上篇推送中分析的,DoubletFinder的使用应该是针对每一个sample进行单独分析,并在用seurat对数据进行任何处理之前进行。同时因为DoubletFinder非常占用内存,我们需要及时保留我们需要的信息,即每个细胞的预测结果以及pANN值,并把不需要的部分删除。其中大写字母的部分需要大家根据自己的需求进行替换。

library(tidyverse)
library(Seurat)
library(DoubletFinder)

sampleV = c('SAMPLE1', 'SAMPLE2', 'SAMPLE3')
dataL = list()
## create seurat object for all samples
for (sample in sampleV) {
  file_name = str_c('DATA_FOLDER', sample, '_cellranger/outs/filtered_feature_bc_matrix')
  seuratObj_counts <- Read10X(data.dir = file_name)
  seuratObj = CreateSeuratObject(counts = seuratObj_counts, project = str_c(sample, "1"))
  dataL[[sample]] <- seuratObj
}

for (idx in 1:length(dataL) {
  ## preprocessing
  data <- NormalizeData(dataL[[idx]])
  data <- ScaleData(data, verbose = FALSE)
  data <- FindVariableFeatures(data, verbose = FALSE)
  data <- RunPCA(data, npcs = 30, verbose = FALSE)
  data <- RunUMAP(data, reduction = "pca", dims = 1:30)
  data <- FindNeighbors(data, reduction = "pca", dims = 1:30)
  data <- FindClusters(data, resolution = 0.5)
  sweep.res.list <- paramSweep_v3(data, PCs = 1:30, sct = FALSE)
  sweep.stats <- summarizeSweep(sweep.res.list, GT = FALSE)
  bcmvn <- find.pK(sweep.stats)

  ## Homotypic Doublet Proportion Estimate -------------------------------------------------------------------------------------
  annotations <- data@meta.data$ClusteringResults
  homotypic.prop <- modelHomotypic(annotations)           ## ex: annotations <- sample14@meta.data$ClusteringResults
  nExp_poi <- round(0.075*nrow(data@meta.data))  ## 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 ----------------------------------------------------------------
  data <- doubletFinder_v3(data, PCs = 1:30, pN = 0.25, pK = 0.09, nExp = nExp_poi, reuse.pANN = FALSE, sct = FALSE)
  ## save results
  dataL[[idx]]$doubFind_res = data@meta.data %>% select(contains('DF.classifications'))
  dataL[[idx]]$doubFind_score = data@meta.data %>% select(contains('pANN'))
}

经过上述处理后,dataL就是由所有sample得到的seurat对象构成的list。其中每一个seurat对象的metadata里包含了pANN值(doubFind_score)以及DoubletFinder的预测结果(doubFind_score)。用dataL,我们就可以根据seurat的教程进行后续的filtering、integration等操作了。需要注意的是,如果sample的数量太多,我们可以通过改for循环里idx的范围,每次只处理3-5个sample并及时保存,来避免内存不够等问题。

scrublet的安装可以在这里找到(https://github.com/swolock/scrublet),使用指南则在这里(https://github.com/swolock/scrublet/blob/master/examples/scrublet_basics.ipynb)。需要注意的是scrublet使用的是解压后的matrix.mtx文件和genes.tsv文件,而cellranger产生的是压缩文件,我们需要先提前解压一下。

在使用指南中更改好matrix.mtx和genes.tsv文件对应的位置之后,我们就可以按照里面的步骤一步步跑下来,得到doublet_scores和predicted_doublets。因为后续小L是用seurat对数据进行处理,所以就将doublet_scores和predicted_doublets单独保存下来,再放进seurat对象的metadata里。其中scrublet_res是scrublet的判断结果而scrublet_score是scurblet得到判断所基于的分数。

## 保存doublet_scores和predicted_doublets (python)
d = {'score': doublet_scores, 'prediction': predicted_doublets}
df = pd.DataFrame(data=d)
df.to_csv(DATA_FOLDER2 + sample + '_doublet.csv', index=False)

## 导入seurat对象的metadata里 (R)
library(tidyverse)
library(Seurat)
library(DoubletFinder)

sampleV = c('SAMPLE1', 'SAMPLE2', 'SAMPLE3')
dataL = list()
## create seurat object for all samples
for (sample in sampleV) {
  file_name = str_c('DATA_FOLDER', sample, '_cellranger/outs/filtered_feature_bc_matrix')
  seuratObj_counts <- Read10X(data.dir = file_name)
  seuratObj = CreateSeuratObject(counts = seuratObj_counts, project = str_c(sample, "1"))
  doublet_info = read.csv(str_c('DATA_FOLDER2', sample, '_doublet.csv'))
  seuratObj$scrublet_res = doublet_info$prediction
  seuratObj$scrublet_score = doublet_info$score
  dataL[[sample]] <- seuratObj
}

不管是使用DoubletFinder还是scrublet,我们最终都获得了一个分数和一个基于分数得到的判断结果。之后的步骤就是大家看一看doublet的分布,和在每一个细胞类型中的比例,来判断doublet的检测结果是否可靠,以及是要移除整一个细胞群还是移除被判定是doublet的细胞了。
祝大家吃好喝好睡好,科研快乐~

DoubletFinder的可重复性问题

作为一个多细胞去除工具,DoubletFinder似乎很流行,另外有文章指出,这个工具相对于其他同类竞品工具具有优势,详见此文:

https://www.sciencedirect.com/science/article/pii/S2405471220304592

真的是不是最优秀的我不知道,我一贯对此类benchmark文章持怀疑态度。但似乎有一点比较肯定就是这个工具还可以。但是有一个问题非常的让我恼怒,就是没有重复性,同样的结果,多跑几次后你会发现,每次的结果都不一样。虽然这个差别在此一步是小差别,但是其影响确实很深远的,单从我的个人使用经验来看,DoubletFinder的不可重复直接导致你的最后的cell cluster形状变个不停。

也许你觉得这不是问题,只要保存一下结果下次继续就可以了,或者设施一下set.seed。保存结果,这虽然是个解决方案,但是在分析数据时,如果常需要调整其他的参数,进而必须要重新进行Doublet 去除,就会带来很大的困扰。如果你尝试过set.seed,你会发现根本不管用,那么问题在哪里呢?今天决定决绝这个问题,一劳永逸。

研究源代码后发现,问题的的关键是,该包的作者在多个函数你使用了sample函数,所以这就解释了问什么结果总是变,另外set.seed之所以没用,也是因为这个原因,作者前后在不同的自定义函数里使用了多个含有sample的函数,且每次使用都没有设定seed,在我看来,这应该算是个bug,虽然该工具只是预测Doublets,但是也不代表你可以没有重复性吧?

解决方案:

如果你想让自己的结果可重复,可以在以下函数加入set.seed(),这些函数分别是paramSweep_v3,parallel_paramSweep_v3,summarizeSweep。最好加在sample函数的前面。实测加了之后,结果就非常的可重复了。希望能帮你省点时间。

作者:POLLIGATOR
链接:https://www.jianshu.com/p/f08df3f8860e
来源:简书
著作权归作者所有。商业转载请联系作者获得授权,非商业转载请注明出处。

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 159,716评论 4 364
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 67,558评论 1 294
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 109,431评论 0 244
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 44,127评论 0 209
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 52,511评论 3 287
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 40,692评论 1 222
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 31,915评论 2 313
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 30,664评论 0 202
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 34,412评论 1 246
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 30,616评论 2 245
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 32,105评论 1 260
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 28,424评论 2 254
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 33,098评论 3 238
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 26,096评论 0 8
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,869评论 0 197
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 35,748评论 2 276
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 35,641评论 2 271

推荐阅读更多精彩内容