软件简介
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数