exome analysis やってみる2

このスクリプトだと多分数時間以上かかると思われるのでもっと並列化を目指す。
ちょうどHGCのサイトに参考になるスクリプトがあった
https://supcom.hgc.jp/internal/mediawiki/アレイジョブでbwaの重い処理を分割

手順としてはfastqファイルを分割
fastqは4行ひとまとめで1エントリーなので、4の倍数の行数で分割する必要がある。
アレイジョブで16ジョブに分割して並列化するならwc -lで行数を読み、4で割ってエントリー数を出し、expr (エントリー数）/ 16 + 1にして割り切れなくても最後のファイルだけ少なくても良しとし、もう一度4倍して1ファイル分の行数を決める。

#!/bin/sh
#$ -S /bin/bash
#$ -cwd
#$ -l s_vmem=8G,mem_req=8G
#$ -t 1:16

#解凍
gzip -d ./${alias}/trimmed/${file1%%.*}.paired.fastq.gz
gzip -d ./${alias}/trimmed/${file2%%.*}.paired.fastq.gz
#ファイルの行数からリード数を求める
line=`wc -l ./${alias}/trimmed/${file1%%.*}.paired.fastq | awk '{print $1;}'`
#4行で1リードであるファイルを16分割
read=`expr ${line} / 4 / 16 + 1`
#1ファイルあたりの行数
line2=`expr ${read} \* 4`

split -l ${line2} -d --verbose ./${alias}/trimmed/${file1%%.*}.paired.fastq read1-
split -l ${line2} -d --verbose ./${alias}/trimmed/${file2%%.*}.paired.fastq read2-

INDEX=$(( ${SGE_TASK_ID} - 1 ))
# 整数をゼロ埋めして2桁の文字列にする
LABEL=`printf %02d ${INDEX}`

READ_CHUNK_1=./${alias}/trimmed/read1-${LABEL}
READ_CHUNK_2=./${alias}/trimmed/read2-${LABEL}

# アライメントする
bwa mem ${ref} ${READ_CHUNK_1} ${READ_CHUNK_2} |samtools view -hbS - > ./result-${LABEL}.bam
samtools sort result-${LABEL}.bam result-${LABEL}

全部のアライメントが終わったら次のスクリプトを投げる
(終わったのを待って次のステップに自動的に移るにはどうしたらいいかわからない。）

もしくはqsub -hold_jid xxxxxxxxxx merge.sh
(上のジョブのIDをqstatで調べてから投げる。終わり次第始まるらしい。)

#!/bin/sh
#$ -S /bin/bash
#$ -cwd
#$ -l s_vmem=8G,mem_req=8G

#分割したbamをまとめる
LIST=`ls result-*.bam`
samtools merge ${alias}_aligned_reads_sorted.bam ${LIST}