新格元的单细胞转录组软件CeleScope实战

生信技能树

发布于 2022-03-03 13:26:45

2.3K0

发布于 2022-03-03 13:26:45

文章被收录于专栏：生信技能树

前因：带领学徒做单细胞数据集图表复现任务的时候，我选择了在《单细胞天地》公众号分享的教程里面的文献，见：哪怕是同一个单细胞数据里面的不同细胞亚群质量控制参数也可以不一样 , 但是学徒对PRJNA764023数据集的fastq原始数据处理失败了，因为原本以为是基于10X的单细胞数据集，但是走cellranger流程失败了，结果看回原文才发现所用的试剂盒是新格元的，其中走cellranger流程失败的日志如下所示：

[error] Pipestance failed. Error log at:
B19T/SC_RNA_COUNTER_CS/SC_MULTI_CORE/MULTI_CHEMISTRY_DETECTOR/_GEM_WELL_CHEMISTRY_DETECTOR/DETECT_COUNT_CHEMISTRY/fork0/chnk0-u5aeb009d49/_errors

Log message:
An extremely low rate of correct barcodes was observed for all the candidate chemistry choices for the input: Sample B19T in "/home/data/bylin/PRJNA764023/02cellranger". Please check your input data.
- 0.1% for chemistry SC3Pv3
- 0.1% for chemistry SC3Pv3HT
- 0.1% for chemistry SC5P-PE
- 0.1% for chemistry SC3Pv2
- 0.0% for chemistry SC3Pv3LT

于是安排学徒去到新格元的官方网站，有对这款试剂盒及其分析软件（celescope）的介绍，在github上有软件的使用说明及下载：https://github.com/singleron-RD/CeleScope

其官方文档也很清晰：https://github.com/singleron-RD/CeleScope/blob/master/docs/quick_start.md

在新格元的GitHub里面下载CeleScope

我们这里直接使用conda安装即可；你的服务器空白用户里面安装自己的conda，每个用户独立操作，安装方法代码如下：

# 首先下载文件，20M/S的话需要几秒钟即可
wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
# 接下来使用bash命令来运行我们下载的文件，记得是一路yes下去
bash Miniconda3-latest-Linux-x86_64.sh 
#  安装成功后需要更新系统环境变量文件
source ~/.bashrc

接下来使用conda安装aspera，后面我们下载PRJNA764023数据集的fastq文件用得着，代码如下：

conda create -n download 
conda activate download 
conda install -y -c hcc aspera-cli
conda install -y -c bioconda sra-tools
which ascp 
## 一定要搞清楚你的软件被conda安装在哪
ls -lh ~/miniconda3/envs/download/etc/asperaweb_id_dsa.openssh

我们已经多次介绍过conda细节了，这里就不再赘述。

conda管理生信软件一文就够
生信技能树B站软件安装视频
https://www.bilibili.com/video/av28836717

但是新格元的CeleScope软件需要重新建立一个conda环境。

git clone https://github.com/singleron-RD/CeleScope.git 
cd CeleScope
# 其GitHub提供了一个  conda_pkgs.txt 文件，所以需要首先 git clone 它
conda install mamba
mamba create -n celescope -y --file conda_pkgs.txt
####这个txt文件中包含了运行这个软件所需要的其他软件的版本信息 
conda activate celescope
pip install celescope

也就是说，本来是应该是很简单的 pip install 安装一个Python模块，结果因为它太多的依赖，我们怕 pip install 失败，所以首先使用conda新建了一个 celescope的小环境，在这个小环境里面安装了大量的依赖，最后才小心翼翼的 pip install celescope

测试CeleScope软件

官方教程提供了一个小数据来测试这款软件是否安装好了：https://github.com/singleron-RD/celescope_test_script

mkdir test_dir
cd test_dir
git clone https://github.com/singleron-RD/celescope_test_data.git
git clone https://github.com/singleron-RD/celescope_test_script.git

pip install pytest

cd celescope_test_script
python -m pytest -s test_multi.py --assays rna
####我们要跑的是rna数据，所以--assays后面写rna

cd celescope_test_script
python -m pytest -s test_multi.py

如果测试成功了，返回success的结果。

正式使用CeleScope软件

首先需要下载参考基因组数据和基因注释文件

我们的物种是人类，这里就很简单了：

mkdir hs_ensembl_99
cd hs_ensembl_99

wget ftp://ftp.ensembl.org/pub/release-99/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
wget ftp://ftp.ensembl.org/pub/release-99/gtf/homo_sapiens/Homo_sapiens.GRCh38.99.gtf.gz

gunzip Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
gunzip Homo_sapiens.GRCh38.99.gtf.gz

conda activate celescope
celescope rna mkref \
 --genome_name Homo_sapiens_ensembl_99 \
 --fasta Homo_sapiens.GRCh38.dna.primary_assembly.fa \
 --gtf Homo_sapiens.GRCh38.99.gtf

