前面的教程里面 能从源头解决数据分析的瑕疵吗 ,我们重现了普通单细胞转录组数据分析的从fastq文件开始的走cellranger的定量流程。接下来,继续,应粉丝要求,跑一个肝癌的单细胞转录组10x数据定量流程!
在 https://www.ebi.ac.uk/ena/browser/view/PRJNA793914 可以看到这个项目详情,而且前些天我们演示了如何下载这个项目的fastq格式的测序数据原始文件,详见:aspera的高速下载确实很快吗。但是从网络下载的单细胞转录组数据文件的样品名字别抹掉了,变成了顺序编号的id,而且呢,文件名字并不符合规则:
7.1G 2月 23 02:24 SRR17418283_1.fastq.gz
17G 2月 23 02:25 SRR17418283_2.fastq.gz
13G 2月 23 02:32 SRR17418284_1.fastq.gz
13G 2月 23 02:39 SRR17418284_2.fastq.gz
16G 2月 23 03:05 SRR17418295_1.fastq.gz
17G 2月 23 03:13 SRR17418295_2.fastq.gz
11G 2月 23 03:19 SRR17418296_1.fastq.gz
12G 2月 23 03:25 SRR17418296_2.fastq.gz
5.0G 2月 23 03:28 SRR17418297_1.fastq.gz
12G 2月 23 03:34 SRR17418297_2.fastq.gz
12G 2月 23 03:40 SRR17418298_1.fastq.gz
12G 2月 23 03:46 SRR17418298_2.fastq.gz
9.4G 2月 23 03:51 SRR17418299_1.fastq.gz
9.7G 2月 23 03:55 SRR17418299_2.fastq.gz
18G 2月 23 04:04 SRR17418300_1.fastq.gz
18G 2月 23 04:13 SRR17418300_2.fastq.gz
5.8G 2月 23 04:16 SRR17418301_1.fastq.gz
13G 2月 23 04:23 SRR17418301_2.fastq.gz
4.9G 2月 23 04:26 SRR17418302_1.fastq.gz
11G 2月 23 04:32 SRR17418302_2.fastq.gz
所以是需要看网页里面的样品信息详情,然后构建一个 id.txt 的文本文件,内容如下所示:
SRR17418283,PT1
SRR17418284,PT2
SRR17418295,PT3
SRR17418296,PT4
SRR17418297,PT5
SRR17418298,TT1
SRR17418299,TT2
SRR17418300,TT3
SRR17418301,TT4
SRR17418302,TT5
然后借助于chatGPT我简单的写一个脚本批量改名,这个shell脚本会读取上面的构建好的 id.txt 的文本文件:
while IFS=',' read -r old_name new_name || [ -n "$old_name" ]; do
if [ -z "$old_name" ] || [ -z "$new_name" ]; then
continue
fi
mv "${old_name}_1.fastq.gz" "${new_name}_S1_L001_R1_001.fastq.gz"
mv "${old_name}_2.fastq.gz" "${new_name}_S1_L001_R2_001.fastq.gz"
echo "Renamed $old_name files to $new_name"
done < id.txt
最后拿到了如下所示的文件名字:
$ ls -lh *gz|cut -d" " -f 5-
4.9G 2月 24 09:19 PT1_S1_L001_R1_001.fastq.gz
11G 2月 24 09:32 PT1_S1_L001_R2_001.fastq.gz
13G 2月 21 03:54 PT2_S1_L001_R1_001.fastq.gz
13G 2月 21 04:22 PT2_S1_L001_R2_001.fastq.gz
16G 2月 21 05:21 PT3_S1_L001_R1_001.fastq.gz
17G 2月 21 05:41 PT3_S1_L001_R2_001.fastq.gz
11G 2月 21 05:56 PT4_S1_L001_R1_001.fastq.gz
12G 2月 21 06:10 PT4_S1_L001_R2_001.fastq.gz
5.0G 2月 21 00:08 PT5_S1_L001_R1_001.fastq.gz
12G 2月 21 00:45 PT5_S1_L001_R2_001.fastq.gz
12G 2月 21 06:25 TT1_S1_L001_R1_001.fastq.gz
12G 2月 21 06:40 TT1_S1_L001_R2_001.fastq.gz
9.4G 2月 21 06:53 TT2_S1_L001_R1_001.fastq.gz
9.7G 2月 21 07:02 TT2_S1_L001_R2_001.fastq.gz
18G 2月 21 01:45 TT3_S1_L001_R1_001.fastq.gz
18G 2月 21 02:18 TT3_S1_L001_R2_001.fastq.gz
5.8G 2月 21 02:29 TT4_S1_L001_R1_001.fastq.gz
13G 2月 21 02:55 TT4_S1_L001_R2_001.fastq.gz
4.9G 2月 21 03:06 TT5_S1_L001_R1_001.fastq.gz
11G 2月 21 03:28 TT5_S1_L001_R2_001.fastq.gz
也就是说,每个样品是2个文件,可以看到,有一些样品它的R1文件是远小于R2文件的,但是有一些样品的1文件和R2文件居然是大小一致的。这些文件名字就是符合要求的了,正常走cellranger的定量流程即可,代码我已经是多次分享了。参考:
- 10X单细胞转录组原始测序数据的Cell Ranger流程(仅需800元)
- 10X的单细胞转录组原始数据也可以在EBI下载
- 一个10x单细胞转录组项目从fastq到细胞亚群
- 一文打通单细胞上游:从软件部署到上游分析
- PRJNA713302这个10x单细胞fastq实战
- 一次曲折且昂贵的单细胞公共数据获取与上游处理
- 只能下载bam文件的10x单细胞转录组项目数据处理
- 不知道10x单细胞转录组样品和fastq文件的对应关系
- 10X单细胞转录组测序数据的 SRA转fastq踩坑那些事
- 10x的单细胞转录组fastq文件的R1和R2不能弄混哦
差不多几个小时就可以完成全部的样品的cellranger的定量流程。基础知识非常重要,我们在单细胞天地多次分享过cellranger流程的笔记(2019年5月),大家可以自行前往学习,如下:
- 单细胞实战(一)数据下载
- 单细胞实战(二) cell ranger使用前注意事项
- 单细胞实战(三) Cell Ranger使用初探
- 单细胞实战(四) Cell Ranger流程概览
- 单细胞实战(五) 理解cellranger count的结果
接下来我们就常规的跑一下它:
bin=/home/jmzeng/x10/pipeline/cellranger-6.1.2/bin/cellranger
# db=/home/bakdata/x10/pipeline/refdata-gex-mm10-2020-A
db=/home/jmzeng/x10/pipeline/refdata-gex-GRCh38-2020-A
ls $bin; ls $db
fq_dir=/home/jmzeng/x10/jmzeng/hcc_patients_PRJNA793914/asp/raw
$bin count --id=$1 \
--localcores=4 \
--transcriptome=$db \
--fastqs=$fq_dir \
--sample=$1 \
--expect-cells=5000
有意思的是,上面的10个样品我同一时间(2024-02-24 09:33:21)开始运行cellranger的定量流程,但是大家的耗时其实是不一样的:
log.cellranger.0.txt:2024-02-24 13:12:44 Shutting down.
log.cellranger.1.txt:2024-02-24 18:19:43 Shutting down.
log.cellranger.2.txt:2024-02-24 19:43:29 Shutting down.
log.cellranger.3.txt:2024-02-24 17:07:41 Shutting down.
log.cellranger.4.txt:2024-02-24 13:20:12 Shutting down.
log.cellranger.5.txt:2024-02-24 19:56:38 Shutting down.
log.cellranger.6.txt:2024-02-24 16:24:52 Shutting down.
log.cellranger.7.txt:2024-02-24 19:09:46 Shutting down.
log.cellranger.8.txt:2024-02-24 13:40:30 Shutting down.
log.cellranger.9.txt:2024-02-24 13:04:43 Shutting down.
也就是说有4个样品是三个半小时就跑完了cellranger的定量流程,但是另外的6个样品基本上耗时10个小时。区别在于有一些样品它的R1文件是远小于R2文件的,但是有一些样品的1文件和R2文件是大小一致的,是因为这些样品的fastq文件并没有在r1文件被trim,我们可以看2022的一个文章的示意图(The cellBC and UMI assigned Nanopore reads (median 28,120 /cell) reflect a median of 2427 genes and 6047 UMIs per cell with a good )详见:https://www.nature.com/articles/s41467-020-17800-6: