「教本」に従ってサンプルファイルをダウンロードしてくる。
諸般の都合でDDBJのサーバではなくhttp://www.ebi.ac.uk/ena/data/view/DRR006760
こっちからsraじゃなく、fastq.gzファイルを落とす。（若干手間省略）
FastQC〜Trimmomaticの流れはこれまでのRNA-seqと何ら変わらないのだけど、HGCスパコン上でやるのは初めてで、色々手間取る。
Javaのバージョンはやっぱり1.8に統一しておくのが良さそうだな。昨日書いたaliasの設定をいじる。
FastQCにもPATHが通ってなかったので改めて.bashrcに追記。
スレッド数12でやろうとしたら
Error occurred during initialization of VM Could not reserve enough space for 3072000KB object heap
というエラーが出て動かず、スレッド数８だと動く。
ここはFastQCの方の問題なのだろうか。FastQCは実行可能なJARファイルとなっているが、これに-Xmxオプションは効くのだろうか？
次のTrimmomaticについても、javaがどうもすぐにcoreを吐いて止まってしまうので困った。どうやらやはりメモリヒープの上限設定がネックになっている様子。いっそのことと、-Xmxオプションを外してやってシステムの制限いっぱいにしてやるとうまく走り出した。
と思ったが、どうも遅い。
cpuスレッド数16を指定して実行させようとしているのにqhostで見る限りはそんなに並列化されている感じがない。
$ qhost |grep sc099i HOSTNAME ARCH NCPU NSOC NCOR NTHR NLOAD MEMTOT MEMUSE SWAPTO SWAPUS ---------------------------------------------------------------------------------------------- sc099i lx-amd64 24 2 24 24 0.19 125.2G 14.4G 2.0G 348.0K
ジョブが利用するスロット数を予め指定してやる必要があるらしい。
16コア使うなら
#$ -pe def_slot 16
このままでは1スロットあたり5.3Gメモリが要求され、トータル84.8G使うことになるらしい。
一応ノードのメモリ容量に収まるわけだな。
要するに1ノードに128Gメモリがあるので1ノードあたり24スロットあるから何も指定しないと1スロットあたり5.3Gのメモリが割り当てられるって寸法かな。
で、ログを見るとこんな感じ
Picked up JAVA_TOOL_OPTIONS: -XX:+UseSerialGC -Xmx64m -Xms32m
TrimmomaticPE: Started with arguments:
-threads 16 -phred33 -trimlog ./DRR006760/trimmed/DRR006760.trimlog ./DRR006760_1.fastq.gz ./DRR006760_2.fastq.gz ./DRR006760/trimmed/DRR006760_1.paired.fastq.gz ./DRR006760/trimmed/DRR006760_1.unpaired.fastq.gz ./DRR006760/trimmed/DRR006760_2.paired.fastq.gz ./DRR006760/trimmed/DRR006760_2.unpaired.fastq.gz TRAILING:20 MINLEN:50
Exception in thread "Thread-4" Exception in thread "Thread-2" java.lang.OutOfMemoryError: Java heap space
at java.util.Arrays.copyOfRange(Arrays.java:2694)
at java.lang.String.<init>(String.java:203)
at java.lang.StringBuilder.toString(StringBuilder.java:405)
at org.usadellab.trimmomatic.fastq.FastqSerializer.writeRecord(FastqSerializer.java:63)
at org.usadellab.trimmomatic.threading.SerializerWorker.run(SerializerWorker.java:45)
at java.lang.Thread.run(Thread.java:745)
java.lang.OutOfMemoryError: Java heap space
at java.lang.AbstractStringBuilder.<init>(AbstractStringBuilder.java:64)
at java.lang.StringBuilder.<init>(StringBuilder.java:97)
at org.usadellab.trimmomatic.fastq.FastqSerializer.writeRecord(FastqSerializer.java:51)
at org.usadellab.trimmomatic.threading.SerializerWorker.run(SerializerWorker.java:45)
at java.lang.Thread.run(Thread.java:745)なんかメモリー絡みにエラーが出ている。
やはり-Xmxオプションは必要ということか。
#!/bin/sh
#$ -S /bin/bash
#$ -cwd
#$ -l s_vmem=8G,mem_req=8G
#$ -pe def_slot 8
#処理するfastqファイルのあるディレクトリをワーキングディレクトリとする
#--------------------------------------------------------------------
#トリミングパラメータ
#Phred quality score "qscore"未満の末端を除去
qscore=20
length=50
#CPUスレッド数
cpu=8
#--------------------------------------------------------------------
#データ格納ディレクトリがなければ作る
mkdir -p ./${alias}/trimmed
#Trimmomaticでトリミング
java -Xms4g -Xmx4g -jar /usr/local/package/trimmomatic/0.36/trimmomatic-0.36.jar PE \
-threads ${cpu} -phred33 -trimlog ./${alias}/trimmed/${alias}.trimlog \
./${file1} \
./${file2} \
./${alias}/trimmed/${file1%%.*}.paired.fastq.gz ./${alias}/trimmed/${file1%%.*}.unpaired.fastq.gz \
./${alias}/trimmed/${file2%%.*}.paired.fastq.gz ./${alias}/trimmed/${file2%%.*}.unpaired.fastq.gz \
TRAILING:${qscore} \
MINLEN:${length}これでちゃんと走りきることをまずは確認しよう。
その上でスロットをもうちょっと増やしてみるかね。
あとtrimlogは無駄に巨大なので、うまく走るようになったら-trimlogのオプションを消してしまってもいいかも。
$ qstat -f |grep yc043 mjobs.q@yc043i BP 0/20/24 0.53 lx-amd64
8スロット予約したけどノードには20スロット走っているのでおそらくなにか相乗りしているようで。
で、cpu使用率が0.53だから結構いっぱいいっぱいで走っているかも。
続いてマッピング。まずは基本的なスクリプトを作ってみる。（bwa.sh）
#!/bin/sh
#$ -S /bin/bash
#$ -cwd
#$ -l s_vmem=8G,mem_req=8G
#$ -pe def_slot 8
#---------------------------------------------------
cpu=8
ref=human_g1k_v37_decoy.fasta
#---------------------------------------------------
bwa mem -t $cpu -M \
-R "@RG\tID:FLOWCELLID\tSM:${alias}\tPL:illumina\tLB:${alias}_library_1" \
${ref} \
./${alias}/trimmed/${file1%%.*}.paired.fastq.gz ./${alias}/trimmed/${file2%%.*}.paired.fastq.gz \
| samtools view -@${cpu} -bS - > ${alias}_aligned_reads.bam
#sorting & index
samtools sort -@${cpu} -m 2G ${alias}_aligned_reads.bam ${alias}_aligned_reads_sorted.bam
samtools index ${alias}_aligned_reads_sorted.bam ${alias}_aligned_reads_sorted.bam.baiこんな感じでどうかな。
なんだかうまくいかない。
[main] unrecognized command 'mem' view: invalid option -- '@' view: invalid option -- '6' [samopen] no @SQ lines in the header. [sam_read1] missing header? Abort! sort: invalid option -- '@' sort: invalid option -- '1' sort: invalid option -- '6' [bam_header_read] EOF marker is absent. The input is probably truncated.
-R "@RG\tID:FLOWCELLID\tSM:${alias}\tPL:illumina\tLB:${alias}_library_1" \オプション-Rってなんなのだ？
-R STR read group header line such as '@RG\tID:foo\tSM:bar' [null]
となっている。
どうやらGATKでの処理のときに必要となるヘッダらしく、あとでpicardで入れることも可能な模様。困ったらあとに回そう。
というかエラーの原因はここではなかった。bwaのバージョンが古くてbwa memが機能してなかった。
同様にsamtoolsも0.1.18になってたので、-@による並列化ができなかった模様。
samtools viewはbwaの標準出力を入力としてパイプで繋いでいるため引数として-（ハイフン）が入っている。
パイプで繋げないと巨大なSAMファイルができ、ディスクを圧迫する。まああとでrmすりゃいいだけだけど。
samtools sortのオプション-m 2Gはもっと大きくてもいいかも。
というかcore吐いてコケた。samtoolsももうちょっとチューニングが必要だな。
#$ -l s_vmem=8G,mem_req=8G
#$ -pe def_slot 8
samtools sort -@${cpu} -m 6G ${alias}_aligned_reads.bam ${alias}_aligned_reads_sorted.bamこの条件でとりあえず走りきったが、メモリの設定値、コア数もうちょっと煮詰めてみようか。
javaメモリのヒープサイズは大きければ大きいほどいいのか？
参考
http://otndnld.oracle.co.jp/document/products/jrockit/geninfo/diagnos/memman.html
64ビットシステムではXms, Xmxとも4GB未満に抑えた方がいいらしい。
まあ、無駄に確保するとシステムの方でメモリが足らなくなって落ちるということか？
