y空间下载_repeat Nature论文:如何从SRA数据库下载数据

48a29aa4a3990e17d74dec0cc1899f86.png
Herman's Blog​zhaohuanan.cc

应用场景:

最近我需要重复一篇Nature Biotechnology的文章,研究了一下怎么从ncbi的sra数据库下载文献测序数据。

SRA数据库介绍

https://www.ncbi.nlm.nih.gov/sra/

SRA(Sequence ReadArchive)数据库是NCBI用于存储二代测序的原始数据,包括 454,Illumina,SOLiD,IonTorrent等。

我们经常会看到文献中给出数据名字为SRA然后后面接一些数字。

我们根据这个SRA的ID就可以进行下载了,然后进行数据的分析,重复文献的分析内容 。

根据SRA数据产生的特点,将SRA数据分为四类:

Studies-- 研究课题
Experiments-- 实验设计
Samples-- 样品信息
Runs-- 测序结果集

这四种分类有一个层次关系。

首先是Studies->Experiments->Samples->Runs。

这也是一个研究项目正常的逻辑关系。

了解这个层级关系,否则找sra数据就会感觉比较混乱。

一个study可能包含多个Experiment,Experiments包含了Sample、DNA source、测序平台、数据处理等信息。

SRA数据库用不同的前缀加以区分:

ERP或SRP表示Studies;
SRS 表示 Samples;
SRX 表示 Experiments;
SRR 表示 Runs。

SRA数据下载

sra数据的下载可以通过网页端下载,但是比较不方便。

  1. NCBI官方提供了SRA Toolkit软件包来进行下载。https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=software 软件的下载比较容易,都是编译好的版本,选择对应的体统,下载之后解压缩就可以使用了。
# 也可以使用bioconda直接进行安装,需要注意的是软件的名字在bioconda中是sra-tools。
conda install sra-tools

软件的用法也比较简单,根据命名我们就可以看出来,它是一个处理sra格式文件的工具包。可以用来管理和操作sra数据库的资源,里面包含了很多工具。可以处理多种测序平台的数据,这些工具大部分根据命名就知道功能。

  • fastq-dump: 最常用的,将SRA数据转换为fastq格;
  • prefetch: 下载sra数据
  • sam-dump: 将 SRA 转换为sam格式,如果原始数据是sam或bam,就需要使用这个工具;
  • sra-pileup: 生成 pileup统计结果,pileup是堆叠的意思,类似于samtools的pileup;
  • 一些不太常用的工具:
    • abi-dump: 处理abi格式数据;
    • sff-dump: 处理454测序数据;
    • 由于abi和454测序数据越来越少,相应的工具也不是特别重要了。
    • illumina-dump: 将sra转换为illumina原始的qseq文件;
    • sra-stat: 统计sra文件
    • vdb-config,vdb-decrypt,vdb-dump,vdb-encrypt,vdb-validate处理vdb格式数据。

这里我们要下载PRJNA553240的数据。

如果想知道数据的具体信息,比如数据是什么样品,采用哪种平台测序的,测序长度是多少,测序深度是多少,都可以到SRA网站上去查。

一般来说需要到sra网站点击点击点击,跳转跳转,贼麻烦,而且这个网站还反爬虫。。真是一口老血……

然后查了一下找到了最简便的方法直接下载各个Runs数据文件

  1. 首先可以使用esearch 搜索sra数据库,看这个项目下都有哪些数据。esearch默认的结果是xml格式,需要使用efetch进行解析。esearch和efetch来自NCBI的一款工具套件:entrez-edirect。conda安装即可
# 安装 
conda install entrez-direct
# 解析路径到文件
esearch -db sra -query PRJNA553240 | efetch -format runinfo >info.csv  

10ad94dbcb5521098bac5ee477f627d5.png

csv内容如下

777796f9947c48d44ab37e46a046fdd1.png

我们下载单一的SRA数据,SRR1972917,直接输入SRA ID即可。 prefetch SRR1972917

检查当前地址,如果没有就在~/ncbi下

046970ae9ef41e882abe79bbe923de73.png

我试了一下应该是prefetch版本问题,mac上的prefetch会直接下载到~/ncbi下,而服务器centos则会下载到当前文件夹

若要批量下载csv中的文件,提供一个思路:

  • 首先输入上面得到的info.csv,使用pandas解析表格并且获得Run名称和原文件名,使用prefetch下载文件并命名为原文件名
  • print ok说明下载正常
    print not found说明没找到
    print Error 说明遇到了其他问题
  • 我写了个小脚本prefetch_download_from_efetch.py
下载地址​github.com
# cd到info.csv的路径下
mkdir sra
cd sra
./prefetch_download_from_efetch.py
# 然后等待下载结束就好了,这个脚本会把info.csv中所有的文件都下载好并且重命名为原文件名
# 如果要下载的文件很多,体积很大,建议使用screen命令或者nohup命令放在后台跑

f6bd5c60aa185cf5187bafb60e7b19aa.png

默认下载的是sra格式数据,可以使用fastq-dump将sra转换为fastq了。

fastq-dump–gzip --split-3SRR1972917.sra

其实,也可以直接使用fastq-dump下载数据,下载之后直接即使fastq格式,不过还是选择prefetch比较好,因为sra数据格式比fastq格式占用空间较小,下载速度快;另一方面,sra也方便断点续传。

#直接利用fastq-dump下载数据  
 fastq-dump --split-files SRR1972917   

下载和合并hg19

  • 可以去ucsu官网下载
  • 这里图方便直接使用它给的tar包
   wget ftp://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz
   tar zxvf chromFa.tar.gz
  • 然后只留下我们需要的chr1~22XY,其他都删掉

01c7dea3265fbf66510b2d175d6b3e79.png


于是。。。。我又写了个脚本

你可以touch一个form_hg19.py然后把下面的代码粘贴进去

import os
ls = list(range(1,23))
ls.extend(["X","Y"])
para = ""
for i in ls:
    para += "chr%s.fa " % str(i)
print(para)
os.system('cat %s > hg_19.fa' % para)
print('success!')

然后运行一下

python form_hg19.py

3fa6db5e01f72a3eacd4b0cafc927841.png

然后不需要的文件全删光只留下hg_19.fa即可,less查看一下

58aaed083fc9141e1614a5612222f31e.png

c56925ca08221aa52274d9a0d43bf48f.png

ok了

后面可以供mapping使用

这样我们就学会了批量下载SRA数据库提供的测序数据了

同时也学了一下如何合并生成一个人类基因组fasta文件

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值