生信公共資料庫下載處理

CASTWJ發表於2024-04-15
資料庫

下載資料

基礎知識

首先了解一下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檔案第二列的資料進行替換實現批次下載

一些小問題

將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 &

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

相關文章