CeleScope — 新格元单细胞多组学分析工具箱

生信菜鸟团

发布于 2024-07-31 19:20:45

2500

发布于 2024-07-31 19:20:45

文章被收录于专栏：生信菜鸟团

工欲善其事必先利其器

CeleScope

CeleScope是一款由新格元生物科技有限公司自主开发的，用于处理新格元单细胞系列产品测序数据的开源生信软件。可从二代测序下机的原始fastq数据开始，经过细胞标签提取、质控与校正，参考基因组比对，完成基因定量，最终得到质控报告和表达矩阵，以用于单细胞下游分析。

GitHub：https://github.com/singleron-RD/CeleScope
官方文档：https://github.com/singleron-RD/CeleScope/blob/master/doc/user_guide.md

如何安装

使用conda 配置所需依赖软件
pip安装最新的celescope

##下载
wget https://gitee.com/singleron-rd/celescope/raw/master/conda_pkgs.txt

##配置环境
mamba create -n celescope -y --file conda_pkgs.txt

## pip安装
conda activate celescope
pip install celescope

## 测试是否安装成功
celescope -h

所需依赖 conda_pkgs.txt

按正常流程安装，由于pip安装不检查依赖，所以会自动安装的最新的的numpy2.0.1版本，这样的话。后面运行会出现以下报错：

运行测试报错

报错信息提示在加载 pandas 库的过程中，初始化 pandas 的 _libs.interval 模块时报错。错误信息 ValueError: numpy.dtype size changed, may indicate binary incompatibility 指出 numpy 的数据类型大小发生了变化，这可能表明存在兼容性问题。

由于pandas 库依赖于特定版本的 numpy 。这个报错明显是安装的pandas和 numpy 版本不匹配的问题。这个时候我们检查安装记录，由于pandas是要求固定版本 1.4.2 ，所以只能考虑降级numpy版本

版本要求

已安装版本检查

如果不确定需要哪个numpy版本，可以先尝试最低要求版本，然后再根据提示调整适配的版本

#先尝试最低要求版本，
pip install  numpy==1.21

#根据提示调整适配的版本
pip install  numpy==1.22.4

测试

调整

测试，安装成功

测试运行通过

功能模块

功能模块

最小化使用

我们还是以单细胞转录组数据定量为例

首先还是构建参考基因组文件

##激活环境
conda activate celescope
##过滤gtf
celescope utils mkgtf Homo_sapiens.GRCh38.99.gtf Homo_sapiens.GRCh38.99.filtered.gtf
##构建参考文件
celescope rna mkref \
 --genome_name Homo_sapiens_ensembl_99_filtered \
 --fasta Homo_sapiens.GRCh38.dna.primary_assembly.fa \
 --gtf Homo_sapiens.GRCh38.99.filtered.gtf \
 --mt_gene_list mt_gene_list.txt

--genome_name ##基因组名称
--fasta ##参考基因组fasta文件
--gtf ##参考基因组过滤后的gtf文件
--mt_gene_list ##线粒体基因列表名
--thread ##调用线程数。默认6

其次是根据数据运行对应的功能模块

用 multi_rna 为每个样本构建 celescope rna 分析的 shell 脚本

multi_rna\
 --mapfile ./rna.mapfile\
 --genomeDir {path to hs_ensembl_99 or mmu_ensembl_99}\
 --thread 8\
 --mod shell

--mapfile ##一个以制表符分隔的文本文件
--genomeDir ##参考基因组文件所在目录
--thread ##线程数，默认4
--mod ##指定生成脚本的类型，可选{sjm,shell}

mapfile是以制表符分隔的文本文件，至少3列，分别为

Fastq文件前缀
Fastq文件目录路径
样本名，也即所有输出文件的前缀

第四和第五列根据不同的分析，会有不同的要求。视具体分析而定。

一个实例

同样我们也找一个实例数据跑跑看。数据集来自文章《CD276-dependent efferocytosis by tumor-associated macrophages promotes immune evasion in bladder cancer》，发表在《Nature Communications》上。文章主要研究了在膀胱癌中，肿瘤相关巨噬细胞（TAMs）通过CD276依赖的凋亡细胞清除作用（efferocytosis）如何促进肿瘤免疫逃逸。

数据集是：CRA008674
链接：：https://ngdc.cncb.ac.cn/gsa/browse/CRA008674
数据下载方法见: 使用人工智能优化一个数据库文件批量下载脚本

准备参考文件(mouse)

