前面的教程里面 能从源头解决数据分析的瑕疵吗 ,我们重现了普通单细胞转录组数据分析的从fastq文件开始的走cellranger的定量流程。使用的是kingfisher进行了一行代码下载数据,但是针对同一个文章的空间单细胞转录组数据就总是失败,既然没办法直接下载fastq文件,就只能说是曲线救国啦,先现在ncbi数据库里面的sra文件然后转为fastq文件即可。
下载sra文件
首先呢,需要熟读文献,在推文:数据分析有错误并不可怕,造假才不可饶恕 提到了这个新鲜出炉( 2023年12月5日)的cell期刊的文章。根据文献里面的空间单细胞转录组样品情况,制作一个文本文件;srr_list.txt ,内容就是下面的3个id即可:
SRR19996402
SRR19996403
SRR19996404
然后使用ncbi提供的SRA Toolkit 工具集里面的命令prefetch
根据上面的3个id去下载,如下所示的很简单的一句话代码即可 : nohup prefetch -O ./ --option-file srr_list.txt --max-size u 1>down.log 2>&1 &
其中prefetch
是 NCBI(National Center for Biotechnology Information)提供的一个命令行工具( SRA Toolkit 工具集),用于下载和缓存 NCBI 数据库中的生物信息学数据。它是 SRA Toolkit 工具集的一部分,专门用于处理 NCBI Sequence Read Archive(SRA)中的高通量测序数据。
上面的Linux命令是使用nohup
命令在后台运行prefetch
程序,并将输出记录到名为down.log
的文件中。以下是每个部分的详细解释:
nohup
: 该命令用于在后台运行其他命令,同时忽略SIGHUP(hangup)信号,这样即使用户退出终端,命令仍会继续执行。prefetch
: 这是 NCBI SRA Toolkit 提供的命令行工具,用于从 NCBI Sequence Read Archive(SRA)下载测序数据。它需要一个或多个参数,这里可能是通过--option-file srr_list.txt
指定一个包含SRR号的文本文件,表示要下载的数据。-O ./
: 指定数据下载的输出目录为当前目录下(./
)。--max-size u
: 这个选项可能是指定下载的数据的大小上限为 unlimited(不受限制)。这意味着,无论数据有多大,都会被下载。1>down.log
: 将标准输出重定向到名为down.log
的文件中。这意味着prefetch
程序的输出将保存在down.log
文件中。2>&1
: 将标准错误(stderr)重定向到与标准输出相同的位置,即down.log
文件。&
: 在命令末尾使用&
符号表示在后台运行该命令。
当然了,如果下载和安装NCBI(National Center for Biotechnology Information)提供的一个命令行工具( SRA Toolkit 工具集),就是很简单的conda啦,在前面的教程里面 能从源头解决数据分析的瑕疵吗 ,给出来了代码,这里就不再赘述。
命令prefetch
根据上面的3个id去下载,得到如下所示的文件:
ls -lh */*sra|cut -d" " -f 5-
36G 1月 15 05:14 SRR19996402/SRR19996402.sra
54G 1月 15 06:36 SRR19996403/SRR19996403.sra
32G 1月 15 11:04 SRR19996404/SRR19996404.sra
需要把这样的sra文件转为fastq格式的测序数据文件才可以跑流程:
ls SRR*/*sra |while read id;do (fasterq-dump -O ./ --split-files -e 8 --include-technical ${id} );done
会得到如下所示的文件哦:
ls -lh SRR1999640*gz|cut -d" " -f 5-
15G 1月 12 18:20 SRR19996402_1.fastq.gz
17G 1月 12 18:31 SRR19996402_2.fastq.gz
21G 1月 14 12:36 SRR19996403_1.fastq.gz
26G 1月 14 12:54 SRR19996403_2.fastq.gz
14G 1月 14 18:23 SRR19996404_1.fastq.gz
15G 1月 14 18:32 SRR19996404_2.fastq.gz
可以看到,每个样品最开始是单个sra的文件,转了之后是两个fastq.gz文件,理论上是足够走spaceranger定量流程啦。前面我们介绍的单细胞转录组的cellranger的定量流程可以作为参考,代码我已经是多次分享了。参考:
- 10X单细胞转录组原始测序数据的Cell Ranger流程(仅需800元)
- 10X的单细胞转录组原始数据也可以在EBI下载
- 一个10x单细胞转录组项目从fastq到细胞亚群
- 一文打通单细胞上游:从软件部署到上游分析
- PRJNA713302这个10x单细胞fastq实战
- 一次曲折且昂贵的单细胞公共数据获取与上游处理
- 只能下载bam文件的10x单细胞转录组项目数据处理
- 不知道10x单细胞转录组样品和fastq文件的对应关系
- 10X单细胞转录组测序数据的 SRA转fastq踩坑那些事
- 10x的单细胞转录组fastq文件的R1和R2不能弄混哦
因为这次是空间单细胞转录组,所以换一个软件,不再是cellranger的定量流程,应该是spaceranger,大同小异哈!
spaceranger的准备工作
主要是需要下载软件以及配套的数据库文件即可,都是需要去10x的官网下载即可,都是解压即可使用!
虽然说Space Ranger跟cellranger软件名字不一样,但是使用方法是很类似的,感兴趣的可以试试看10x官网提供的空间单细胞转录组案例数据:
基本上也是常规的fq文件命名格式即可,因为不同文件命名有不同含义,所以一定是不能弄混哦 :
- 28bp read 1 (16bp Visium spatial barcode, 12bp UMI),
- 50bp read2 (transcript)
- 10bp i7 sample barcode
- 10bp i5 sample barcode.
我们这里演示的是上面的sra文件变成的fq文件,所以是需要合理的改名的,如下所示::
ls -lh GSM*ST*gz|cut -d" " -f 5-
22 1月 13 23:41 GSM6294362_ST1_S1_L001_R1_001.fastq.gz -> SRR19996404_1.fastq.gz
22 1月 13 23:41 GSM6294362_ST1_S1_L001_R2_001.fastq.gz -> SRR19996404_2.fastq.gz
22 1月 13 23:41 GSM6294363_ST2_S1_L001_R1_001.fastq.gz -> SRR19996403_1.fastq.gz
22 1月 13 23:41 GSM6294363_ST2_S1_L001_R2_001.fastq.gz -> SRR19996403_2.fastq.gz
22 1月 13 23:41 GSM6294364_ST3_S1_L001_R1_001.fastq.gz -> SRR19996402_1.fastq.gz
22 1月 13 23:41 GSM6294364_ST3_S1_L001_R2_001.fastq.gz -> SRR19996402_2.fastq.gz
可以看到,其实每个样品只需要是有r1和r2即可,并不一定要四个文件或者三个文件哈。但是文件名字一定要符合规则哦!我上面演示了通过软连接的方式合理的给空间单细胞转录组测序数据文件修改名字哦!
运行spaceranger进行空间单细胞的上游定量
写一个脚本文本文件( run_sapceranger.sh ),如下所示:
#! /bin/bash -xe
##
bin=/home/jmzeng/x10/pipeline/spaceranger-2.1.1/spaceranger
db=/home/jmzeng/x10/pipeline/refdata-gex-mm10-2020-A
fq_dir=${2}
/usr/bin/time -v $bin count --id=${1} \
--transcriptome=$db \
--fastqs=$fq_dir \
--sample=${1} \
--image=${3} \
--unknown-slide visium-1 \
--localcores=2 \
--localmem=80
上面的脚本文本文件( run_sapceranger.sh ),可以接受3个参数,包括样品的fastq格式的测序数据文件前缀,样品的fastq格式的测序数据文件所在的文件夹路径,以及样品的图片。也就是说,是需要每个样品有fastq格式的测序数据文件的同时还需要有图片。。。。前提是作者确实是提供了,幸好我们在文章对应的geo数据库里面看到了压缩包里面是有图片的 :
ls -lh GSE207546_RAW/ |cut -d" " -f 5-
13M 1月 18 18:37 GSM6294362_ST_treatment1_spatial.tar.gz
12M 1月 18 18:37 GSM6294363_ST_treatment2_spatial.tar.gz
13M 1月 18 18:37 GSM6294364_ST_treatment3_spatial.tar.gz
我们拿其中一个样品举例,解压后可以看到 GSE207546_RAW/ST_treatment3/aligned_fiducials.jpg 就可以后面的走流程啦!所以,单个样品跑上面的脚本文本文件( run_sapceranger.sh )的代码就蛮复杂啦;
nohup bash run_sapceranger.sh GSM6294364_ST3 \
./rawdata \
./GSE207546_RAW/ST_treatment3/aligned_fiducials.jpg \
1>log_GSM6294364_ST3.txt 2>&1 &
值得注意是上面的Linux命令非常依赖于文件夹路径的准确性,文件的存在与否。也就是说要准备好3个参数,包括样品的fastq格式的测序数据文件前缀,样品的fastq格式的测序数据文件所在的文件夹路径,以及样品的图片。
让我们伤心的是很多空间单细胞转录组文献他们公开了样品的测序数据,但是缺乏对应的 图片文件,所以很难走这个定量流程。