Archives

ENA和SRA数据预处理: SRA Toolkit 和 seqtk_sra

前面一篇文章讨论了如何使用 aspera 上传/下载 神器 Aspera:

数据高速上传下载的利器:Aspera及其在生命科学中的应用

告诉我们如何高效(其实最重要的还是本地带宽, wget照样可以高效!)下载数据,遗留点下载数据后如何处理的问题,下面可以继续探讨。

SRA和ENA下载的数据,一般非标准的fastq文件格式,有时需要预处理才能进行下一步数据分析,比如:

1. 从sra文件提取fastq/sff 文件
2. 格式化fastq文件

一、数据下载

按照先前套路下载一个测试文件:http://www.ebi.ac.uk/ena/data/view/ERS980542

ENA

ascp_ena_download  ftp://ftp.sra.ebi.ac.uk/vol1/ERA536/ERA536205/fastq/RRCI25.31.fastq.gz
ascp_ena_download  ftp://ftp.sra.ebi.ac.uk/vol1/ERA536/ERA536205/fastq/RRCI25.31.fastq.gz

SRA

ascp_sra ERR1145336

二、数据格式化

2.1 安装 NCBI SRA Toolkit

下载地址: https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.8.2/

ascp_ncbi_download  ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.8.2/sratoolkit.2.8.2-centos_linux64.tar.gz
cd sratoolkit.2.8.2-centos_linux64/bin
export PATH=$PWD:$PATH

2.2 提取fastq文件
ENA下载的就是fastq文件,所以不需要额外处理,现在针对 sra 文件, 使用 fastq-dump 程序

fastq-dump的命令行参数:

$fastq-dump

Usage:
  fastq-dump [options] <path> [<path>...]
  fastq-dump [options] <accession>

Use option --help for more information

fastq-dump : 2.8.2

zhangl@localhost:/project/zhangl/sra-utils/SRA/ERR1145336$ fastq-dump  --help

Usage:
  fastq-dump [options] <path> [<path>...]
  fastq-dump [options] <accession>

INPUT
  -A|--accession <accession>       Replaces accession derived from <path> in
                                   filename(s) and deflines (only for single
                                   table dump)
  --table <table-name>             Table name within cSRA object, default is
                                   "SEQUENCE
PROCESSING

Read Splitting                     Sequence data may be used in raw form or
                                     split into individual reads
  --split-spot                     Split spots into individual reads

...

参数很多,不过大部分情况只需要如下操作:

fastq-dump  --split-3  --defline-seq '@$sn/$ri' --defline-qual '+'   ERR1145336.sra  -O  .

这样生成的fastq文件兼容性比较好了,比如(为展示方便格式进行了调整,fastq文件一个记录只有4行):

@CI25.31_9/1
GTGCCAGCAGCCGCGGTAATACGGAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTTAAGCTAGATGTGAAAGCCCC
GCGCTCAACGTGGGAGGGTCATTTAGAACTGGCAGACTAGAGTCTTGGAGAGGGGAGTGGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGG
AACATCAATGGCGAAGGCAACTCCCTGGCCAAAGACTGACGCTCATGTGCG
+
CFFGGGGGGGGGGGGGGGGHHHHGGGGGGGGGGGGGGGGGHHGGGGGHHHHHHHHGGGGHHHGGGGGGGGGGGGCGGHHHHHGGGHHHHHHHHHHHHHHG
GGGGGGGHHGHHGGHGGGGHHHHHHHHGGGGGGGGGGGGGGGGGGGGGFGFGGGGFFFFFFFFFFFFFFBFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
FFAFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFBBDDDFFFFFFFF

2.3 格式化序列

如果下载的fastq数据格式比较乱: 比如这种格式(为展示方便格式进行了调整,fastq文件一个记录只有4行):

@CI25.31_9 926R.7_206 HWI-M02034:5:000000000-A3P1T:1:1101:12893:1645 1:N:0: orig_bc=GATCTG new_bc=GATCTG \ 
      bc_diffs=0 orig_bc=GCTAC new_bc=GCTAC bc_diffs=0
GTGCCAGCAGCCGCGGTAATACGGAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTTAAGCTAGATGTGAAAGCCCCGC
GCTCAACGTGGGAGGGTCATTTAGAACTGGCAGACTAGAGTCTTGGAGAGGGGAGTGGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAACA
TCAATGGCGAAGGCAACTCCCTGGCCAAAGACTGACGCTCATGTGCG
+
CFFGGGGGGGGGGGGGGGGHHHHGGGGGGGGGGGGGGGGGHHGGGGGHHHHHHHHGGGGHHHGGGGGGGGGGGGCGGHHHHHGGGHHHHHHHHHHHHHHGGG
GGGGGHHGHHGGHGGGGHHHHHHHHGGGGGGGGGGGGGGGGGGGGGFGFGGGGFFFFFFFFFFFFFFBFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFAF
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFBBDDDFFFFFFFF

虽然格式没有问题,但是兼容性可能有问题,我们还是希望序列的名称是这样的:

1. @CI25.31_9/1
2. HWI-M02034:5:000000000-A3P1T:1:1101:12893:1645/1

使用 seqtk_sra 可以进行格式化:

seqtk_sra命令行接口:

seqtk_sra
Usage: seqtk_subseq [options] <in.fa>
Options:
  -f INT     reads name option: 0 reads name, 1: the first elements in comment fields, 2: the second ..., default: 0
  -o INT     reads orientation option: 1: forward, 2: reverse default: 1
  -v         print version number

通过下面两种模式操作就可以获得期望序列:

seqtk_sra  -f 0  -o 1 CI25.31.fastq.gz 
seqtk_sra  -f 1  -o 1 CI25.31.fastq.gz 

本文材料为 BASE (Biostack Applied bioinformatic SEies ) 课程 Linux Command Line Tools for Life Scientists 材料, 版权归 上海逻捷信息科技有限公司 所有

Last Update: 2017-09-20 3:02 PM

Comments are closed.