RNA-seqその２、データのダウンロードと変換

シークエンスのデータの入手方法には

公共のデータベースからシークエンスデータをダウンロードする
シークエンスを行う

などがある。

今回は公共のデータベースのデータを使う
公共データベースによって登録されているファイル形式は異なる
NCBIはSRA、DDBJ/EMBLEはfastq形式である
シークエンスの結果は一般にはfastq形式で
SRAはFastqに情報を追加して圧縮した形式である

解析にはFastq形式を使用するので、ツールを用いてSRAからFastqに変換をする必要がある。
解析に使用するシークエンスデータが複数ある場合
ファイル名を記入したリストファイルの作成してあると一括で処理できる。

まずはダウンロード。
Home - SRA - NCBI

例題としてHeLaのRNAseqを探す
検索窓に"HeLa RNAseq"と入れて検索すると312個引っかかった

試しにこのデータを使ってみる
RNA-seq of HeLa cells - SRA - NCBI
しかし、データの評価の仕方がまだよくわからない
これを選んだのはなんとなく。。。
データのアクセッションナンバーはSRR6799791で、
アドレスはここで探せる
SRA Explorer

実際に探したアドレスを使ってファイルをダウンロード

$ cd $HOME/usr/data/
$ mkdir ./sra/
$ cd ./sra/
$ wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR679/SRR6799791/SRR6799791.sra

サイズが2.9 GBあるので時間がそれなりに掛かる
"$HOME"は、チルダで書くのがいいんだろうけど
スクリプト内で"$HOME"を使っているのでこの表記で。

SRA Toolkitを入れてればコマンドも使えるようだけど
使ったことがないので割愛
SRA Toolkit 使い方公開データのダウンロードとsra fastq変換 – バイオインフォ道場 [bioinfo-Dojo]

ダウンロードしたファイルはSRA Toolkitsに含まれる
fastq-dumpを使ってfastqに変換する
SRA Toolkitは前回インストール済み
RNA-seqその１、インストール - mecobalamin’s diary

16 April 2020 追記
SRA Explorerで生成したアドレスではダウンロードできなくなっているっぽい
SRA toolkitのprefetchでダウンロードする
prefetchを実行する前にvdb-confingを実行する
Toolkit Documentation : Software : Sequence Read Archive : NCBI/NLM/NIH
vdb-configのあるディレクトリで以下のように実行する

./vdb-config -i

RAMのサイズ以外は変更していない

追記ここまで

fastq-dumpで変換しただけでは動かないコマンドがあったりする
実際Trinityが動かなかった
NCBI SRA Toolkitの使い方 - アメリエフのブログ
Trinityはアセンブルを行うツール
詳しくはまた今度。

で、sraファイルを複数置換できるようにスクリプトを書いた
まずsraファイルのリストを作成

$ ls $HOME/usr/data/sra/ > $HOME/usr/scripts/list_sra.txt

以下のスクリプトこの名前"RNAseq01_SratoFastq.sh"で保存
保存場所は$HOME/usr/scripts/

#!/bin/bash

echo "change directory to $HOME/usr/scripts/"
cd $HOME/usr/scripts/

w_dir="$HOME/usr/data"
scripts_dir="$HOME/usr/scripts"
sra_dir="${w_dir}/sra"
out_dir="${w_dir}/raw"

list_name="list_sra.txt"
list_dir="$HOME/usr/scripts"

echo "making directory"
mkdir -p ${out_dir}

echo "run fastq-dump"
for i in $(cat ${list_dir}/${list_name} | sed 's/\.sra//g')
  do
    file=${list_dir}/${i}.sra
    fastq-dump \
      --split-files \
      --defline-seq '@$sn[_$rn]/$ri' \
      --outdir ${out_dir} \
      ${sra_dir}/${i}.sra
    sed -i -e 's/_forward//' ${out_dir}/${i}_1.fastq
    sed -i -e 's/_reverse//' ${out_dir}/${i}_2.fastq
  done

"--defline-seq"でRead1とRead2の区別がつくように、
配列名の最後に/1または/2をつけている
NCBI SRA Toolkitの使い方 - アメリエフのブログ
元データがペアエンドのデータなので"--split-files"を指定する

シェルスクリプトの実行方法はこちら

$ cd $HOME/usr/scripts/
$ bash RNAseq01_SratoFastq.sh

ファイルに実行権限を与えて実行しても良い

$ cd $HOME/usr/scripts/
$ chmod 755 RNAseq01_SratoFastq.sh
$ ./RNAseq01_SratoFastq.sh

chmodの実行は最初の一回だけ行う
【 chmod 】コマンド――ファイル／ディレクトリのパーミッション（許可属性）を変更する：Linux基本コマンドTips（14） - ＠IT

結果はrawに保存される
ファイル名はSRR6799791_1.fastq及びSRR6799791_2.fastqとなる
1、2には対になるリードがそれぞれ保存されている

十分にメモリとストレージがあれば、このままでいいのだが、
使っている環境に対してファイルサイズが大きすぎるのでファイルを分割する
fastq形式は1つのリードを4行で表示するので
4行単位で分割する
ここを参考にしてスクリプトを書いた
大きいfastqファイルを分割マッピング(tophat, linux) - script of bioinformatics
fastqを分割するツールもあった
fastq / fastaの操作ツール seqkit - macでインフォマティクス

とりあえずファイルを10分割する
ファイル名は"RNAseq02_divideFastqFiles.sh"
このスクリプトを実行すると10個に分割されたfastqと
ファイルリストのファイル"list_div_SRR6799791"が作成される

#!/bin/bash

echo "change directory to $HOME/usr/data/raw/"
cd $HOME/usr/data/raw/
scripts_dir="$HOME/usr/scripts"

w_dir="$HOME/usr/data"
out_dir="${w_dir}/raw"

list_name="list_sra.txt"
list_dir="$HOME/usr/scripts"

n=10

for i in $(cat ${scripts_dir}/list_sra.txt | sed 's/\.sra//g')
  do
    div_file="div_${i}"
    echo "count row number"
    read=$(grep "\+${i}" $HOME/usr/data/raw/${i}_1.fastq | wc -l)

    sp=$(expr \( \( ${read} \/ ${n} \) + 1 \) \* 4)

    echo "making directory"
    mkdir -p $HOME/usr/data/raw/${div_file}/
    echo "read=${read}, sp=${sp}"
    echo "dividing fastq file"

    for j in {1..2}
      do
        split --numeric-suffixes=01 -l ${sp} $HOME/usr/data/raw/${i}_${j}.fastq ${div_file}_${j}_
        mv ${div_file}_${j}_* ./${div_file}
      done
    ls ${out_dir}/${div_file}/ > ${out_dir}/${div_file}/list_${div_file}.txt
  done