生信公共数据库下载处理
下载数据
基础知识
首先了解一下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号
NCBI中SRA Normalized Format和SRA Lite Format的区别
SRA Lite采用了简化的quality scores，具体它是如何简化的请参考原文
SRA Lite files are produced from SRA Normalized Format by assessing overall read quality, setting a per-read quality flag (Read_Filter), and removing base quality scores from the file. In the resulting files, all reads have a Read_Filter flag with value pass or reject. Importantly, it is still possible to produce fastq formatted files from SRA Lite format using the SRA toolkit. In this case, each read will have a constant quality score set to 30 for reads with Read_Filter value "pass" or 3 for reads with a value "reject".
Illumina fastq and sam/bam specifications support a quality bit that is set by the sequencing instrument and SRA Lite stores this as a "pass"/"reject" Read_Filter value. If this bit is set in the submitted fastq or bam file, the value is retained. If it is not, SRA will set a pass/reject value based on the quality score distribution within each read. Reads that have more than half of quality score values <20 are flagged "reject". Reads that begin or end with a run of more than 10 quality scores <20 are also flagged "reject". Reads that pass these quality checks are flagged "pass". When dumping data using the fastq-dump, fasterq-dump, or sam-dump utilities in the SRA toolkit, all reads are included by default. However, the fastq-dump tool has an option to include only passed or only rejected reads:
fastq-dump --read-filter <[pass|reject]>
In order to interact with these files and set your preference for SRA Lite files, please use SRA Toolkit version 2.11.2 or later.
简单来说，SRA Lite文件将碱基质量得分分为了pass和reject两种，pass统一给分为30，而reject统一给分为3。这里我们要注意，在后面去接头等trim的分析步骤中，有些参数设置为去除碱基质量低于36的片段，在处理SRA Lite时要小心这一点。
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文件第二列的数据进行替换实现批量下载
一些小问题
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 &
ascp下载数据
在ENA上下载数据，以SRR为单位下载，下载的数据为fastq.gz格式，下载数据的地址等信息已经提供
- 准备工作
# 安装prefetch
conda install sra-tools=2.9.6 -y
prefetch -h
# 安装ascp
conda install -c hcc aspera-cli -y
conda install -c hcc aspera-cli=3.7.7 # 密钥地址要对应的修改为：-i ~/miniconda3/envs/wgbs/etc/asperaweb_id_dsa.openssh
- 下载数据
# 下载数据
## 未使用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 文件作为私钥。使用conda下载的私钥可能在/miniconda3/pkgs/aspera-cli-3.7.7-0/etc/asperaweb_id_dsa.openssh或者/miniconda3/envs/wgbs/etc/asperaweb_id_dsa.openssh，可以使用ls -R | grep openssh查看是否存在，然后使用find . -iname 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地址：要在最后留一行空行，这是因为如果最后一行没有换行符，那么最后一行的内容就不会被读取，所以最后一行的内容就不会被下载
原因详解：【Linux】Shell脚本：while read line无法读取最后一行？？？
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
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 > aspera.log 2>&1 &
批量下载脚本/检测并下载
#!/usr/bin/env bash
# 假定要下载的文件列表为
# files=(
# "SRR13289578_1.fastq.gz"
# "SRR13289578_2.fastq.gz"
# ...
# ...
# "SRR13289606_1.fastq.gz"
# "SRR13289606_2.fastq.gz"
# )
# 定义文件名前缀和范围
prefix="SRR13289"
start=578
end=606
# 初始化文件列表
files=()
# 循环生成文件名列表
for id in $(seq $start $end); do
base="${prefix}${id}"
files+=("${base}_1.fastq.gz" "${base}_2.fastq.gz")
done
# 定义下载命令
download_cmd="ascp -QT -l 300m -P33001 -i ~/miniconda3/envs/rnaseq/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:vol1/fastq"
# ascp -QT -l 300m -P33001 -i ~/miniconda3/envs/rnaseq/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:vol1/fastq/SRR132/078/SRR13289578/SRR13289578_1.fastq.gz .
# 循环检查每个文件是否存在
echo "#!/usr/bin/env bash" > down.sh
for file in "${files[@]}"
do
if [ ! -f "$file" ]; then
# 如果文件不存在，则将下载命令写入脚本文件
# SRR为file的前六个字符
SRR=$(echo "${file:0:6}")
# srr_dir为0加上file的第十位和第十一位
srr_dir="0${file:9:2}"
# base为file的前十一个字符
base=$(echo "${file:0:11}")
${download_cmd}/$SRR/$srr_dir/$base/$file ./
fi
done
数据处理
获取样本名和md5校验
从ENA上下载的数据，需要进行md5校验，以保证数据的完整性
# 查看下载的tsv文件
# less sample.txt | column -t | less -S
# study_accession sample_accession experiment_accession run_accession tax_id scientific_name base_count fastq_md5 fastq_ftp submitted_ftp sra_ftp bam_ftp bam_bytes
# PRJNA672964 SAMN16582408 SRX9390297 SRR12926135 9823 Sus scrofa 13951863300 aaae352258023e73f1dc2634f56cda59;ce7c1558bb360b383a8833661493f1cb ftp.sra.ebi.ac.uk/vol1/fastq/SRR129/035/SRR12926135/SRR12926135_1.fastq.gz;ftp.sra.ebi.ac.uk/vol1/fastq/SRR129/035/SRR12926135/SRR12926135_2.fastq.gz ftp.sra.ebi.ac.uk/vol1/srr/SRR129/035/SRR12926135
# 删除第一行
sed -i '1d' sample.txt
# 保留fastq_md5列和fastq_ftp列，也即第8列和第9列
cat sample.txt | cut -f 8,9 > sampleID.txt
# 将;替换成/t
sed -i 's/;/\t/g' sampleID.txt
# 依次输出第一列，第三列，换行符，第二列，第四列
awk '{print $1"\t"$3"\n"$2"\t"$4}' sampleID.txt > md5.txt
# md5.txt第二列依据/分割，保留分割后的最后一列
awk '{n=split($2, a, "/"); print $1 " " a[n]}' md5.txt > md5.txt.1 && mv md5.txt.1 md5.txt
使用url获取样本名
# ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR132/072/SRR13281272/SRR13281272_1.fastq.gz
# 将文本按照/分割，保留分割后的最后一列
cat sample.txt | cut -d "/" -f 9 > sample.txt.1 && mv sample.txt.1 sample.txt
生成md5.txt以便于校验
import pandas as pd
# 读取文本文件
data = pd.read_csv('md5.txt', delimiter='\t', header=None)
# 删除第二列最后一个斜杠前面的内容
data[1] = data[1].apply(lambda x: x.rsplit('/', 1)[-1])
# 保存修改后的数据
data.to_csv('md5.txt', sep='\t', header=False, index=False)
# 校验md5值
# conda install md5sum -y
md5sum -c md5.txt