##使用前先激活环境
conda activate celescope
##过滤gtf文件
celescope utils mkgtf ../Mus_musculus.GRCm39.111.gtf  Mus_musculus.GRCm39.111.filtered.gtf
##创建索引文件
celescope rna mkref \
 --genome_name Mus_musculus_ensembl_111_filtered \
 --fasta ../Mus_musculus.GRCm39.dna.toplevel.fa \
 --gtf ./Mus_musculus.GRCm39.111.filtered.gtf \
 --mt_gene_list mt_gene_list.txt

输出文件大小

  95 7月  22 16:58 mt_gene_list.txt
1.2G 7月  22 16:58 Mus_musculus.GRCm39.111.filtered.gtf

参考基因组文件

 272 7月  22 19:38 ./celescope_genome.config
 459 7月  22 19:10 ./chrLength.txt
 943 7月  22 19:10 ./chrNameLength.txt
 484 7月  22 19:10 ./chrName.txt
 667 7月  22 19:10 ./chrStart.txt
 28M 7月  22 19:10 ./exonGeTrInfo.tab
 12M 7月  22 19:10 ./exonInfo.tab
1.3M 7月  22 19:10 ./geneInfo.tab
2.7G 7月  22 19:38 ./Genome
 780 7月  22 19:38 ./genomeParameters.txt
 13K 7月  22 19:38 ./Log.out
  95 7月  22 16:58 ./mt_gene_list.txt
1.2G 7月  22 16:58 ./Mus_musculus.GRCm39.111.filtered.gtf
 21G 7月  22 19:38 ./SA
1.5G 7月  22 19:38 ./SAindex
7.9M 7月  22 19:34 ./sjdbInfo.txt
7.8M 7月  22 19:10 ./sjdbList.fromGTF.out.tab
6.2M 7月  22 19:34 ./sjdbList.out.tab
8.2M 7月  22 19:10 ./transcriptInfo.tab

制作mapfile文件

下载数据集信息 CRA008674.xlsx ，可以得到CRR编号和样本的对应信息，保存为crr_sample.txt

$cat crr_sample.txt 
CRR596006       Blca_WT1
CRR596007       Blca_WT2
CRR596008       Blca_wKO1
CRR596009       Blca_wKO2
CRR596010       Blca_WT3
CRR596011       Blca_WT4
CRR596012       Blca_cKO1
CRR596013       Blca_cKO2

## 根据信息制作mapfile文件
$cat crr_sample.txt |awk -v OFS="\t" '{print $1,"/home/yjzhang/celescope_test/st1_data",$2}'
CRR596006       /home/yjzhang/celescope_test/st1_data   Blca_WT1
CRR596007       /home/yjzhang/celescope_test/st1_data   Blca_WT2
CRR596008       /home/yjzhang/celescope_test/st1_data   Blca_wKO1
CRR596009       /home/yjzhang/celescope_test/st1_data   Blca_wKO2
CRR596010       /home/yjzhang/celescope_test/st1_data   Blca_WT3
CRR596011       /home/yjzhang/celescope_test/st1_data   Blca_WT4
CRR596012       /home/yjzhang/celescope_test/st1_data   Blca_cKO1
CRR596013       /home/yjzhang/celescope_test/st1_data   Blca_cKO2

cat crr_sample.txt |awk -v OFS="\t" '{print $1,"/home/yjzhang/celescope_test/st1_data",$2}' > cra008674.mapfile

为每个样本生构建 celescope rna 分析的 shell 脚本

multi_rna \
 --mapfile ./rna.mapfile \
 --genomeDir ~/reference/mouse/ensembl/refdata_celescope \
 --thread 8 \
 --mod shell

此时需要注意一下文原始fq文件的文件命名

由于GSA的数据下载下来文件名为：

27G 7月  15 00:31 CRR596006/CRR596006_f1.fq.gz
27G 7月  15 00:12 CRR596006/CRR596006_r2.fq.gz
41G 7月  14 18:46 CRR596007/CRR596007_f1.fq.gz
44G 7月  14 19:17 CRR596007/CRR596007_r2.fq.gz
33G 7月  14 22:29 CRR596008/CRR596008_f1.fq.gz
34G 7月  14 22:06 CRR596008/CRR596008_r2.fq.gz
42G 7月  14 20:34 CRR596009/CRR596009_f1.fq.gz
44G 7月  14 21:06 CRR596009/CRR596009_r2.fq.gz
32G 7月  14 19:40 CRR596010/CRR596010_f1.fq.gz
33G 7月  14 20:04 CRR596010/CRR596010_r2.fq.gz
28G 7月  14 23:53 CRR596011/CRR596011_f1.fq.gz
28G 7月  14 23:33 CRR596011/CRR596011_r2.fq.gz
29G 7月  14 23:12 CRR596012/CRR596012_f1.fq.gz
31G 7月  14 22:52 CRR596012/CRR596012_r2.fq.gz
24G 7月  14 21:23 CRR596013/CRR596013_f1.fq.gz
26G 7月  14 21:41 CRR596013/CRR596013_r2.fq.gz

