單細胞轉錄組數據三五年前（2018）的10x技術一般來說都是3000左右的細胞數量每個樣品，費用也是三萬多，差不多了10塊錢人民幣一個細胞，那個時候基本上無需擔心雙細胞問題。（下面是2018的價格）
但隨著10x技術的迭代，目前一般來說都可以大多8000個左右的細胞費用降低到1.7萬人民幣左右，相當于2塊錢一個細胞。這個時候就很多人喜歡要求細胞數量越多越好，我看到很多數據集都是1.5萬個細胞甚至3萬個細胞的數量每個10x單細胞轉錄組樣品，實在是讓人很無語。
這樣的話，我們就不得不對每次得到的10x單細胞轉錄組表達量矩陣進行雙細胞的去除了，有一個基于Python的Scrublet工具去除雙細胞速度非常快，值得推薦。
一般來說單細胞轉錄組測序數據走完cellranger的定量流程，每個樣品就會得到3個表達量矩陣文件（barcodes.tsv.gz，matrix.mtx.gz，genes.tsv.gz或者features.tsv.gz），然后就可以走seurat流程進行單細胞降維聚類分群。但是因為要插入Python的Scrublet工具去除雙細胞流程，就變得復雜了一點。
首先是Python的Scrublet工具安裝
很簡單的 pip install 即可，一般來說大家的服務器或者蘋果電腦都是有Python環境的，也有pip這個安裝器，就是需要注意一下選擇合適的鏡像：
pip install -i https://pypi.tuna./simple scrublet
然后寫一個腳本( run_scrublet.py )：
有意思的是這個腳本名字居然還有講究，之前命名是（scrublet.py）一直報錯，后來把腳本的名字修改為 run_scrublet.py :
import sys
print(sys.path)#查看python安裝位置
import os, sys
os.getcwd()
os.listdir(os.getcwd())
# 1.加載Scrublet需要的包
import scrublet as scr
import scipy.io
import matplotlib.pyplot as plt
import numpy as np
import os
import pandas as pd
# 1.讀取輸入文件（輸入文件為10X的標準腳本）
# input_dir = '/home/bakdata/x10/jmzeng/matrix'
input_dir = sys.argv[1]
counts_matrix = scipy.io.mmread(input_dir + '/matrix.mtx.gz').T.tocsc()
out_df = pd.read_csv(input_dir + '/barcodes.tsv.gz', header = None, index_col=None, names=['barcode'])
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)
# doublets占比
print (scrub.detected_doublet_rate_)
out_df['doublet_scores'] = doublet_scores
out_df['predicted_doublets'] = predicted_doublets
out_df.to_csv(input_dir + '/doublet.txt', index=False,header=True)
out_df.head()
如果你有一個cellranger的定量流程的outputs文件夾，每個樣品就會得到3個表達量矩陣文件（barcodes.tsv.gz，matrix.mtx.gz，genes.tsv.gz或者features.tsv.gz），很容易對單個樣品跑上面的流程。
python run_scrublet.py your_out_by_cellranger
但是我們通常情況下都不會僅僅是一個樣品。
批量運行這個腳本
很簡單的，比如我的授課項目（PRJNA768891-ccRCC）路徑下面有4個樣品：
$ ls -d /home/bakdata/x10/jmzeng/2022-PRJNA768891-ccRCC/rna/outputs/matrix/*
/home/bakdata/x10/jmzeng/2022-PRJNA768891-ccRCC/rna/outputs/matrix/SRR16213611
/home/bakdata/x10/jmzeng/2022-PRJNA768891-ccRCC/rna/outputs/matrix/SRR16213612
/home/bakdata/x10/jmzeng/2022-PRJNA768891-ccRCC/rna/outputs/matrix/SRR16213613
/home/bakdata/x10/jmzeng/2022-PRJNA768891-ccRCC/rna/outputs/matrix/SRR16213614
那么，我們就批量運行這個腳本，使用下面的代碼：
ls -d /home/bakdata/x10/jmzeng/2022-PRJNA768891-ccRCC/rna/outputs/matrix/* |while read id;do(nohup python run_scrublet.py $id & );done
一般來說單細胞轉錄組測序數據走完cellranger的定量流程，每個樣品就會得到3個表達量矩陣文件（barcodes.tsv.gz，matrix.mtx.gz，genes.tsv.gz或者features.tsv.gz），但是我這個腳本并不需要基因信息，僅僅是讀取表達量矩陣即可判斷細胞是否是雙細胞，然后讀取barcodes.tsv.gz后把細胞是否是雙細胞這樣的信息寫出來即可。
每個樣品對應的文件夾里面都會輸出一個文本文件 (doublet.txt )：
$ head doublet.txt
barcode,doublet_scores,predicted_doublets
AAACCCAAGCAATTCC-1,0.020104244229337306,False
AAACCCAAGCATCGAG-1,0.10305343511450382,False
AAACCCAAGCTGGTGA-1,0.05794460641399417,False
AAACCCAAGGTGGGTT-1,0.014997115939242453,False
AAACCCAAGTGTTGTC-1,0.05228981544771018,False
AAACCCACAAATGGTA-1,0.04285714285714286,False
AAACCCACAACGGCTC-1,0.007255723960012898,False
AAACCCACACAGTCAT-1,0.03284527107241788,False
AAACCCACACGCTGCA-1,0.015420424789060229,False
當然了，其實還可以出一個umap去可視化被Python的Scrublet工具預測的雙細胞：
需要注意的是Python的Scrublet工具并不是幫你去除雙細胞，僅僅是預測給你。后續我們在進行降維聚類分群的時候就需要讀入這個文本文件 (doublet.txt )，自己進行去除雙細胞操作。我的示例代碼是：
###### step1:導入數據 ######
library(data.table)
dir='matrix/'
samples=list.files( dir )
samples
sceList = lapply(samples,function(pro){
# pro=samples[1]
folder=file.path( dir ,pro)
print(pro)
print(folder)
print(list.files(folder))
sce=CreateSeuratObject(counts = Read10X(folder),
project = pro ,
min.cells = 5,
min.features = 500)
doub = fread(file.path( dir ,pro,'doublet.txt'))
print(table(doub$predicted_doublets))
rm_barcodes = doub$barcode[doub$predicted_doublets]
print(sce)
sce=sce[,!colnames(sce) %in% rm_barcodes]
print(sce)
return(sce)
})
names(sceList)
后面就是走走seurat流程進行單細胞降維聚類分群。這樣的基礎分析，也可以看基礎10講：
