课题需要某一篇已发表文章的基于10x Genomics平台的scATAC-seq数据,于是乎,开始了踩坑指南。
No.1 数据下载
众所周知,公共数据的下载一般主要有GEO和ENA两个大的数据库,且由于ENA支持(或者说更便于)Aspera(特点:下载速度非常快,通常能达到百M/s),于是一般情况下,利用ENA的Aspera链接进行数据下载是我的首选。
▲ 然而,对于10x scATAC-seq的数据,这个方法会出现什么坑呢?
那就是网上很多人已经说过的,这一途径下载的文件不全,导致无法进一步用cellRanger-atac进行分析,因为总是会报错,有一个啥啥啥文件没找到。
(这一部分由于前期的试错和下载数据都删了,贴一下别人的:ENA和SRA数据库存放的单细胞转录组测序数据并不一致)
于是乎,我又开始了检索之路,那么如何快速从GEO下载呢?
▲ 首先GEO下载的SRA数据进行下一步分析的主要限速步骤有两个“
(1) 一个是下载:prefetch的下载速度大概是5-6M/s,一般下载一个scATAC-seq的SRA用时至少半个小时(可能有误,因为已经忘记了)
(2) 另一个是SRA转fastq:采用以往大家熟知的fastq-dump进行SRA到fastq的转化单个文件大约需要10-20min(也是大概时间,具体我已经忘了,至少我在今天早上8:46提交的fastq-dump,到9:34拆分完全)
然后我找到了这个生信 | scATAC-Seq数据获取(超快sra转fq)
主要参考了这一部分:
No.2 SRA数据拆解为fastq文件
嗯。。速度确实很快,得到的fastq也是完整的4个。。一个缺点(在原本的文章中也提到了)就是fasterq-dump拆解的fastq不是.gz的压缩形式,这样需要较大的内存空间。
目前看起来似乎不错,于是乎,我开始了美滋滋的往下:
cellranger-atac
cellranger-atac count \ --id=$i \ --fastqs=scATAC/rawfastq \ --sample=$i \ --reference=refdata-cellranger-arc-GRCh38-2020-A-2.0.0 \ --localcores=8 \ --localmem=64
比对输出也没问题。
但是下一步我要利用cellranger-atac得到的bam文件进行进一步比对之后,出现了问题了:
继续查找原因,终于看到了一篇有用的文章:bam转fastq一直WARNING,并且两个fastq结果文件都是空的
这个文章大概的解释是:”之所以这样报错,是因为SRR3286802-24999的两条reads(/1, /2)只有一条reads比对成功了,而另一条reads没有比对成功,所以在bam文件中没有mate的记录。”
由于没有在fasterq-dump中找到如上所述的--defline-seq定义reads name的参数,于是我不得不重新再转回SRA,再次利用fastq-dump来转换。
No.3 Fatsq-dump拆解SRA生成I1,R1,R2,R3
又一个众所周知的事情是,在运用fastq-dump时,我们需要加入--split-3 这个参数,否则拆解出来只会有一个被merge到一起的fastq.gz文件。
然而,--split-3确实只能拆出SRRxxxxxx_1.fastq.gz和SRRxxxxxx_2.fastq.gz两个文件,很显然,这样得到的fastq连ENA下载的3个文件都比不上。
那么到底怎么样才能得到有I1,R1,R2,R3同时又不会被改了reads name的fastq呢?
查看fastq-dump --help的帮助文件,我发现出了以往常用的--split-3这个参数,还有一个--split-files ,这一参数的描述为:Dump each read into separate file.Files will receive suffix corresponding to read number,于是我将fastq-dump的拆解命令写为:
fastq-dump --gzip --split-files --defline-seq '@$ac-$si/$ri' -O <outfile.path> SRRxxxxxxx.sra
最终it worked, 成功得到了.gz形式的4个fastq文件。
开始测试下一步。如果再有坑再更新。