单细胞drop-seq数据的分析流程以及debug过程
前言
单细胞数据目前除了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模块包括：
fastqc
umi_create_whitelist
whitelist_for_solo
align
index_bam
collect_rna_metrics
make_report
软件安装：
1git clone git@github.com:aselewa/dropseqRunner.git
2cd dropseqRunner
3#假设已安装conda
4conda env create -f environment.yaml
5source 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.建库：
1python ~/soft/dropseqRunner-master/makeref.py --fasta species.fa --gtf species.gtf --outDir species
这一步可能会报STAR建库内存不足的问题，比如：
1EXITING because of FATAL PARAMETER ERROR: limitGenomeGenerateRAM=31000000000is too small for your genome
解决方法为vim makeref.py,在STAR的建库代码后面添加--limitGenomeGenerateRAM 139855002325 参数，这个数值可以根据自己的物种进行调整。
2.主程序运行:
1python ~/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"去命名则不会得到分析结果，所以需要对样本名进行修改：
1python ~/soft/dropseqRunner-master/dropRunner.py --R1 SRR1_R1.fastq.gz --R2 SRR1_R2.fastq.gz --indices ~/species --protocol drop --sample SRR1
第二个bug为STAR运行时会报错，报错代码为：
1EXITING 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.批量跑样本：
该流程提供了批量跑样本的功能，使用方法为：
1R1=$(ls *_R1.fastq.gz | paste -sd,)
2R2=$(ls *_R2.fastq.gz | paste -sd,)
3
4python ~/soft/dropseqRunner-master/dropRunner.py --R1 $R1 --R2 $R2 --indices ~/species --protocol drop --sample SRR
这里也会报错，显示：
1“Please provide a gzipped fastq file for read 1”
其原因为dropRunner.py 的第107行和第108行会对R1和R2样本是否存在进行检测，但是输入多个样本时会检测失败，导致程序报错退出。解决办法是用“#”注释掉这两行即可。
结语
以上就是drop-seq的主线分析流程以及bug解决方案，流程分析的结果在形式上类似10x的分析结果，所以可以直接用seurat的Read10X()方法导入进行下一步的分析。如果是多个样本同时输入运行，不建议太多样本，因为STAR运行需要较高的内存，如果同时并行多个STAR有一定可能导致内存爆满导致卡机。