Python社区  »  Python

批量使用基于Python的Scrublet工具去除双细胞

生信技能树 • 3 周前 • 93 次点击  

单细胞转录组数据三五年前(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.tsinghua.edu.cn/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讲:


Python社区是高质量的Python/Django开发社区
本文地址:http://www.python88.com/topic/147237
 
93 次点击  
分享到微博