软件简介
Celescope可从二代测序下机的原始fastq数据开始处理，经过细胞标签的提取、质控与校正，测序数据质控，参考基因组比对，基因定量，UMI纠错与计数后确定细胞数，最终得到数据的质控报告和细胞的表达矩阵，用于后续分析，具有灵活、准确、全面的特点，是非常有力的单细胞转录组测序数据处理软件。
环境配置
conda
linux
minimum 32GB RAM(to run STAR aligner)
下载安装celescope
编写运行如下代码进行下载安装：
git clone https://github.com/zhouyiqi91/CeleScope.git
cd CeleScope
source setup.sh
如果没有报错，就说明celescope安装成功。
下载参考基因组并生成index文件
不管用什么软件，做什么分析，参考基因组都是必不可少的。
从ensembl官网下载人类基因组的参考序列文件（.fa）和基因组注释文件（.gtf）：
wgetftp://ftp.ensembl.org/pub/release-99/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
wgetftp://ftp.ensembl.org/pub/release-99/gtf/homo_sapiens/Homo_sapiens.GRCh38.99.gtf.gz
解压参考基因组文件到指定文件夹：
mkdir -preferences/Homo_sapiens/Ensembl/GRCh38
gzip -c -d Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz> references/Homo_sapiens/Ensembl/GRCh38/Homo_sapiens.GRCh38.fa
gzip -c -d Homo_sapiens.GRCh38.99.gtf.gz >
references/Homo_sapiens/Ensembl/GRCh38/Homo_sapiens.GRCh38.99.gtf
Note：运行celescope要激活conda环境。
调用STAR生成参考基因组的index文件。
conda activate celescope
gtfToGenePred -genePredExt -geneNameAsName2references/Homo_sapiens/Ensembl/GRCh38/Homo_sapiens.GRCh38.99.gtf /dev/stdout |\
awk '{print$12"\t"$1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}'> references/Homo_sapiens/Ensembl/GRCh38/Homo_sapiens.GRCh38.99.refFlat
STAR --runMode genomeGenerate\
--runThreadN 6\
--genomeDir references/Homo_sapiens/Ensembl/GRCh38 \
--genomeFastaFiles references/Homo_sapiens/Ensembl/GRCh38/Homo_sapiens.GRCh38.fa \
--sjdbGTFfile references/Homo_sapiens/Ensembl/GRCh38/Homo_sapiens.GRCh38.99.gtf \
--sjdbOverhang 100
STAR结果中会生成reference文件夹，里面包含人类基因组的index信息，如染色体、外显子等。
$ ls celescope_test/references/Homo_sapiens/Ensembl/GRCh38/
chrLength.txt exonGeTrInfo.tab genomeParameters.txt SA sjdbList.out.tab
chrNameLength.txt exonInfo.tab Homo_sapiens.GRCh38.99.gtf SAindex transcriptInfo.tab
chrName.txt geneInfo.tab Homo_sapiens.GRCh38.99.refFlat sjdbInfo.txt
chrStart.txt Genome Homo_sapiens.GRCh38.fa sjdbList.fromGTF.out.tab
小鼠及其他物种的参考基因组下载和index文件生成方法同理。
至此celescope分析的前期准备工作已经差不多完成了，下面开始正式分析。
Celescope可以用于Single cell RNA-seq，Single cell VDJ和Single cell Multiplexing。
Single cell RNA-seq（单细胞转录组分析）
激活conda环境：
conda activate celescope
编写如下脚本进行单样本分析：
celescope rna run\
--fq1/SGRNJ/DATA_PROJ/2003/20200710/LC20062911_2_L1_1.fq.gz\
--fq2/SGRNJ/DATA_PROJ/2003/20200710/LC20062911_2_L1_2.fq.gz\
--genomeDir /SGR/references/Homo_sapiens/Ensembl/GRCh38\
--sample BEPM\
--thread 4\
--chemistry auto
NOTE：运行之前要先下载好fastqc软件
输入：
--fq1 双端测序FASTQ read 1的路径
--fq2 双端测序FASTQ read 2的路径
--genomeDIR 参考基因组的路径
--sample 样本名
--thread 分析使用的线程数。在RNA-seq分析中最好不要超过8个，否则容易报错
Single cell RNA-seq还支持多样本运行模式，接口为multi_{assay}。
编写并运行如下脚本：
multi_rna\
--mapfile /SGRNJ02/RandD4/test/20200713.mapfile\
--chemistry scopeV2.1.1\
--genomeDir/SGRNJ/Public/Database/genome/homo_mus\
--thread 4\
--modshell
输入：
--mapfile：包含三列, 每列之间用tab分割；每一行是一个样本。
第一列：fastq前缀
第二列：fastq所在文件夹
第三列：{sample}（即生成文件的前缀）
第四列：可选，期望细胞数（scRNA-Seq）或者match_dir(scVDJ) 注意：当一个样本有多个fastq，且这些fastq不在同一个文件夹下时，每个fastq占一行，sample名称相同即可。
mapfile示例：
$ cat /SGRNJ02/RandD4/test/20200713.mapfile
LC20062911 /SGRNJ/DATA_PROJ/2003/20200710 S062907-3
$ ll/SGRNJ/DATA_PROJ/2003/20200710
total 26181688
-rw-r--r--. 1download ssh.bioinfo 3056870860 Jul 10 13:52 LC20062911_2_L1_1.fq.gz
-rw-r--r--. 1download ssh.bioinfo 3105319350 Jul 10 14:04 LC20062911_2_L1_2.fq.gz
运行后会在当前目录下生成一个shell文件夹，里面包含一个与sample名相同的shell脚本
$ ls -l
-rw-r-----. 1zhouxin ssh.randd 1504 Dec 9 14:20S062907-3.sh
在根目录下运行该脚本即可开始RNA-seq分析：
$ shshell/S062907-3.sh
输出：
Single cell RNA-seq输出文件包括每一步的结果文件和HTML报告。HTML报告可以查看每一步分析结果的统计，点击标题旁的？可以查看每个分析指标的具体含义。基因表达矩阵在05.count目录下。
$ ll S062907-3
total 5028
drwxr-xr-x. 2zhouxin ssh.randd 30 Dec 9 14:53 00.sample
drwxr-xr-x. 2zhouxin ssh.randd 128 Dec 9 15:31 01.barcode
drwxr-xr-x. 2zhouxin ssh.randd 89 Dec 9 15:54 02.cutadapt
drwxr-xr-x. 3zhouxin ssh.randd 4096 Dec 9 17:50 03.STAR
drwxr-xr-x. 2zhouxin ssh.randd 186 Dec 9 17:54 04.featureCounts
drwxr-xr-x. 3zhouxin ssh.randd 243 Dec 9 18:12 05.count
drwxr-xr-x. 2zhouxin ssh.randd 99 Dec 9 18:17 06.analysis
-rw-r--r--. 1zhouxin ssh.randd 5141357 Dec 9 18:17S062907-3_report.html
$ ll 05.count/
total 367216
-rw-r--r--. 1 zhouxin ssh.randd 20154 Dec 9 18:02 barcode_filter_magnitude.pdf
-rw-r--r--. 1 zhouxin ssh.randd 355985586Dec 9 17:57 S062907-3_count_detail.txt
-rw-r--r--. 1 zhouxin ssh.randd 17730119 Dec 9 18:02 S062907-3_counts.txt
-rw-r--r--. 1 zhouxin ssh.randd 229 Dec 9 18:12 S062907-3_downsample.txt
drwxr-xr-x. 2 zhouxin ssh.randd 77 Dec 9 18:03 S062907-3_matrix_10X
-rw-r--r--. 1 zhouxin ssh.randd 2279564 Dec 9 18:05 S062907-3_matrix.tsv.gz
-rw-r--r--. 1 zhouxin ssh.randd 176 Dec 9 18:12 stat.txt
NOTE：
S062907-3_matrix.tsv.gz：tab分隔的表达矩阵，行名是HGNC gene symbol, 列名是cell barcode。
S062907-3_matrix_10X：10X格式的表达矩阵，可以用Seurat的Read10X函数读入。
Single cell VDJ（单细胞免疫组库分析）
激活conda环境：
conda activatecelescope
编写并运行如下脚本：
celescope vdjrun
--fq1/SGRNJ/DATA_PROJ/2003/20200710/LC20062911_2_L1_1.fq.gz\
--fq2/SGRNJ/DATA_PROJ/2003/20200710/LC20062911_2_L1_2.fq.gz\
--sample S062907-3\
--thread 8\
--type TCR\
--chemistry auto
输入：
--fq1 FASTQ read1文件
--fq2 FASTQ read2文件
-- sample 样本名
-- thread 线程数
--type T细胞或B细胞受体
--match_dir 与RNA-seq目录匹配
输出：
$ ll S062907-3/
total 2888
drwxr-xr-x. 2zhouxin ssh.randd 30 Dec 9 15:08 00.sample
drwxr-xr-x. 2zhouxin ssh.randd 128 Dec 9 16:46 01.barcode
drwxr-xr-x. 2zhouxin ssh.randd 89 Dec 9 17:05 02.cutadapt
drwxr-xr-x. 2zhouxin ssh.randd 235 Dec 9 18:41 03.mapping_vdj
drwxr-xr-x. 2zhouxin ssh.randd 194 Dec 9 18:41 04.count_vdj
-rw-r--r--. 1zhouxin ssh.randd 2956241 Dec 9 18:41S062907-3_report.html
04.count_vdj目录包含如下文件：
S062907-3_cell_confident.tsv：VDJ cell barcode的克隆型，每一条链（TRA和TRB）分别占一行。
S062907-3_cell_confident_count.tsv：VDJ cell barcode的克隆型，每个细胞占一行。
S062907-3_clonetypes.tsv：VDJ cell barcode每种克隆型的计数和百分比。
S062907-3_match_clonetypes.tsv：VDJ cell barcode与sc-RNA-Seq cell
barcode交集的每种克隆型的计数和百分比。当提供了match_dir参数时，才会生产该文件。
Single cell Multiplexing（拆分）
激活conda环境：
conda activatecelescope
编写并运行如下代码：
celescope smkrun\
--fq1 {smk fq1.gz}\
--fq2 {smk fq2.gz}\
--sample {sample name}\
--SMK_pattern L25C45\
--SMK_barcode {SMK barcode fasta}\
--SMK_linker {SMK linker fasta}\
--match_dir {match_dir}\
--dim 2\
--combine_cluster {combine_cluster.tsv}
输入：
--SMK_pattern 必需参数。L25C45 指25 bp linker + 45 bpcell barcode
C: cell barcode
U: UMI
T: polyT
L: linker
--SMK_barcode 必需参数，标签的fa文件
--SMK_linker 必需参数，linker的fa文件
--match_dir 必需参数，运行完celescope后与scRNA-seq目录进行匹配
--dim 必需参数，规定纬度
--combine_cluster 可选参数，整合cluster文件
第一列：原始cluster数
第二列：整合后的cluster数
网友评论