【原】10X单细胞转录组测序数据的 SRA转fastq踩坑那些事

健明 2021-11-09

展开全文

前言

笔者在学习单细胞数据分析之前，以掌握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./single-cell-gene-expression/software/downloads/latest

注：服务器下载慢的话可以在windows环境下载，然后上传压缩文件至服务器

## 1.cellranger软件6.2.1
wget -O cellranger-6.1.2.tar.gz "https://cf./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./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./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

image-20211106194423075

二. 下载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个小时后的结果）：

image-20211107205119579

生信技能树提示要理解这三个文件：

四. 使用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

image-20211107210443193

可以看到fasterq-dump的优点是转的快，问题是不能重命名及压缩。

现在将fasterq-dump用于10x单细胞的SRR，先转一个试试：

fasterq-dump -O ./ --split-3 -e 40 ./SRR7722937
ls -lh

虽然速度极快，但是fasterq-dump没有把SRR分为-1，-2，-3三个部分，这可能会影响CellRanger的后续分析。

image-20211107210816699

换一个参数也是不行：

fasterq-dump -O ./ --split-files -e 40 ./SRR7722937
ls -lh

image-20211107211127357

查了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输出的文件大小一致。

image-20211107211851372

接下来批量处理，这里有两种办法：

## 新建文件夹
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