查看原文
其他

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

生信技能树 生信技能树 2022-08-10


考虑到咱们生信技能树粉丝对单细胞数据挖掘的需求,我开通了一个专栏《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
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
## 大概耗时2分钟完成
ls -lh
image-20211107212059503

多折腾一会,两分钟即可完成三五个小时的程序。

对于我这种每隔一会就要去看一下进度的强迫症患者,借用群友的一句话来说,节约时间就是节约生命。

写在文末

虽然说上面的代码都是复制粘贴即可运行,但是如果要更好地完成上面的图表,通常是需要掌握5个R包,分别是: scater,monocle,Seurat,scran,M3Drop,需要熟练掌握它们的对象,:一些单细胞转录组R包的对象 而且分析流程也大同小异:

  • step1: 创建对象
  • step2: 质量控制
  • step3: 表达量的标准化和归一化
  • step4: 去除干扰因素(多个样本整合)
  • step5: 判断重要的基因
  • step6: 多种降维算法
  • step7: 可视化降维结果
  • step8: 多种聚类算法
  • step9: 聚类后找每个细胞亚群的标志基因
  • step10: 继续分类

单细胞转录组数据分析的标准降维聚类分群,并且进行生物学注释后的结果。可以参考前面的例子:人人都能学会的单细胞聚类分群注释 ,我们演示了第一层次的分群。

如果你对单细胞数据分析还没有基础认知,可以看基础10讲:

咱们现在这个专栏《100个单细胞转录组数据降维聚类分群图表复现》分享的代码是到此为止,但是一般来说单细胞文章数据分析还有很多进阶图表制作,比如inferCNV看肿瘤拷贝数变异,monocle看拟时序等等。如果你也需要,可以加入我们这个专栏《100个单细胞转录组数据降维聚类分群图表复现》创作团队,获取进阶指引哦!见:急!计划招募100个单细胞爱好者,免费学全套单细胞降维聚类分群和生物学亚群注释


您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存