得到的参考基因组后续用得上

接下来需要下载PRJNA764023数据集的fastq文件

制作如下所示文件名 fq.txt 是的文本文件，内容如下所示：

fasp.sra.ebi.ac.uk:/vol1/fastq/SRR159/076/SRR15927876/SRR15927876_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR159/076/SRR15927876/SRR15927876_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR159/077/SRR15927877/SRR15927877_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR159/077/SRR15927877/SRR15927877_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR159/078/SRR15927878/SRR15927878_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR159/078/SRR15927878/SRR15927878_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR159/079/SRR15927879/SRR15927879_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR159/079/SRR15927879/SRR15927879_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR159/080/SRR15927880/SRR15927880_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR159/080/SRR15927880/SRR15927880_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR159/081/SRR15927881/SRR15927881_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR159/081/SRR15927881/SRR15927881_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR159/082/SRR15927882/SRR15927882_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR159/082/SRR15927882/SRR15927882_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR159/083/SRR15927883/SRR15927883_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR159/083/SRR15927883/SRR15927883_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR159/084/SRR15927884/SRR15927884_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR159/084/SRR15927884/SRR15927884_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR159/085/SRR15927885/SRR15927885_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR159/085/SRR15927885/SRR15927885_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR159/086/SRR15927886/SRR15927886_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR159/086/SRR15927886/SRR15927886_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR159/087/SRR15927887/SRR15927887_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR159/087/SRR15927887/SRR15927887_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR159/088/SRR15927888/SRR15927888_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR159/088/SRR15927888/SRR15927888_2.fastq.gz

一个简单的脚本就可以批量下载上面文本文件里面的全部的fastq文件啦：

cat fq.txt |while read id
do
        ascp -QT -l 300m -P33001  \
                -i /home/data/bylin/miniconda3/envs/download/etc/asperaweb_id_dsa.openssh   \
                era-fasp@$id  .
        done

得到如下所示的文件：

$ ls -lh SRR159278*gz|cut -d" " -f 5-
3.5G 2月   6 23:23 SRR15927876_1.fastq.gz
3.5G 2月   6 23:14 SRR15927876_2.fastq.gz
2.6G 2月   6 23:25 SRR15927877_1.fastq.gz
2.7G 2月   6 23:17 SRR15927877_2.fastq.gz
2.6G 2月   6 23:14 SRR15927878_1.fastq.gz
2.6G 2月   6 23:16 SRR15927878_2.fastq.gz
2.8G 2月   6 23:19 SRR15927879_1.fastq.gz
2.8G 2月   6 23:15 SRR15927879_2.fastq.gz
2.4G 2月   6 23:24 SRR15927880_1.fastq.gz
2.4G 2月   6 23:20 SRR15927880_2.fastq.gz
2.8G 2月   6 23:17 SRR15927881_1.fastq.gz
2.9G 2月   6 23:22 SRR15927881_2.fastq.gz
2.3G 2月   6 23:21 SRR15927882_1.fastq.gz
2.2G 2月   6 23:22 SRR15927882_2.fastq.gz
2.4G 2月   6 23:20 SRR15927883_1.fastq.gz
2.3G 2月   6 23:15 SRR15927883_2.fastq.gz
2.7G 2月   6 23:19 SRR15927884_1.fastq.gz
2.6G 2月   6 23:24 SRR15927884_2.fastq.gz
2.6G 2月   6 23:25 SRR15927885_1.fastq.gz
2.4G 2月   6 23:23 SRR15927885_2.fastq.gz
3.4G 2月   6 23:18 SRR15927886_1.fastq.gz
3.3G 2月   6 23:22 SRR15927886_2.fastq.gz
2.9G 2月   6 23:21 SRR15927887_1.fastq.gz
2.9G 2月   6 23:16 SRR15927887_2.fastq.gz
2.2G 2月   6 23:18 SRR15927888_1.fastq.gz
2.2G 2月   6 23:26 SRR15927888_2.fastq.gz

可以看到，绝大部分样品的数据量都很小，并不是10x那样的动辄几十个G的文件。

制作mapfile

按照官网要求，我们要先做一个mapfile，告诉celescope那些文件属于哪些sample, 必须要有3列：

第一列：Fastq 文件前缀第二列：Fastq 文件目录路径第三列：样本名称，是所有输出文件的前缀

下面是官网示例：

$cat ./my.mapfile
fastq_prefix1 fastq_dir1 sample1
fastq_prefix2 fastq_dir2 sample1
fastq_prefix3 fastq_dir1 sample2

$ls fastq_dir1
fastq_prefix1_1.fq.gz fastq_prefix1_2.fq.gz
fastq_prefix3_1.fq.gz fastq_prefix3_2.fq.gz

$ls fastq_dir2
fastq_prefix2_1.fq.gz fastq_prefix2_2.fq.gz

