用Scrublet推断单细胞实验中的doublet

基于液体的单细胞捕获技术的理想情况是一个液滴捕获一个细胞，但通常都会有捕获到双细胞甚至多细胞的情况。一般来说捕获的细胞数越多，概率上越可能遇到两个或多个细胞被包装进同一个droplet中，比如下表所示（来自：http://grcf.jhmi.edu/wp-content/uploads/2016/08/Chromium_SingleCell3.pdf）。

Screen Shot 2021-01-10 at 10.38.18 PM.png

Scrublet是个基于python的预测双细胞的工具，原理如文章中示意图所示。首先定义两种doublet，植入型（ embedded）和新表型（ neotypic），前者通常是同一细胞类型双细胞，后者是不同细胞类型双细胞，Scrublet算法只用于寻找后者。算法随机抽取成对的barcodes模拟出双细胞加到原表达矩阵中，并对所有细胞（包括模拟出的双细胞）进行聚类。最后根据聚类结果对细胞进行打分（doublet score），可以从示意图上看出细胞的分值和其有关联的模拟双细胞数目成正比，且分值越高越可能是真正的新表型双细胞。

Screen Shot 2021-01-10 at 11.06.32 PM.png

接下来用示例数据跑一下流程。首先加载scrublet及其他需要的python modules。

%matplotlib inline
import scrublet as scr
import scipy.io
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import os

读入示例数据，是一个有8381个细胞的矩阵。

input_dir = './test_data/pbmc8k_10X/filtered_gene_bc_matrices/GRCh38/'
counts_matrix = scipy.io.mmread(input_dir + '/matrix.mtx').T.tocsc()
genes = np.array(scr.load_genes(input_dir + 'genes.tsv', delimiter='\t', column=1))

print('Counts matrix shape: {} rows, {} columns'.format(counts_matrix.shape[0], counts_matrix.shape[1]))
print('Number of genes in gene list: {}'.format(len(genes)))

Counts matrix shape: 8381 rows, 33694 columns
Number of genes in gene list: 33694

建立scrublet对象并计算分数（doublet score）。这里设置doublet rate为6%，参考数值见文章开头的表。

scrub = scr.Scrublet(counts_matrix, expected_doublet_rate=0.06)

doublet_scores, predicted_doublets = scrub.scrub_doublets(min_counts=2, 
                                                          min_cells=3, 
                                                          min_gene_variability_pctl=85, 
                                                          n_prin_comps=30)

scrub.plot_histogram();

Preprocessing...
Simulating doublets...
Embedding transcriptomes using PCA...
Calculating doublet scores...
Automatically set threshold at doublet score = 0.23
Detected doublet rate = 4.4%
Estimated detectable doublet fraction = 61.0%
Overall doublet rate:
Expected = 6.0%
Estimated = 7.2%
Elapsed time: 9.4 seconds

image.png

上面左边的直方图显示实际细胞的分数，右边直方图显示模拟的双细胞的分数。右图中可以看出模拟的双细胞的分数成二项分布，左边的峰代表植入型双细胞分数较低，右边的峰代表新表型双细胞分数较高。利用右图的分布选定分数的阈值，就可以对实际的细胞进行双细胞的标注了。

看一眼效果。看起来有两个cluster是doublet。

def get_tsne(X, angle=0.5, perplexity=30, random_state=0, verbose=False):
    from sklearn.manifold import TSNE
    return TSNE(angle=angle, perplexity=perplexity, random_state=random_state, verbose=verbose).fit_transform(X)

scrub.set_embedding('TSNE', scr.get_tsne(scrub.manifold_obs_, 0.5,10))
scrub.plot_embedding('TSNE', order_points=True);

tsne_v2.png

把结果输出，然后就可以整合到下游分析中去啦。

barcodes = np.array(pd.read_csv(input_dir + '/barcodes.tsv', header = None, index_col=None))
a= np.array([barcodes[:,0], doublet_scores, predicted_doublets])
data = pd.DataFrame({'barcodes': a[0, :], 'score': a[1, :], 'prediction': a[2, :]})
data.to_csv('./doublet_prediction.txt', index=False,header=True)

Reference:
https://www.cell.com/cell-systems/pdf/S2405-4712(18)30474-5.pdf
https://github.com/AllonKleinLab/scrublet