二、数据下载
2021-02-24 本文已影响0人
白米饭睡不醒
1.准备工作
- 工作目录管理
## 示例如下：tree -L 1
├── biosoft #软件安装目录
├── database #数据库存放目录，包括参考基因组，注释文件，公共数据库等
├── pipline #流程搭建目录
├── project #项目分析目录
└── project_backup #数据备份目录
- 详细命令
# 进入到个人目录
cd ~
## 1.建立数据库目录
mkdir database
# 在数据库下建立参考基因组数据库，注意命名习惯：数据库名字/物种名字/参考基因组版本信息
mkdir -p database/genome/Ensembl/Homo_sapiens/GRCh38_release95
## 2.建立软件安装目录
mkdir biosoft
## 3.建立项目分析目录
mkdir project
cd project
# 注意项目命名习惯：物种-样本数-疾病-分析流程
mkdir Human-16-Asthma-Trans
mkdir Human-16-Asthma-sRNA Human-16-Asthma-lncRNA Human-16-Asthma-RRBS
## 4.建立数据备份目录
mkdir project_backup
# 哪天备份就在那个日期文件夹，非常有利于项目溯源
mkdir 20200622
# 在project_backup下放一个所有备份数据路径的txt文件，方便快速查找项目
touch zhangj_project_backup_all_20200622_update.xls
2.数据下载
示例文章及其流程
第一步
在文章中找到GSE编号 → 在GEO页面的accession搜索GSE编号 → 在结果页面寻找 BioProject：和 SRA：两个编号
第二步：
用 conda 安装 aspera 软件（记得激活小环境再安装）
#激活小环境
conda activate rnaseq
#安装两种方法
conda install -c rpetit3 aspera-connect
conda install -c hcc aspera-cli
第三步：
🍀 两种下载方法
2.1
2.2
🌱 方法一：
2.3
得到Accession list文件后将文件上传到服务器 ~/project/Human-16-Asthma-Trans/data/rawdata/sra 位置
2.4
# 使用prefetch命令下载单个文件：如SRR1039510
prefetch SRR1039510
# 批量下载：建立循环，并查看
# =号两边没有空格
outputdir=/teach/project/Human-16-Asthma-Trans/data/rawdata/sra
cat sampleId.txt | while read id
do
echo "prefetch ${id} -O ${outputdir} "
done >download.sh
# 运行脚本，由于服务器资源有限，运行不报错就好了，ctrl+C强制退出
nohup sh download.sh >download.log &
# 验证数据的完整性
vdb-validate SRR1039510
* 小技巧：如何查看自己需要的信息在第几列：
less -S filereport_read_run_PRJNA* |head -n 1 | awk '{for(i=1;i<=NF;i=i+1){a[NR,i]=$i}}END{for(j=1;j<=NF;j++){str=a[1,j];for(i=2;i<=NR;i++){str=str " " a[i,j]}print str}}' |less -SN
🌿 方法二：
打开网址：https://www.ebi.ac.uk/ena/browser/home
2.5
2.6
在路径下：~/project/Human-16-Asthma-Trans/data/rawdata/sra
#从老师那调取文件（自己下载的话直接传到以上路径就好）
ln -s /teach/data/airway/sra/* .
#查看
less -S filereport_read_run_PRJNA229998_tsv.txt
#查看自己需要的信息在第几列(该命令会输出所有列名，找到需要的信息在第几列）
less -S filereport_read_run_PRJNA* |head -n 1 | awk '{for(i=1;i<=NF;i=i+1){a[NR,i]=$i}}END{for(j=1;j<=NF;j++){str=a[1,j];for(i=2;i<=NR;i++){str=str " " a[i,j]}print str}}' |less -SN
#根据需要信息在13列，修改$后的数值
cat filereport_read_run_PRJNA229998_tsv.txt |awk 'NR>1{print $13}' >sra.url
#可用以下查看
cat filereport_read_run_PRJNA229998_tsv.txt |awk 'NR>1{print $13}'|less
#查看是否有特殊字符
cat -A sra.url
#去除特殊字符
sed -i "s/\s*$//g" sra.url
2.7
下载
#查找密钥
find ~ -name asperaweb_id_dsa.openssh
#激活小环境
conda activate rnaseq
#下载sra格式数据（自己还没尝试成功，~后加的是自己的密钥）
ascp -k 1 -QT -l 300m -P33001 -i ~/miniconda3/envs/rnaseq/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:/vol1/srr/SRR103/008/SRR1039508 .
#下载gz格式
ascp -k 1 -QT -l 300m -P33001 -i ~/biosoft/miniconda3/envs/rna/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR103/000/SRR1039510/SRR1039510_1.fastq.gz .
#转后台
Ctrl+z → 输入命令 bg → jobs 查看（top)
#转前台
fg %1
#提交到后台并行下载
cat sra.url |while read id
do
nohup ascp -k 1 -QT -l 300m -P33001 -i ~/miniconda3/envs/rnaseq/etc/asperaweb_id_dsa.openssh era-fasp@${id} ./ &
done
#查看
ps fx | less
#杀掉一个任务
ps fx | grep 'ascp' | awk '{print "kill -9" $1}
##完整性检验
#(删除中间文件)
rm *.partial
rm *aspera-ckpt
#查看md5值在第几列（run_accession , sra_md5)
less -S filereport_read_run_PRJNA* |head -n 1 | awk '{for(i=1;i<=NF;i=i+1){a[NR,i]=$i}}END{for(j=1;j<=NF;j++){str=a[1,j];for(i=2;i<=NR;i++){str=str " " a[i,j]}print str}}' |less -SN
# 得到md5值，$11与$4之间为两个空格，检验格式要求(看生成的数据名字是否有sra，有的话就写，没有就不写）
awk 'NR>1{print $11" "$4".sra"}' filereport_read_run_PRJNA229998_tsv.txt >md5.txt
# md5值检验（挂后台检验）
md5sum -c md5.txt >check &
#查看检验结果
cat check
#生成md5值
md5sum sampleId.txt（文件名） > md5_test.txt
#sra转换为fq格式（可选）
# 定义文件夹（路径改成自己的）
fqdir=/teach/project/Human-16-Asthma-Trans/data/rawdata/fastq
# 单个转换(用fastq-dump -h查看split-几，改在下面）
fastq-dump --gzip --split-e -X 25000 -O ${fqdir} SRR1039510.sra
#fasterq-dump --split-files SRR11180057.sra
# 批量转换
ls ~/project/Human-16-Asthma-Trans/data/rawdata/sra/*sra | while read id
do
echo "fastq-dump --gzip --split-e -X 25000 -O ${fqdir} ${id}"
done >sra2fq.sh
# 提交后台运行命令
nohup sh sra2fq.sh >sra2fq.log &
#查看
jobs
2.8
-
注意：
2.9