下面是我自己做的mapfile文件：

SRR15927876     /home/data/bylin/PRJNA764023/01download B19T
SRR15927877     /home/data/bylin/PRJNA764023/01download B11T
SRR15927878     /home/data/bylin/PRJNA764023/01download B11T
SRR15927879     /home/data/bylin/PRJNA764023/01download B9T
SRR15927880     /home/data/bylin/PRJNA764023/01download B9T
SRR15927881     /home/data/bylin/PRJNA764023/01download B9T
SRR15927882     /home/data/bylin/PRJNA764023/01download B4T
SRR15927883     /home/data/bylin/PRJNA764023/01download B4T
SRR15927884     /home/data/bylin/PRJNA764023/01download B13T
SRR15927885     /home/data/bylin/PRJNA764023/01download B13T
SRR15927886     /home/data/bylin/PRJNA764023/01download B19T
SRR15927887     /home/data/bylin/PRJNA764023/01download B2T
SRR15927888     /home/data/bylin/PRJNA764023/01download B2T

可以看到，其实是6个样品，但是每个样品会有多个测序数据。

然后运行下面的代码，会产生每个样本的.sh文件

conda activate celescope
multi_rna \
 --mapfile ./rna.mapfile\###mapfile文件
 --genomeDir /home/data/bylin/hs_ensembl_99\####参考基因组位置
 --thread 8\###运行线程数
 --mod shell###输出的文件格式

这个multi_rna命令并不会耗费计算资源，因为它就是帮你输出每个样品一个shell脚本，后面仍然是需要运行shell脚本的。

打开其中一个shell脚本文件看一下里面有什么

celescope rna sample --outdir .//B11T/00.sample --sample B11T
celescope rna barcode --outdir .//B11T/01.barcode --sample B11T
celescope rna cutadapt --outdir .//B11T/02.cutadapt --sample B11T
celescope rna star --outdir .//B11T/03.star --sample B11T
celescope rna featureCounts --outdir .//B11T/04.featureCounts --sample B11T
celescope rna count --outdir .//B11T/05.count --sample B11T
celescope rna analysis --outdir .//B11T/06.analysis --sample B11T

蛮简单的，是一个普通的 star加上featureCounts的定量流程。可以看到，就是celescope软件的multi_rna目录自动根据每个样品的fastq文件去产生的运行比对流程的命令，然后运行这些脚本。

批量提交每个样品的shell脚本

因为这个时候，其实就是6个病人的数据，所以每个病人的8个线程，还是在我的服务器承受范围内。

nohup sh ./shell/*.sh &

官网还特别强调了：必须在工作目录下运行，不应该在shell目录下运行它们。（应该是考虑到很多shell脚本不熟悉的小伙伴)

要注意的一个问题，那就是不同时间发表的文章，运用的试剂可能是不一样的，官方教程里面的脚本默认是2代以上的（scopeV2），但是这篇文章用的是1代的试剂，所以在前面设置的时候，要指定试剂，否则会报错：

multi_rna\
    --mapfile ./rna.mapfile\
    --genomeDir /SGRNJ/Public/Database/genome/homo_mus\
    --thread 8\
    --mod shell
    --chemistry scopeV1###在这里指定试剂

其他的具体参数详见官方教程：https://github.com/singleron-RD/CeleScope/blob/master/docs/rna/multi_rna.md

比对结果

比对成功后，会产生对应sample的结果，放在对应的文件夹里面

$ ls B11T/
00.sample   02.cutadapt  04.featureCounts  06.analysis
01.barcode  03.star      05.count          B11T_report.html

同样的，他会对每个样品产生相应的矩阵和流程报告，跟10x类似，供下游seurat分析的数据在每个样品的 05.count/B11T_matrix_10X 文件夹里面；

$ ls B11T/05.count/*
B11T/05.count/B11T_count_detail.txt  B11T/05.count/B11T_downsample.txt
B11T/05.count/B11T_counts.txt        B11T/05.count/stat.txt

B11T/05.count/B11T_all_matrix:
barcodes.tsv  genes.tsv  matrix.mtx

B11T/05.count/B11T_matrix_10X:
barcodes.tsv  genes.tsv  matrix.mtx

可以看到，这个时候的 genes.tsv 还是比较古老的文件名字。

总体来说，报表也是很容易理解：

报表也是很容易理解

对这个数据集的6个病人，进行降维聚类分群和生物学命名也很合理：

降维聚类分群和生物学命名

如果你确实觉得我的教程对你的科研课题有帮助，让你茅塞顿开，或者说你的课题大量使用我的技能，烦请日后在发表自己的成果的时候，加上一个简短的致谢，如下所示：

We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.

十年后我环游世界各地的高校以及科研院所（当然包括中国大陆）的时候，如果有这样的情谊，我会优先见你。

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2022-02-15，如有侵权请联系 cloudcommunity@tencent.com 删除

shell