GEXSCOPE®单细胞(核)转录组数据分析
单细胞转录组测序可一次性解析成千上万个单细胞的基因表达情况，已经广泛应用于细胞图谱构建、细胞发育轨迹刻画、疾病机理研究。新格元自主开发的GEXSCOPE®单细胞转录组试剂盒可用于快速构建单细胞测序文库，可完成从样本前处理至单细胞测序文库构建全部流程。GEXSCOPE®单细胞测序试剂盒基于微流控微孔芯片技术，可一次性完成数百至数万个细胞的分离，技术稳定可靠，无样本类型偏好性，可适用于超400余种样本的单细胞测序建库需求。
在使用新格元GEXSCOPE®单细胞转录组试剂盒完成文库构建，上机测序得到测序序列之后，可以使用 celescope rna 管线进行分析。
celescope rna pipeline (单细胞转录组分析) 包含八个主要指令，可以通过 celescope rna {指令} --help 查看：
conda activate celescope
celescope rna --help
usage: celescope rna [-h]
{mkref,sample,barcode,cutadapt,star,featureCounts,count,analysis}
...
Single-cell rna
positional arguments:
{mkref,sample,barcode,cutadapt,star,featureCounts,count,analysis}
optional arguments:
-h, --help show this help message and exit
# 在运行之前检查一下软件版本
celescope -v
# 1.14.1 我们测试使用的版本
- mkref：创建基因组参考目录.
- sample：提供样本信息，包括样本ID，分析类型，scope 试剂版本和 celescope 软件版本.
- barcode：基于 read1 序列信息过滤，提取并矫正 barcode，将矫正后的 barcode 和原始的 UMI 序列添加到 read2 的 ID 中.
- cutadapt：对数据进行质控.
- star：调用 STAR 将 reads2 序列比对到基因组上.
- featureCounts：FeatureCounts 根据 bam 文件将序列比对的结果计数到细胞内.
- count：进行UMI计数以及细胞数目评估 (cell calling), 最终输出表达矩阵.
- analysis：分析表达矩阵，包括线粒体统计、tSNE 分群等.
1. 用 mkref 指令创建一个参考基因组文件
mkref 用于创建基因组参考目录，每个物种构建一次参考基因组即可，以后运行直接读取，不用多次构建。
- 需要基因组序列文件与注释信息：fasta 文件和 gtf 文件
下面我们以Homo sapiens为例来演示基因组构建过程。
第一步，下载基因组序列文件与注释信息
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
第二步，使用 celescope rna mkref 构建用于单细胞比对的基因组索引文件
conda activate celescope
celescope rna mkref \
--genome_name \#设置参考基因组的名称，这个名称会最终展示到celescope生成的报告里，默认为None
--fasta \#接fasta文件
--gtf \#接gtf文件
CeleScope 下 mkref 指令的功能是创建一个包含比对索引信息的基因组参考目录:
|-- Genome
|-- Homo_sapiens_ensembl_99.refFlat
|-- Log.out
|-- SA
|-- SAindex
|-- celescope_genome.config
|-- chrLength.txt
|-- chrName.txt
|-- chrNameLength.txt
|-- chrStart.txt
|-- exonGeTrInfo.tab
|-- exonInfo.tab
|-- geneInfo.tab
|-- genomeParameters.txt
|-- sjdbInfo.txt
|-- sjdbList.fromGTF.out.tab
|-- sjdbList.out.tab
`-- transcriptInfo.tab
2. 下载测试数据与脚本
为了便于测试软件，我们在github上托管了测试数据（请注意，这些数据仅供测试，部分数据是人为生成的）。
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
当然，也可以在gitee上下载：
mkdir test_dir
cd test_dir
git clone https://gitee.com/singleron-rd/celescope_test_data.git
git clone https://gitee.com/singleron-rd/celescope_test_script.git
让我们看一看下载的测试数据
cd celescope_test_data
tree -L 2
.
├── capture_virus
│ ├── EBV_genome
│ ├── fake_match_dir
│ └── fastq
├── citeseq
│ ├── fake_match_dir
│ └── fastqs
├── dynaseq
│ ├── a1.vcf
│ ├── control.PosTag.csv
│ ├── fastqs
│ └── gene.strandedness.csv
├── flv
│ ├── fake_match_dir
│ ├── fltest_R1.fq.gz
│ └── fltest_R2.fq.gz
├── fusion
│ ├── fastqs
│ └── match_dir
├── README.md
├── rna
│ └── fastqs
├── snp
│ ├── fastqs
│ ├── gene_list.tsv
│ └── snp_match_dir
├── sweetseq
│ ├── fake_match_dir
│ ├── fastqs
│ ├── sweet_tag_barcode.fasta
│ └── sweet_tag_linker.fasta
├── tag
│ ├── fake_match_dir
│ └── fastqs
└── vdj
├── fake_match_dir
└── fastqs
下面是测试脚本：
cd celescope_test_script
tree -L 2
.
├── capture_virus
│ ├── capture_virus.mapfile
│ ├── old.sh
│ └── run_shell.sh
├── citeseq
│ ├── CLindex_TAG.fasta
│ ├── run_shell.sh
│ └── test.mapfile
├── conftest.py
├── dynaseq
│ ├── case.mapfile
│ └── run_shell.sh
├── flv_CR
│ ├── run_shell.sh
│ └── test.mapfile
├── flv_trust4
│ ├── run_shell.sh
│ └── test.mapfile
├── fusion
│ ├── fusion_test.mapfile
│ └── run_shell.sh
├── README.md
├── rna
│ ├── rna.mapfile
│ └── run_shell.sh
├── snp
│ ├── annovar
│ ├── annovar.config
│ ├── run_gene_list.sh
│ ├── run_shell.sh
│ └── snp.mapfile
├── sweetseq
│ ├── run_shell.sh
│ └── sweetseq.mapfile
├── tag
│ ├── run_shell.sh
│ ├── smk_barcode.fa
│ └── tag.mapfile
├── test_multi.py
└── vdj
├── run_shell.sh
└── tcr.mapfile
跑通了测试数据之后，拿到我们真实的数据，其实换一下路径就可以啦。
3. 用 multi_rna 构建 celescope rna 分析的 shell 脚本
用文本编辑器编辑两个文件，一个是包含测序数据路径信息的tab分割的mapfile文件：
rna pathto/celescope_test_data/rna/fastqs rna_test
其中， 第一列：Fastq文件前缀，即序列名称 第二列：Fastq文件的文件夹路径 第三列：样品名称，将作为所有输出文件的前缀。
示例：
Sample1有2个 paired-end fastq 文件位于两个不同的文件夹 (fastq_dir1和fastq_dir2)。Sample2 有一个 paired-end fastq 文件在 fastq_dir1 中。
$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
另一个是 shell 脚本文件：test1.sh
conda activate celescope # 如果环境已经激活，可以省略
multi_rna\
--mapfile ./mapfile \# 提供mapfile文件位置
--genomeDir pathto/genome \ # mkref创建的参考基因组目录(路径)
--thread 4 \# 线程数
--mod shell \# 生成的脚本类型
--outdir ./ # 可以指定文件的输出路径
- multi_rna 是生成 CeleScope rna 分析脚本的软件，整合了八个指令中的大部分参数，该指令生成的脚本用 sh 运行后，会生成一个包含 {sample}.sh 文件的文件夹 (文件夹名为shell) ，{sample} 对应 mapfile 第三列.
- 默认--chemistry 参数可以自动识别试剂盒版本，也可以根据需要调整
- auto 默认值。适用于≥scopeV2的版本并且自动检测.
scopeV1 适用于经典的 GEXSCOPE scopeV1 文库.
customized 适用于用户自定义组合. 需要同时提供 pattern, whitelist 和 linker 的参数.
chemistry 数据在https://github.com/singleron-RD/CeleScope/tree/master/celescope/data
运行刚编辑好的shell脚本：
sh rna_test.sh
multi_rna 参数生成的脚本用 sh 运行后，会生成一个包含 rna_tes.sh 文件，是流程每一步的运行过程，可以按顺序单独运行，当然直接nohup sh ./shell/rna_test.sh & 一起也可以。花括号里面换成项目样本名。如果对每一步做了什么感兴趣，可以单独运行查看，./shell/rna_test.sh 里面是：
cat ./shell/rna_test.sh
celescope rna sample --outdir .//rna_test/00.sample --sample rna_test --assay rna --thread 4 --chemistry auto --fq1 /mnt/sdd/training/resources/test_data/rna_fastq/rna_1.fq.gz
celescope rna barcode --outdir .//rna_test/01.barcode --sample rna_test --assay rna --thread 4 --chemistry auto --lowNum 2 --fq1 /mnt/sdd/training/resources/test_data/rna_fastq/rna_1.fq.gz --fq2 /mnt/sdd/training/resources/test_data/rna_fastq/rna_2.fq.gz
celescope rna cutadapt --outdir .//rna_test/02.cutadapt --sample rna_test --assay rna --thread 4 --minimum_length 20 --nextseq_trim 20 --overlap 10 --insert 150 --fq .//rna_test/01.barcode/rna_test_2.fq
celescope rna star --outdir .//rna_test/03.star --sample rna_test --assay rna --thread 4 --genomeDir /mnt/sdd/training/Yannng/mkref/ --outFilterMultimapNmax 1 --starMem 30 --fq .//rna_test/02.cutadapt/rna_test_clean_2.fq
celescope rna featureCounts --outdir .//rna_test/04.featureCounts --sample rna_test --assay rna --thread 4 --gtf_type exon --genomeDir /mnt/sdd/training/Yannng/mkref/
当我们运行完以后，就可以得到一个单细胞转录组的网页报告。
好啦，以上就是一个完整的新格元单细胞转录组分析过程，接下来就可以进行数据分析了。另外透漏一下，得到的表达矩阵是可以上传到新格元单细胞临床数据库SynEcoSys(https://www.synecosys.cn/)进行细胞注释和分析的哟~
本文仅做软件安装测试使用，更多软件更新信息参见：https://github.com/singleron-RD/CeleScope
微信公众号
售后公众号