一次曲折且昂贵的单细胞公共数据获取与上游处理

很久以前分享了：10X单细胞转录组原始测序数据的Cell Ranger流程（仅需800元）以及一个10x单细胞转录组项目从fastq到细胞亚群，但是它缺乏NCBI的SRA数据库下载方式，因为ebi的ena数据库首先是不稳定，其次是部分单细胞数据集的样品在ena上面并不是R1,R2,I1的3个fastq文件形式。所以我们补充了"赵小明"的笔记：一文打通单细胞上游：从软件部署到上游分析，而且在：PRJNA713302这个10x单细胞fastq实战演示一下全部流程。

欢迎大家跟着我们的教程，完成自己的数据分析，并且也同步记录和分享哦！

下面是粉丝投稿的单细胞实战的分享

看文献过程中看到一篇挺好的单细胞文章，思路比较新颖，测序数据量也比较大，因此想用这批数据自己做一做。

首先从文章中找到数据Accession number（PRJNA489757 ），很奇怪，作者只建立了SRA的BioProject，未创建GSE记录，看到共有46个测序数据。

没有多想，按常规流程，走生信技能树的一站式教程：一文打通单细胞上游：从软件部署到上游分析

(1) 数据下载

## 写入需要下载的文件名
cat >download_file
SRR7904860
SRR7904861
SRR7904862
SRR7904863
SRR7904864
.......##省略全部46个SRR

# 批量下载
cat download_file |while read id;do (prefetch $id &);done  # 后台下载

(2) fasterq-dump SRA转fastq

cd 2.raw_fastq
ln -s ../1.sra/SRR* ./

### 方法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

此时检查转出的fastq发现问题：fastq-dump只能转成1个fastq。

而文章采用常规的10x scRNA-seq方案，应是双端测序，理应至少转出R1和R2两个fastq。

于是审查原始数据类型：

发现作者并未上传fastq文件，而是上传的bam，可能问题就是出在这里。

于是检查是否可以SRA转bam，从SRA-tools官网找到sam-dump函数

(3) sam-dump SRA转bam

首先尝试其中一个：

sam-dump SRR7904879.sra | samtools view -bS - > SRR7904879.bam

得到bam文件，检查标签是否正确：

samtools view SRR7904879.bam | less -SN
samtools view SRR7904879.bam | head -3 | tr "\t" "\n" | cat -n
    10    TTATATTTTCTTTTAAGGTCTGATATAACTCTGCACTAAACCCATCTGGTCCTTGGCTTTTTTTTTTTTTTTTTTTTTTTTTTGGGTGGGAGACTATT
    11    F::,,:,F:,:F:,,,:,F,,,,F,,:F::,:,::,,FFF:FF,:,:,F,,,,,F:,FFFFFFF:FFF:FFFF,FFF:FFFFFFFFFFFFFFFFFFFF
    12    RG:Z:Z-10-G:MissingLibrary:1:H5NT3DSXX:4
    13    CB:Z:GTCATTTAGGCGATAC-1
    14    UB:Z:TATACGAGGT
    15    NH:i:1
    16    NM:i:8

似乎CB、UB等标签都存在，于是利用cellranger进行下一步bam转fastq：

(4) cellranger bamtofastq bam转fastq

然而出现报错：

cellranger bamtofastq --nthreads 30 --traceback SRR7904879.bam ~/2.raw_fastq/SRR7904879
    bamtofastqerror: BAM record missing tag: "CR" on read "1". You do not appear to have the original 10x BAM file.
    If you downloaded this BAM file from SRA, you likely need to download an 'unstripped' version of the BAM.

检索报错原因，在Github上找到类似提问：'Invalid BAM record: read: "1" is missing tag: "CR"'

基于该提问，怀疑是SRA转成的bam不是原始bam文件，于是再次核查原始数据：

发现SRA下载的文件大小与Original format的文件大小明显不一致，因此考虑之前下载的SRA文件均无效，重头再来。

(5) 寻找下载原始bam的途径

除了NCBI之外，部分测序数据可能同步上传ENA，因此再去ENA上碰碰运气，可惜该篇测序数据没有上传至ENA。

难道没有能够获得原始文件的方法了么？

于是再次检索，发现除SRA和ENA外下载原始文件的新方法，依然是生信技能树Jimmy老师的教程：初步尝试从AWS下载SRA原始数据

即从Amazon Cloud中转存及下载数据。

完全按照该教程注册账户和设置存储桶，开始传递数据：

可以发现原始文件共2T+，而SRA仅不到1T，因此选择传递原始TenX（10x？）数据。

第一次传递还收到邮件，说发生了错误：

毫不犹豫直接回复邮件询问，对方也不多说直接解决：

