PRJNA713302这个10x单细胞fastq实战
》 很久以前分享了：10X单细胞转录组原始测序数据的Cell Ranger流程（仅需800元）以及一个10x单细胞转录组项目从fastq到细胞亚群，但是它缺乏NCBI的SRA数据库下载方式，因为ebi的ena数据库首先是不稳定，其次是部分单细胞数据集的样品在ena上面并不是R1,R2,I1的3个fastq文件形式。所以我们补充了"赵小明"的笔记：一文打通单细胞上游：从软件部署到上游分析，现在跟着这个笔记演示一下全部流程：
针对的是文献：《Cell-specific alterations in Pitx1 regulatory landscape activation caused by the loss of a single enhancer》，它 单细胞数据集是 GSE168632 ，但是里面并没有给表达量矩阵文件下载，所以是需要去它对应的NCBI的SRA数据库PRJNA713302这个10x单细胞fastq实战。
如果是从ENA数据库下载 PRJNA713302 项目 ，我们很容易拿到地址：
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR139/017/SRR13924917/SRR13924917.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR139/018/SRR13924918/SRR13924918.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR139/019/SRR13924919/SRR13924919.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR139/020/SRR13924920/SRR13924920.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR139/021/SRR13924921/SRR13924921.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR139/022/SRR13924922/SRR13924922.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR139/023/SRR13924923/SRR13924923.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR139/024/SRR13924924/SRR13924924.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR139/025/SRR13924925/SRR13924925.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR139/026/SRR13924926/SRR13924926.fastq.gz
把前面的文字内容存储为 fq.txt 的文件，一般来说我们服务器里面自己安装好了ascp高速下载软件，然后使用脚本批量下载，脚本内容如下：
$ cat step1-aspera.sh
cat fq.txt |while read id
do
ascp -QT -l 300m -P33001 \
-i ~/miniconda3/envs/download/etc/asperaweb_id_dsa.openssh \
era-fasp@$id .
done
# nohup bash step1-aspera.sh 1>step1-aspera.log 2>&1 &
得到文件如下所示：
6.3G 18:45 SRR13924917.fastq.gz
6.3G 19:15 SRR13924918.fastq.gz
8.2G 17:01 SRR13924919.fastq.gz
11G 17:10 SRR13924920.fastq.gz
9.1G 17:20 SRR13924921.fastq.gz
9.8G 17:31 SRR13924922.fastq.gz
13G 17:42 SRR13924923.fastq.gz
6.4G 17:46 SRR13924924.fastq.gz
6.4G 17:51 SRR13924925.fastq.gz
8.5G 19:20 SRR13924926.fastq.gz
可以看到，确实是高速下载，两个小时就完成了全部的fastq文件下载，但是它是4个样品，每个样品有多次单细胞建库测序，如下所示：
也就是说， 前面的ENA数据库下载的10个fastq文件是远远不够的， 理论上是需要10个样品的R1,R2,I1的3个fastq文件形式，但是它ENA数据库仅仅是针对每个样品提供了R2这一个fastq文件，确实了另外两个文件。
所以我们需要从NCBI的SRA数据库下载，也是很简单的prefetch命令即可，最后下载的sra文件 如下所示：
$ ls -lh */*sra |cut -d" " -f 5-
5.8G 10:10 SRR13924917/SRR13924917.sra
5.8G 09:20 SRR13924918/SRR13924918.sra
13G 09:40 SRR13924919/SRR13924919.sra
15G 09:54 SRR13924920/SRR13924920.sra
13G 09:43 SRR13924921/SRR13924921.sra
15G 10:57 SRR13924922/SRR13924922.sra
19G 10:03 SRR13924923/SRR13924923.sra
5.9G 09:23 SRR13924924/SRR13924924.sra
5.9G 09:24 SRR13924925/SRR13924925.sra
13G 09:45 SRR13924926/SRR13924926.sra
每个sra文件都需要解压，命令是 fasterq-dump -O ./ --split-files -e 10 --include-technical ，所以写脚本如下所示：
ls */*sra |while read id;do ( fasterq-dump -O ./ --split-files -e 10 --include-technical $id);done
每个sra文件都会解压出来3个fq文件，如下所示：
$ ls -lh *gz |cut -d" " -f 5-
985M 11:25 SRR13924917_1.fastq.gz
2.2G 11:25 SRR13924917_2.fastq.gz
6.7G 11:25 SRR13924917_3.fastq.gz
987M 11:39 SRR13924918_1.fastq.gz
2.2G 11:39 SRR13924918_2.fastq.gz
6.7G 11:39 SRR13924918_3.fastq.gz
这3个fq文件的大小就可以看得出来它们的格式， 分别是I1，R1，和R2，所以需要写脚本批量改名哦！
ls *_1.fastq.gz |while read id;do (pre=`basename $id|cut -d"_" -f 1`;echo $pre; ln -s $id ${pre}_S1_L001_I1_001.fastq.gz);done
ls *_2.fastq.gz |while read id;do (pre=`basename $id|cut -d"_" -f 1`;echo $pre; ln -s $id ${pre}_S1_L001_R1_001.fastq.gz);done
ls *_3.fastq.gz |while read id;do (pre=`basename $id|cut -d"_" -f 1`;echo $pre; ln -s $id ${pre}_S1_L001_R2_001.fastq.gz);done
接下来就很容易了，需要自己下载 cellranger 软件和小鼠对应的数据库文件，并且解压后，把路径给下面的代码 ：
bin=/x10/pipeline/cellranger-6.1.2/bin/cellranger
db=/pipeline/refdata-gex-mm10-2020-A
ls $bin; ls $db
fq_dir=/mice_GSE168632/input/
$bin count --id=$1 \
--localcores=8 \
--transcriptome=$db \
--fastqs=$fq_dir \
--sample=$1 \
--expect-cells=5000
每个样品独立云运行即可。
4个样品的10次单细胞转录组建库测序，每个都会独立成为一个表达量矩阵。然后批量读取后，进行简单的降维聚类分群，如下所示：
基本上跟原文是一模一样的：
可以看到，占大头的是Mesenchyme，然后小部分的上皮细胞，免疫细胞，内皮细胞以及平滑肌细胞。
不过，我比它多了一个独立的未命名的5群，它的 特异性基因如下所示：
蛮有意思的，感兴趣的可以自己去读一下原文，文献：《Cell-specific alterations in Pitx1 regulatory landscape activation caused by the loss of a single enhancer》。
学徒作业
完成这个笔记提到的PRJNA713302这个10x单细胞fastq实战，就是cellranger流程 处理fastq数据拿到表达量矩阵，然后进行基本的降维聚类分群，并且注释，搞清楚这个多了一个独立的未命名的5群是什么！
最基础的往往是降维聚类分群，参考前面的例子：人人都能学会的单细胞聚类分群注释 ，这样的基础认知，也可以看基础10讲：
写在文末
我在《生信技能树》，《生信菜鸟团》，《单细胞天地》的大量推文教程里面共享的代码都是复制粘贴即可使用的， 有任何疑问欢迎留言讨论，也可以发邮件给我，详细描述你遇到的困难的前因后果给我，我的邮箱地址是 jmzeng1314@163.com
如果你确实觉得我的教程对你的科研课题有帮助，让你茅塞顿开，或者说你的课题大量使用我的技能，烦请日后在发表自己的成果的时候，加上一个简短的致谢，如下所示：
We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.
十年后我环游世界各地的高校以及科研院所（当然包括中国大陆）的时候，如果有这样的情谊，我会优先见你。