10X单细胞转录组测序数据的 SRA转fastq踩坑那些事
考虑到咱们生信技能树粉丝对单细胞数据挖掘的需求，我开通了一个专栏《100个单细胞转录组数据降维聚类分群图表复现》，也亲自示范了几个，不过自己带娃，读博，时间精力有限，所以把剩余的90多个任务安排了学徒，实习生，学员。真的是太棒了，群策群力！
前言
笔者在学习单细胞数据分析之前，以掌握R语言和Linux基本操作。跟着生信技能树B站的视频，依次学习了GEO数据挖掘、Linux公益课、RNA-seq上下游分析、甲基化数据以及CHIP-seq上下游分析。
https://space.bilibili.com/338686099
在学习和实践锻炼了这些技能之后，开始了单细胞学习之旅。单细胞下游分析的初步走了一遍剑桥的课程 《Analysis of single cell RNA-seq data》和Seurat官方文档。对上游分析有一个初步的认识之后，笔者开始学习生信技能树B站教程 《使用10X单细胞转录组数据探索免疫治疗》，计划跟着Jimmy实战一遍单细胞的上下游分析。
一. CellRanger及参考基因组下载
这里跳过了conda安装各类软件的操作，简单介绍一下cellranger及参考基因组的下载和配置：
2021年11月的最新版本是cellranger软件6.2.1，hg38及mm10参考基因组，下载网站：https://support.10xgenomics.com/single-cell-gene-expression/software/downloads/latest
注：服务器下载慢的话可以在windows环境下载，然后上传压缩文件至服务器
## 1.cellranger软件6.2.1
wget -O cellranger-6.1.2.tar.gz "https://cf.10xgenomics.com/releases/cell-exp/cellranger-6.1.2.tar.gz?Expires=1636331862&Policy=eyJTdGF0ZW1lbnQiOlt7IlJlc291cmNlIjoiaHR0cHM6Ly9jZi4xMHhnZW5vbWljcy5jb20vcmVsZWFzZXMvY2VsbC1leHAvY2VsbHJhbmdlci02LjEuMi50YXIuZ3oiLCJDb25kaXRpb24iOnsiRGF0ZUxlc3NUaGFuIjp7IkFXUzpFcG9jaFRpbWUiOjE2MzYzMzE4NjJ9fX1dfQ__&Signature=LVCM9NLNJQ4BgDFmIdLG7jYI~hoT4QB4nNIT9U6krn~gVRHQwLP25YDI~msH6cq6fpd8OC4MTNlLNUrqVFTlBK0~qq9Cq~MA-ofUBxYt74yECf4vgGd5uRTUKsiKy2af~6~kicBaj5ltCoIjmguKVQPof4M2E81IdkMUHsVBp6NoCfQCI68sRtIMyYTx7~fpvY1MZv6PipPUXMctB85usluXc~vMQT6f-W4q-gYKG31y8wqS~4Idh4l1oQpMU4T-I16E-EO04~aLQP9lJSMuswZsUxmOP8~wkR93ULmW16UmCH6YYAXAsVZipc6xe1k12yPZ--ioAT9x4JwZW8EIUQ__&Key-Pair-Id=APKAI7S6A5RYOXBWRPDA"
md5sum cellranger-6.1.2.tar
#310d4453acacf0eec52e76aded14024c cellranger-6.1.2.tar
tar -xzvf cellranger-6.1.2.tar
## 2.mm10
wget https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-mm10-2020-A.tar.gz
md5sum refdata-gex-mm10-2020-A.tar.gz
#886eeddde8731ffb58552d0bb81f533d refdata-gex-mm10-2020-A.tar.gz
tar -xzvf refdata-gex-mm10-2020-A.tar.gz
## 3.hg38
wget https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-GRCh38-2020-A.tar.gz
md5sum refdata-gex-GRCh38-2020-A.tar.gz
#dfd654de39bff23917471e7fcc7a00cd refdata-gex-GRCh38-2020-A.tar.gz
tar -xzvf refdata-gex-GRCh38-2020-A.tar.gz
## 4.环境变量配置：PATH的路径依据自己实际的情况而定
vim ~/.bashrc
# 添加路径
PATH=/home/data/vip10t17/software_install/cellranger-6.1.2:$PATH
# 成功
cellranger
二. 下载SRA数据
视频教程用到的是GSE117888双端数据。我用服务器下载速度有点慢，因此也是用windows下载，然后上传至服务器。但是这里还是贴一下下载的脚步：
以前用ascp可以下载NCBI的SRA数据，现在似乎不行了，NCBI不提供下载的FTP链接。这点有待确认。
## 新建文件夹
mkdir 1.sra && cd 1.sra
## 构建下载list
cat >download_file
SRR7722937
SRR7722938
SRR7722939
SRR7722940
SRR7722941
SRR7722942
## 批量下载
cat download_file |while read id;do (prefetch $id &);done # 后台下载
## 情况不对就kill
## ps -ef | grep prefetch | awk '{print $2}' | while read id;do kill $id;done #批量Kill
三. fastq-dump SRA转fastq
重点来了，SRA转fastq这个过程是我写这篇帖子的原因。
先看一下常规的fastq-dump SRA转fastq，由于fastq-dump是单线程的，转的过程非常非常慢，我从下午五点到晚上九点4个小时了，还在运行当中：
## 新建文件夹
mkdir 2.fastq_raw && cd 2.fastq_raw
## 做一个软连接文件
ln -s ../1.sra/SRR* ../2.raw_fastq/
## 创建文件转换fastq.sh脚本文件
cat > fastq.sh
ls SRR* | while read id;do ( nohup fastq-dump --gzip --split-files -A $id -O ./ $id & );done
bash fastq.sh
## 查看任务进度
ps -ef |grep fastq
htop
## 情况不对就kill
## ps -ef | grep fastq-dump | awk '{print $2}' | while read id;do kill $id;done #批量Kill
经过漫长的等待，一个SRR文件才慢慢生成三个fastq文件（下图是运行4个小时后的结果）：
生信技能树提示要理解这三个文件：
四. 使用fasterq-dump加快SRR转fastq
fastq-dump固然是好使的，最大的问题就是慢，非常慢。笔者想到之前试用过的一款软件fasterq-dump，速度非常快，详见《fasterq快速转换sra文件到fastq测序数据》。但是实践过才知道，对于10X的SRR转化，fasterq-dump有些小坑。
先看看笔者之前使用fasterq-dump转双端测序的bulk RNA-seq SRR数据的效果：
fasterq-dump -O ./ --split-3 -e 40 ./SRR8980084
ls -lh
可以看到fasterq-dump的优点是转的快，问题是不能重命名及压缩。
现在将fasterq-dump用于10x单细胞的SRR，先转一个试试：
fasterq-dump -O ./ --split-3 -e 40 ./SRR7722937
ls -lh
虽然速度极快，但是fasterq-dump没有把SRR分为-1，-2，-3三个部分，这可能会影响CellRanger的后续分析。
换一个参数也是不行：
fasterq-dump -O ./ --split-files -e 40 ./SRR7722937
ls -lh
查了github的issues讨论才知道对于10x的SRR需要加--include-technical参数，原文说到：
SRR9169172 has 3 fragments per spot. they are labeled as this: technical - biological - technical you can see this yourself if you run: 'vdb-dump SRR9169172 -R1 -C READ_TYPE' fasterq-dump ignores by default the technical reads. you can force the technical reads to be written out by 'fasterq-dump SRR9169172 --include-technical'
我重写了代码：
# 转一个看看，速度很快。10x SRR需要加--split-file和--include-technical参数
fasterq-dump -O ./ --split-files -e 40 ./SRR7722937 --include-technical
ls -lh
可以看到SRR成功split为三个fastq文件。注意_3.fastq为11G，和fasterq-dump -O ./ --split-files -e 40 ./SRR7722937输出的文件大小一致。
接下来批量处理，这里有两种办法：
## 新建文件夹
mkdir 2.fastq_fasterq && cd 2.fastq_fasterq
## 软链接
ln -s ../1.sra/SRR* ../2.fastq_fasterq/
## 方法1：多个文件批量做
cat > fastq.sh
ls SRR* | while read id;do ( nohup fasterq-dump -O ./ --split-files -e 6 ./$id --include-technical & );done
# 运行脚本
bash fastq.sh
## 方法2：for循环一个一个做
for i in `ls SRR*`
do
i=$i
echo "fasterq-dump -O ./ --split-files -e 40 ./$i --include-technical"
done >fastq.sh
# 运行脚本
bash fastq.sh
## 情况不对就kill
##ps -ef | grep fasterq-dump | awk '{print $2}' | while read id;do kill $id;done #批量Kill
## 大概耗时2分钟完成
ls -lh
多折腾一会，两分钟即可完成三五个小时的程序。
对于我这种每隔一会就要去看一下进度的强迫症患者，借用群友的一句话来说，节约时间就是节约生命。
写在文末
虽然说上面的代码都是复制粘贴即可运行，但是如果要更好地完成上面的图表，通常是需要掌握5个R包，分别是: scater,monocle,Seurat,scran,M3Drop，需要熟练掌握它们的对象，：一些单细胞转录组R包的对象 而且分析流程也大同小异：
step1: 创建对象 step2: 质量控制 step3: 表达量的标准化和归一化 step4: 去除干扰因素(多个样本整合) step5: 判断重要的基因 step6: 多种降维算法 step7: 可视化降维结果 step8: 多种聚类算法 step9: 聚类后找每个细胞亚群的标志基因 step10: 继续分类
单细胞转录组数据分析的标准降维聚类分群，并且进行生物学注释后的结果。可以参考前面的例子：人人都能学会的单细胞聚类分群注释 ，我们演示了第一层次的分群。
如果你对单细胞数据分析还没有基础认知，可以看基础10讲：
咱们现在这个专栏《100个单细胞转录组数据降维聚类分群图表复现》分享的代码是到此为止，但是一般来说单细胞文章数据分析还有很多进阶图表制作，比如inferCNV看肿瘤拷贝数变异，monocle看拟时序等等。如果你也需要，可以加入我们这个专栏《100个单细胞转录组数据降维聚类分群图表复现》创作团队，获取进阶指引哦！见：急！计划招募100个单细胞爱好者，免费学全套单细胞降维聚类分群和生物学亚群注释