Ji Wang

Ji Wang

No pain, no gain.

生信公共数据库数据下载流程

2 minute read

Published:

了解公共数据库

基础知识

首先了解一下SRA数据库的架构：

SRP(项目 Project)—>SRS(样本 Sample)—>SRX(数据产生 Experiment)—>SRR(数据本身)

国际上的三大生物数据库：SRA, ENA or DDBJ，分别在美国、欧洲、日本，它们之间的数据是同步的，所以可以在任意一个数据库中下载数据，而EBI数据库能够直接下载fq文件，而NCBI则需要先下载sra文件，然后再转换成fq文件（左边方框），所以推荐使用ENA数据库下载数据。或者可以直接使用sra-explorer，它是一个为了让SRA更易检索、更易下载的网页端应用。

  • 生物项目（BioProjects）是最顶层的，根据不同的数据库，它的前缀是PRJ 或者 SRP/ERP/DRP；
  • BioProjects中包含一个或多个的生物样本（BioSamples），它的前缀是SAMN 或者SRS/ERS/DRS；
  • 一个BioSample虽然只是一个样本，但它可以使用多种实验处理，也就是Experiments，前缀是SRX/ERX/DRX；
  • 每个实验都会有一个数据产出Run，它的前缀是SRR/ERR/DRR。因此，一个SRS或许会包含多个实验产生的多个数据，也就可能对应多个SRR号

NCBI找到SRR List

首先在NCBI搜BioSample的ID(DLC15103121701)， 可以进入对应的BioProject(PRJNA511588)， 然后在BioProject的页面中找到SRA Experiments的Link，在这里可以

  • 点send to Run Selector，然后就会跳转到Run Selector的页面，可以在这里点击Accession List获取所有的SRR号
  • 点send to File RunInfo，会获取一个csv文件，文件中会有所有sra的下载链接，可使用wget进行下载

如果我们任意点进一个SRX号，进去后：

  • PRJNA511588：点了就又回到了进入对应的BioProject界面
  • SRP174371：这个界面是一个项目的总览，包括了标识符、研究类型、简要介绍和对应的文章链接
  • All experiments：回到了SRA Experiments的界面
  • All runs：跳转到Run Selector，点击Accession List获取所有的SRR号

EBI找到SRR List

首先在网站ENA数据库中搜索SRR号，可以找到对应的BioProject(PRJNA511588)，然后在BioProject的页面中就可以下载全部的fq文件

获取SRR download URL后使用wget下载

链接获取上述已说

获取SRR List后使用prefetch+ascp进行数据下载

来吧，加速你的下载

conda install sra-tools=2.9.6 -y

SRA数据库如何直接下载fq文件： 使用SRA-toolkit中的工具fastq-dump直接下载SRR文件，并转换为FASTQ格式，–split-3参数表示如果是双端测序就自动拆分，如果是单端不受影响。–gzip转换fastq为压缩文件，节省空间

nohup fastq-dump -v --split-3 --gzip SRR5908360 &

使用kingfisher下载数据

转录组数据下载+文件名整理

ftp下载

下载地址的规律：ftp://ftp-trace.ncbi.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR836/SRR8366764/SRR8366764.sra 后面地址分别为：

  • SRR
  • SRR+前三个数
  • SRR号
  • SSR号.sra

GSA数据库下载

GSA数据库的申请及数据下载

不止是NCBI的SRA可以下载测序数据

使用Edge turbo下载CNCB数据

对于wget的下载方式，暂时只发现在CRR一级的能够直接下载：https://ngdc.cncb.ac.cn/gsa/browse/CRA000005/CRR007231

wget ftp://download.big.ac.cn/gsa/CRA000005/SAMC000237/CRX000077/CRR000071/CRD000132.gz

通过观察可以发现，可以使用md5文件第二列的数据进行替换实现批量下载

一些小问题

将SRS号转换成SRR号 单细胞测序的文件必须下载sra

SRA数据转换成fastq格式

