单细胞drop-seq数据的分析流程以及debug过程---生信技能树jimmy

前言

单细胞数据目前除了10x的测序数据，还有相当一部分是drop-seq的测序数据。笔者在GEO上下载了一批drop-seq的数据，在网上查找了一下没有找到详细的分析流程，想到有些大神封装好的分析流程可能放在github上，果然在上面找到了好几个流程。笔者试了其中几个，有一个名为dropseqRunner的流程可以跑通，但是有些bug。笔者便在此将这个跑通的github流程的使用方法以及出现的4个bug解决方法进行说明，方便大家后续的使用。

该流程github地址为：https://github.com/aselewa/dropseqRunner

分析流程：

dropseqRunner使用Python和Snakemake封装了drop-seq的分析流程，Snakemake drop文件包含的rule模块包括：

fastqcumi_create_whitelistwhitelist_for_soloalignindex_bamcollect_rna_metricsmake_report软件安装：git clone git@github.com:aselewa/dropseqRunner.git cd dropseqRunner #假设已安装conda conda env create -f environment.yaml source activate dropRunner

安装完成后，软件安装目录里包含以下主要文件，其中后续的debug部分需要修改makeref.py 、 dropRunner.py和Snakefile_drop.smk 这三个文件的部分代码： dropRunner.py makeref.py environment.yaml README.md Scripts Snakefile_10x.smk Snakefile_drop.smk

软件使用以及debug：

1.建库：

python ~/soft/dropseqRunner-master/makeref.py --fasta species.fa --gtf species.gtf --outDir species

这一步可能会报STAR建库内存不足的问题，比如：

EXITING because of FATAL PARAMETER ERROR: limitGenomeGenerateRAM=31000000000is too small for your genome

解决方法为vim makeref.py,在STAR的建库代码后面添加--limitGenomeGenerateRAM 139855002325 参数，这个数值可以根据自己的物种进行调整。

2.主程序运行:

python ~/soft/dropseqRunner-master/dropRunner.py --R1 SRR1.R1.fastq.gz --R2 SRR1.R2.fastq.gz --indices ~/species --protocol drop --sample SRR1

这里存在两个bug：

第一个bug输入的样本名称规范有问题，github的官方作者介绍为{}.R1.fastq.gz 格式，但这个名称格式实际上是错误的，在官方作者的Snakefile_drop.smk文件里，可以查到{samples}_R1.fastq.gz的代码，也就是说Snakefile文件里能输入的是"_R1"而不是".R1"的文件，如果按照作者的".R1"去命名则不会得到分析结果，所以需要对样本名进行修改： python ~/soft/dropseqRunner-master/dropRunner.py --R1 SRR1_R1.fastq.gz --R2 SRR1_R2.fastq.gz --indices ~/species --protocol drop --sample SRR1 第二个bug为STAR运行时会报错，报错代码为： EXITING because of FATAL ERROR in input read file: the total length of barcode sequence is 32 not equal to expected 20

这里的原因是，R1的前12个碱基为cell barcode,第13-20个碱基为UMI，加起来为20个碱基，但是实际上R1的长度为32碱基，与实际不符，所以STAR报错退出。笔者发现有些样本的R1文件为20bp，则不会报此错误。解决办法为，在Snakefile_drop.smk的STAR命令后面添加参数--soloBarcodeReadLength 0 ，该参数的作用是即使两个长度不一致，也不会报错，顺利跑完程序。

3.批量跑样本：

该流程提供了批量跑样本的功能，使用方法为：

R1=$(ls *_R1.fastq.gz | paste -sd,) R2=$(ls *_R2.fastq.gz | paste -sd,) python ~/soft/dropseqRunner-master/dropRunner.py --R1 $R1 --R2 $R2 --indices ~/species --protocol drop --sample SRR

这里也会报错，显示：

“Please provide a gzipped fastq file for read 1”

其原因为dropRunner.py 的第107行和第108行会对R1和R2样本是否存在进行检测，但是输入多个样本时会检测失败，导致程序报错退出。解决办法是用“#”注释掉这两行即可。

结语

以上就是drop-seq的主线分析流程以及bug解决方案，流程分析的结果在形式上类似10x的分析结果，所以可以直接用seurat的Read10X()方法导入进行下一步的分析。如果是多个样本同时输入运行，不建议太多样本，因为STAR运行需要较高的内存，如果同时并行多个STAR有一定可能导致内存爆满导致卡机。

---来自腾讯云社区的---生信技能树jimmy