如果直接使用这个文件名作为输入文件，celescope则会报错无法识别R1、R2文件

可识别的文件名后缀

因此需要对原fq文件重命名

mkdir st1_data

ls -d /home/jmzeng/downloads/CRR*/*f1.fq.gz |while read id ;do (crr=`echo ${id}|cut -d "/" -f 5`;ln -s ${id} ~/celescope_test/st1_data/${crr}_R1.fq.gz);done

ls -d /home/jmzeng/downloads/CRR*/*r2.fq.gz |while read id ;do (crr=`echo ${id}|cut -d "/" -f 5`;ln -s ${id} ~/celescope_test/st1_data/${crr}_R2.fq.gz);done

重命名后

为每个样本构建 celescope rna 分析的 shell 脚本

862 7月  23 00:34 ./shell/Blca_cKO1.sh
862 7月  23 00:34 ./shell/Blca_cKO2.sh
862 7月  23 00:34 ./shell/Blca_wKO1.sh
862 7月  23 00:34 ./shell/Blca_wKO2.sh
855 7月  23 00:34 ./shell/Blca_WT1.sh
855 7月  23 00:34 ./shell/Blca_WT2.sh
855 7月  23 00:34 ./shell/Blca_WT3.sh
855 7月  23 00:34 ./shell/Blca_WT4.sh

我们也可以查看构建的样本脚本的具体内容，可以看到有三个步骤，其中每一步也可以有选择的分步运行。

celescope rna sample --outdir .//Blca_cKO1/00.sample --sample Blca_cKO1 --thread 8 --chemistry auto --wells 384  --fq1 /home/yjzhang/celescope_test/st1_data/CRR596012_R1.fq.gz 
celescope rna starsolo --outdir .//Blca_cKO1/01.starsolo --sample Blca_cKO1 --thread 8 --chemistry auto --adapter_3p AAAAAAAAAAAA --genomeDir /home/yjzhang/reference/mouse/ensembl/refdata_celescope --outFilterMatchNmin 50 --soloCellFilter "EmptyDrops_CR 3000 0.99 10 45000 90000 500 0.01 20000 0.001 10000" --starMem 32 --soloFeatures "Gene GeneFull_Ex50pAS"  --fq1 /home/yjzhang/celescope_test/st1_data/CRR596012_R1.fq.gz --fq2 /home/yjzhang/celescope_test/st1_data/CRR596012_R2.fq.gz 
celescope rna analysis --outdir .//Blca_cKO1/02.analysis --sample Blca_cKO1 --thread 8 --genomeDir /home/yjzhang/reference/mouse/ensembl/refdata_celescope  --matrix_file .//Blca_cKO1/outs/filtered

提交 celescope rna 分析任务

根据自己的服务器和资源使用情况可以一次提交运行，也可分批并行提交任务

## 全部提交
nohup bash ./shell/*.sh 

##分批提交，每次运行4个
screen -r cele
conda activate celescope
cat crr_sample.txt |cut -f 2|xargs -iF -P 4 sh -c 'bash ./shell/F.sh 1>log.F 2>&1 '

结果文件

每个样本的输出结果如下。我们通常需要的关注的也就是html报告和定量矩阵文件。

结果文件目录结构

查看一下结果报告

report.html

整理我们需要的结果

pro_folder=/home/jmzeng/testceles
mkdir -p $pro_folder
ls -lh  $pro_folder
mkdir -p $pro_folder/html 
ls  */*.html |while read id;do (cp $id $pro_folder/html/${id%%/*}.html );done 
ls */outs/
mkdir -p $pro_folder/matrix 
ls -d */outs/filtered |while read id;do (cp -r  $id $pro_folder/matrix/${id%%/*} );done

整理后的结果

拿到以上文件，后面我们就可以正常走单细胞下游seurat分析流程。见：胃癌单细胞数据集GSE163558复现(二)：Seurat V5标准流程

更多用法详见：https://github.com/singleron-RD/CeleScope/blob/master/doc/user_guide.md 看其介绍，目前在准备更新nextflow流程，可以期待一下：https://github.com/singleron-RD/nf-celescope

参考：