工欲善其事必先利其器
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版本，可以先尝试最低要求版本，然后再根据提示调整适配的版本
代码语言：javascript复制#先尝试最低要求版本，
pip install numpy==1.21
#根据提示调整适配的版本
pip install numpy==1.22.4
测试
调整
测试，安装成功
测试运行通过
功能模块
功能模块
最小化使用
我们还是以单细胞转录组数据定量为例
首先还是构建参考基因组文件
代码语言：javascript复制##激活环境
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)
代码语言：javascript复制##使用前先激活环境
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
输出文件大小
代码语言：javascript复制 95 7月 22 16:58 mt_gene_list.txt
1.2G 7月 22 16:58 Mus_musculus.GRCm39.111.filtered.gtf
参考基因组文件
代码语言：javascript复制 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 脚本
代码语言：javascript复制multi_rna
--mapfile ./rna.mapfile
--genomeDir ~/reference/mouse/ensembl/refdata_celescope
--thread 8
--mod shell
此时需要注意一下文原始fq文件的文件命名
由于GSA的数据下载下来文件名为：
代码语言：javascript复制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文件重命名
代码语言：javascript复制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
我们也可以查看构建的样本脚本的具体内容，可以看到有三个步骤，其中每一步也可以有选择的分步运行。
代码语言：javascript复制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 分析任务
根据自己的服务器和资源使用情况可以一次提交运行，也可分批并行提交任务
代码语言：javascript复制## 全部提交
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
整理我们需要的结果
代码语言：javascript复制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
参考：
- https://cn.singleronbio.com/product/detail-26.html
- https://mp.weixin.qq.com/s/DMRVuVc635XyPmDz1lSQDA
- https://mp.weixin.qq.com/s/sU4Rl_USalAPsy2xCdTSrA