而后成功收到传递完成的邮件：

检查存储桶中的数据并下载，速度确实如教程所说在1-2M/s。

## 在终端进行
### 查看储存桶中的数据
aws s3 ls s3://folder
    PRE SRR7904860/
    PRE SRR7904861/
    PRE SRR7904862/
    PRE SRR7904863/
    PRE SRR7904864/
    PRE SRR7904865/
    PRE SRR7904866/
    PRE SRR7904867/
    PRE SRR7904868/
...........
### 下载数据到本地
aws s3 sync s3://folder/ /Users/data/

(6) original bam转fastq

首先尝试其中一个，先检查标签类型：

samtools view SRR7904864.possorted_genome_bam.bam | head -3 | tr "\t" "\n" | cat -n
    10    TTAGTATAGCCTTTCATAGTAGAATCTGATGATGTTTTTGATATCCTCATGTTCTGTTGTTATGTCTCCTTTTTCATTTCTGATTTTGTTAATTATAG
    11    FFFA<AAF777-FJAFA<AJJJ<F-AA<A7-<FAJJFJ<AAFJJFFJFFJAJFAJFFJ<F<FJJJJJAJJJJJJFAJJFJJJFJJJJJJAA7<AAA--
    12    NH:i:1
    13    HI:i:1
    14    AS:i:92
    15    nM:i:2
    16    RE:A:I
    17    CR:Z:CAAGATCCAAGACACG
    18    CY:Z:AAAFAAJJJJJFFAJ<
    19    CB:Z:CAAGATCCAAGACACG-1
    20    UR:Z:TATTCTTTAG
    21    UY:Z:-J<AF<AJJF
    22    UB:Z:TATTCTTTAG
    23    RG:Z:3:MissingLibrary:1:HLG3NCCXY:5
    24    E00511:295:HLG3NCCXY:5:1117:9039:15531

CR标签存在，满足转fastq的要求，继续转fastq。

cellranger bamtofastq --nthreads 30 --traceback SRR7904864.possorted_genome_bam.bam ~/2.raw_fastq/SRR7904864
ls

成功转成标准的I1、R1、R2三个fastq文件，可以下一步运行cellranger：

cellranger count --id=SRR7904864 \
                   --transcriptome= ~/refdata-gex-mm10-2020-A \
                   --fastqs= ~/SRR7904864 \
                   --sample=bamtofastq \
                   --localcores=20 \
                   --localmem=64 \
                   --nosecondary
...........
Pipestance completed successfully!

后续就是写脚本批量处理啦~

下载得到的bam数据结构：

文件名不能体现样本，因此按照上级文件夹名按照SRA号重命名bam文件：

cat >filename.list
SRR7904860
SRR7904861
SRR7904862
SRR7904863
SRR7904864
.....#省略全部46个文件名，从SRA的Run Slector中下载txt文件然后复制即可
#批量改名，改成上级文件夹名字
cat filename.list |while read id
do
mv /mnt/JG/Rawdata/mouse/MaXin_NSR/1.bam/${id}/*.bam ${id}_possorted_genome_bam.bam
done
bash rename.sh

改名完成后还会自动剪切到上级目录（意外之喜）：

批量bam转fastq

cat >bam.list
SRR7904860
SRR7904861
SRR7904862
SRR7904863
......#省略
#转fastq
cat > bamtofastq.sh
cat bam.list |while read id
do
cellranger bamtofastq --nthreads 30 --traceback ${id}_possorted_genome_bam.bam /mnt/JG/Rawdata/mouse/MaXin_NSR/2.raw_fastq/${id}
done
export PATH=/mnt/JG/Cellranger_install/cellranger-6.1.2:$PATH
nohup bash bamtofastq.sh &

代码很简单， 就是需要学习和理解Linux命令而已。

(7) 花费

因为AWS存在免费存储限制，这一批2T+数据的传输超过免费额度，最终花费近1500元，明细如下：

AWS已支持银联信用卡，下载完成后记得及时删除存储桶~（我该如何向老板报销这笔钱TT）

最后，多亏技能树的详细教程，本次公共数据库获取虽一波三折但总算成功，以后还需好好学习~

学徒作业

从这个笔记里面的 PRJNA489757 项目里面 挑选任意一个样品，进行上游处理拿到10x格式的表达量矩阵3个文件，然后走单细胞降维聚类分群流程并且合理的生物学命名。仅仅是一个样品即可，节省成本哈

文末友情宣传

强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：

• 数据挖掘（GEO,TCGA,单细胞）2022年6月场，快速了解一些生物信息学应用图表
• 生信入门课-2022年6月场，你的生物信息学第一课