RNA-seqその４、Fastqファイルのマージ

クリーニングしたFastqファイルを一つにまとめる
mecobalamin.hatenablog.com

plinseq-liteで評価の良かったリードだけをマージする
ファイル名に"gd"を含むファイルを使うことになる
read 1とread 2を区別する

マージしたファイルは、マージ前と同じディレクトリに保存される

#!/bin/bash

dir_scripts="$HOME/usr/scripts/"
dir_working="$HOME/usr/data/"

dir_list_sra=${dir_scripts}
name_sra="list_sra.txt"
list_filename=$(cat ${dir_list_sra}/${name_sra} | sed 's/\.sra//g')

for i in ${list_filename}
  do
    dir_div_file="${dir_working}/raw/div_${i}"
    list_div_file="list_div_${i}.txt"
    list_read_1=$(cat ${dir_div_file}/${list_div_file} | grep _1_)
    list_read_2=$(cat ${dir_div_file}/${list_div_file} | grep _2_)

    echo "merge read 1"
    :> ${dir_working}/trim_${i}/ps/ps_${i}_gd_1.fastq
    for j in ${list_read_1}
      do
        cat ${dir_working}/trim_${i}/ps/ps_${j/_1_/_}_gd_1.fastq >> \
          ${dir_working}/trim_${i}/ps/ps_${i}_gd_1.fastq
        #echo ${i}, ${j/_1_/_}
      done

    echo "merge read 2"
    :> ${dir_working}/trim_${i}/ps/ps_${i}_gd_2.fastq
    for j in ${list_read_2}
      do
        cat ${dir_working}/trim_${i}/ps/ps_${j/_2_/_}_gd_2.fastq >> \
          ${dir_working}/trim_${i}/ps/ps_${i}_gd_2.fastq
      done
  done

14 April 2019追記
複数のファイルを処理していたら
trimmomaticのあとディスクの容量が足りなくなった
いくつか中間ファイルを消していたら
rawのディレクトリを間違って消してしまって
一緒に分割ファイルのリストも無くなった。。。

次に進めず。。。。

分割ファイルを再作成するスクリプト

mkdir $HOME/usr/data/raw/
f="$HOME/usr/scripts/list_sra.txt"
fn=$(cat ${f} | sed 's/\.sra//g')
for i in ${fn}; do mkdir $HOME/usr/data/raw/div_${i}; done
for i in ${fn}; do for j in {1..2}; do for k in 01 02 03 04 05 06 07 08 09 10; do :> $HOME/usr/data/raw/div_${i}/div_${i}_${j}_${k}; done; done; done
for i in ${fn}; do ls $HOME/usr/data/raw/div_${i}/ > $HOME/usr/data/raw/div_${i}/list_div_${i}.txt; done

空のファイル群を作成raw/に作成したあとに
ファイルのリストを作成する
あまりきれいな書き方ではないが一応大丈夫っぽい
コマンドライン上で動かしていたのでfor文も一行で書いている
このままシェルスクリプトファイルにしても動くかも

追記ここまで

20 April 2019追記
別の処理中にペアエンドのリードの数が1と2で違っている雰囲気があったので確認
前にも確認してたけど念のためもう一度
あとスクリプトの記録

$ cat hogehoge_gd_1.fastq | wc -l
61570060
$ cat hogehoge_gd_2.fastq | wc -l
61570060

リードの数は同じ
4の倍数かどうかも確認
割り算だと整数で返してくるので余りを求める

$ echo $((61570060%4))
0

いくつかマージ後のファイルを調べたけど大丈夫そう
マージ前は？

$ for i in 01 02 03 04 05 06 07 08 09 10; do echo $(($(cat hogehoge_${i}_gd_1.fastq | wc -l)%4)); done

こっちも問題なさそう

追記ここまで