​对于sra数据转换，fasterq-dump可不止快了亿点点

  • sra转换为fastq，使用的转换工具，推荐顺序为：fasterq-dump ＞ pfastq-dump ＞ fastq-dump
  • 需先下载并安装好sratoolkit，fasterq-dump和 fastq-dump都包含在sratoolkit里，而pfastq-dump需要单独下载
  • fasterq-dump最快最棒，转化速度是其他两个比不了的。但它的参数有些特殊，需要注意。

下载好数据后的处理：

cat Accession List | while read id ; do mv ./${id}/* ./ ; done # 从文件夹中将数据取出
cat Accession List | while read id; do rm -r ${id}; done # 去除文件夹
# 需要安装pigz
# conda install pigz -y
cat Accession List | while read id;do echo "fasterq-dump -e 8 --split-files -O ./ --outfile ${id}.fastq ${id}.sra";echo "pigz -p 8 -f ./${id}_1.fastq";echo "pigz -p 8 -f ./${id}_2.fastq";done > sra2fq.sh
nohup bash sra2fq.sh &

下载数据

在ENA上下载数据，以SRR为单位下载，下载的数据为fastq.gz格式，下载数据的地址等信息已经提供

  • 准备工作 ```bash

    安装prefetch

    conda install sra-tools=2.9.6 -y prefetch -h

安装ascp

conda install -c hcc aspera-cli=3.7.7 # 密钥地址要对应的修改为：-i ~/miniconda3/envs/wgbs/etc/asperaweb_id_dsa.openssh


- 下载数据
```bash
# 下载数据
## 未使用conda环境
ascp -QT -l 300m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:vol1/fastq/SRR645/002/SRR6457892/SRR6457892_1.fastq.gz ./
## 使用conda环境：设置-m可以抢点网速
ascp -QT -l 1000m -m 300m -k 1 -P33001 -i ~/miniconda3/envs/wgbs/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR558/006/SRR5582006/SRR5582006_1.fastq.gz ./

-i string输入私钥，安装 aspera 后有在目录 ~/.aspera/connect/etc/ 下有几个私钥，使用 linux 服务器的时候一般使用 asperaweb_id_dsa.openssh 文件作为私钥 -l string设置最大传输速度，比如设置为 200M 则表示最大传输速度为 200m/s。若不设置该参数，则一般可达到10m/s的速度，而设置了，传输速度可以更高。 -m 设置最小传输速度 -T 不进行加密。若不添加此参数，可能会下载不了。 –host=stringftp的host名，NCBI的为http://ftp-private.ncbi.nlm.nih.gov；EBI的为http://fasp.sra.ebi.ac.uk –user=string用户名，NCBI的为anonftp，EBI的为era-fasp –mode=string选择模式，上传为 send，下载为 recv

  • 下载多个文件 fastq_aspera地址
    fasp.sra.ebi.ac.uk:/vol1/fastq/SRR558/005/SRR5582015/SRR5582015_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR558/005/SRR5582015/SRR5582015_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR645/008/SRR6457888/SRR6457888_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR645/008/SRR6457888/SRR6457888_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR645/001/SRR6457891/SRR6457891_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR645/001/SRR6457891/SRR6457891_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR645/002/SRR6457892/SRR6457892_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR645/002/SRR6457892/SRR6457892_2.fastq.gz

aspera.sh脚本命令

cat fq.txt | while read id
do
ascp -QT -l 300m -P33001  \
-i ~/miniconda3/envs/wgbs/etc/asperaweb_id_dsa.openssh   \
era-fasp@$id  .
done

运行脚本：

nohup bash aspera.sh > step1-aspera.log 2>&1 &

数据处理

md5校验

从ENA上下载的数据，需要进行md5校验，以保证数据的完整性

# 查看下载的tsv文件
less filereport_read_run_PRJNA610526_tsv.txt | column -t | less -S
# 获取md5值，使用分号隔开了两个fastq文件
awk 'NR>2{print $9"\t"$4}' filereport_read_run_PRJNA610526_tsv.txt

# 校验md5值
# conda install md5sum -y
md5sum -c md5.txt

数据清洗

Trim galore可以自动检测adapter对NGS数据进行质量过滤，还可以查看质量分布图。也就是等于fastqc+trimmomatic。trim_galore适用于所有高通量测序，包括RRBS(Reduced Representation Bisulfite-Seq ), Illumina、Nextera和smallRNA测序平台的双端和单端数据 Trim_galore的工作流程： 第一步 首先去除低质量碱基，然后去除3’末端的adapter(如果没有指定具体的adapter，程序会自动检测前1million的序列） 第二步 对比前12-13bp的序列是否符合以下类型的adapter:

  • Illumina: AGATCGGAAGAGC # 如果输入参数 –Illumina,就会默认trimmed前13bp的adapter
  • Small RNA: TGGAATTCTCGG # 同上 如果输入参数 –small_rna,就会默认trimmed前12bp的adapter
  • Nextera: CTGTCTCTTATA # 同上 如果输入参数 –nextera,就会默认trimmed前12bp的adapter

在trim_galore运行中，会先生成trimmed文件，随后生成val文件，可以看官方讲解这二者的区别：Can you explain the difference between _val_1/_val_2 and _trimmed?

# 安装Trim galore
conda install -c bioconda trim-galore -y

# 双端测序数据
cleandata=./cleandata
rawdata=./rawdata
cat sampleID.txt | while read id
do
echo "trim_galore --phred33 -j 8 -q 20 --length 36 --stringency 3 --fastqc --paired --max_n 3 -o ${cleandata} ${rawdata}/${id}_1.fastq.gz ${rawdata}/${id}_2.fastq.gz" 
done > trim_galore.sh

nohup sh trim_galore.sh >trim_galore.log &

参数详解：

  • –phred33：指定测序数据的Phred质量值，33或者64，默认为33。如何判断呢，可以见下面两个文章，但是最近的数据基本上都是33了：

    fastq格式文件及phred33的判断 Fastq 格式说明 & (Phred33 or Phred64)

  • -j 8：指定线程数，8为最大。它的使用基于python3版本，如果报错，可以尝试更新python环境
  • -q 20：指定质量值，切除质量得分低于20的序列
  • –length：设定输出reads长度阈值小于设定值会被抛弃，默认值20bp；即小于20bp的被去除。注意，在pe150下，不要设置太高，可以50或36
  • –stringency 3：可以忍受的前后adapter重叠的碱基数，默认为1（非常苛刻）。可以适度放宽，因为后一个adapter几乎不可能被测序仪读到
  • –fastqc：当分析结束后，使用默认选项对结果文件进行fastqc分析
  • –paired：对于双端测序结果，一对reads中，如果有一个被剔除，那么另一个会被同样抛弃，而不管是否达到标准
  • –max_n 3：read中允许包含的N的总数（整数），否则将完全删除该read。 在配对末端设置中，任一read超过此限制将导致整个配对从修剪后的输出文件中删除。
  • -o：指定输出目录

数据质控

# 安装MultiQc
conda install -c bioconda multiqc -y
# 整合FastQC结果
multiqc *.zip

第四篇博客？

less than 1 minute read

Published:

This is a sample blog post. Lorem ipsum I can’t remember the rest of lorem ipsum and don’t have an internet connection right now. Testing testing testing this blog post. Blog posts are cool.

Blog Post number 3

less than 1 minute read

Published:

This is a sample blog post. Lorem ipsum I can’t remember the rest of lorem ipsum and don’t have an internet connection right now. Testing testing testing this blog post. Blog posts are cool.

博客第二篇

less than 1 minute read

Published:

This is a sample blog post. Lorem ipsum I can’t remember the rest of lorem ipsum and don’t have an internet connection right now. Testing testing testing this blog post. Blog posts are cool.