ABySSの実行

Trans-ABySSは、ライブラリーごとにdirectoryを分けてそれぞれの中に k[数字]というデイレクトリーが作られることを期待している。とりあえず、全ライブラリーをひとまとめに解析する予定で、ABySS/liballというディレクトリを作る。

$ mkdir ABySS
$ mkdir ABySS/liball

#!/usr/bin/ruby

#で始まる行はコメントでありプログラムとしては意味を持たない。しかし、最初の行に#!でインタープリタのパス名を書いておき実行許可を出しておくとそのファイル名で実行すると、当該インタープリタで解釈して実行させる事が可能であるので、一行目はいらなくてもとりあえず#!/usr/bin/rubyとしておく

name='Lib'
petargets = [
['lib1', '../../../s_8_2.sequence.txt'],
['lib2', '../../../s_8_4.sequence.txt'],
['lib3', '../../../s_8_5.sequence.txt'],
['lib4', '../../../s_8_6.sequence.txt'],
]

qcut=20

ここでライブラリーの名前と実際のシーケンスのファイルを書いておく。ここを、目的に合わせて書き換え

$ ruby submit.rb

とすれば一応submitできるのではないかと期待している。一応その先まで解説を書いておくこととする。../の回数は１つ下の??ディレクトリから相対指定なので、カレントディレクトリからより１組多くなる。qcut=20はquality cutoffの指定である。

上で設定した、petargetsのデータ構造をもとに、コマンドラインのlib1=../../../s_8_2.sequence.txt lib2=../../../s_8_4.sequence.txt lib3=../../../s_8_5.sequence.txt lib4=../../../s_8_6.sequence.txt というのを組み立てる。その組み立て部が

libs=[]
petargets.each do |ar|
  libs << ar[0]
end
lib=libs.join(" ")
libs=[]
petargets.each do |ar|
  libs << ar.join('=')
end

kvals=(26..64)

実行するkの範囲を指定する。(上限はABySSのコンパイル条件により制約が有る)

mpihead = '#!/bin/bash
#$ -l vf=2.5G
#$ -cwd
#$ -S /bin/bash
#$ -pe openmpi 48
#$ -q normal

source $HOME/.bashrc
export TMPDIR=/work/scratch/
echo starting
date
'

最初のsingle end assembly phaseを実行するjobのheaderを設定する。#!/bin/bashはbashのプログラムである旨の宣言。#$で始まる行はSGEのjob管理システムへの指令である。-cwdはjobを発行した時のカレントディレクトリにおいて実行すべしという旨の指令。#$ -S /bin/bashはjob fileのインタープリタとしては/bin/bashを使えと言う指令。#$ -pe openmpi 48で並列に48 core占有を宣言する。これによりopenmpiによる並行実行が可能になる。#$ -l vf=2.5Gはそれぞれのプロセスが占める事を予定する仮想記憶の量を指定している。計算nodeは24GBの主記憶を8Coreで使うので平均3GB弱であり、2.5というのは、8 process/nodeで入れる事ができるようにということになる。利用するメモリーの総量はcore数×2.5GBになるのでメモリーが足りなくて失敗する場合は-pe openmpi の後の数字を増やす。-q normalはどのキューに送るかの設定。source $HOME/.bashrcは環境変数等の初期化を不断通りにさせる宣言。export TMPDIR=/work/scratch/この環境変数を設定することによりテンポラリーファイルの作成場所が/tmpから/work/scratchに変わる。/tmpは通常容量が不足するが、/work/scratchはもう少し余裕がある。ただし十分とは限らない。/work/scratchが不足する場合はnfsでマウントしているディスクを使うことになるが、local diskでできる方が望ましい。echo starting dateによって開始時刻を記録する。

singlehead = '#!/bin/bash
#$ -l vf=11.5G
#$ -cwd
#$ -S /bin/bash
##$ -pe make 8
#$ -q normal

source $HOME/.bashrc
export TMPDIR=/work/scratch/
echo starting
date
'

ほぼ同様であるが、-pe makeを有効にしていないので、１つのslotのみを占めるjobとして扱われる。-l vf=11.5Gと大きめの容量を宣言することにより、このプロセス二つで全仮想記憶を要求し、他のプロセスが割り当てられないようにする。実際に必要なメモリーがさらに大きい場合は22Gなどの要求をしないとエラーが発生する可能性がある。逆に、もし実際のメモリー使用量が少ないようなら、vf=7.5Gなどと減らすことにより、１つのノードで同時に3~4個のjobを処理出来る。

normalhead = '#!/bin/bash
#$ -l vf=2.5G
#$ -cwd
#$ -S /bin/bash
###$ -pe openmpi 4
#$ -q normal

source $HOME/.bashrc
export TMPDIR=/work/scratch/
echo starting
date
'

こちらは、空間占有率もあまり大きくせず普通程度の2.5GBとする指定

curdir=Dir.pwd
kvals.each do |k|
  Dir.mkdir("k#{k.to_s}")
  Dir.chdir("k#{k.to_s}")

  ...

  Dir.chdir(curdir)
end

現在のカレントディレクトリを記録してから、それぞれのk値に対応するsubdirectoryに入って作業をし、記録したディレクトリに戻る事を繰り返す。

  do_se = false
  do_postse = false
  do_ali = false
  do_pe = false

  alignjobarraycount=0
  open("logk#{k}","w") do |log|
  open("postsek#{k}job", "w") do |postsejob|
  postsejob.puts normalhead
  open("alignk#{k}job", "w") do |alignjob|
  alignjob.puts singlehead
  open("pek#{k}job", "w") do |pejob|
  pejob.puts singlehead

下準備。４フェーズのjobファイルを開き、先のヘッダーを書き込む。

cmd = "abyss-pe -n lib='#{lib}' #{libs.join(' ')} OVERLAP_OPTIONS=--no-scaffold SIMPLEGRAPH_OPTIONS=--no-scaffold E=0 n=10 v=-v name=#{name} k=#{k.to_s} np=48 q=#{qcut}"

ここでabyss-peに渡す引数を全て結合した文字列を用意する。

  log.puts cmd

  IO.popen(cmd) do |io|
    io.each_line do |line|
    ...
    ...
    end
  end

これでabyss-pe を-nオプションつきで呼び出すことによりそれぞれのサブコマンドを書き出すが実行はしない。

      if line =~ /mpirun/
        open("sek#{k}job", "w")do |sejob|
          sejob.puts mpihead
          sejob.puts line
          sejob.puts "echo ending"
          sejob.puts "date"
        end
        do_se = true
      end

mpirunというのが含まれている行についてはsejobに書き出す。通常mpirunの行は１つしかないので、これはループの中でsejobのファイルを開いてヘッダからコマンドを書き出して閉じる所まで行う。

      if line =~ /^(AdjList|PopBubbles|MergeContigs|awk)/
        postsejob.puts line
        do_postse = true
      end

AdjList, PopBubbles, MergeContigs, awkのいずれかで始まる行はsejobが終了したら行う一連のしょりであるのでpostsejobに記録する。

      if line =~ /^KAligner/
        alignjobarraycount+=1
        alignjob.puts "if [ $SGE_TASK_ID -eq #{alignjobarraycount} ]; then"
        alignjob.puts line
        while line =~ /\\$/
          line = io.gets
          alignjob.puts line.sub(/-j48/,"-j8")
        end
        alignjob.puts "fi"
        do_ali = true
      end

KAlignerはもとのライブラリーをscaffoldにアラインメントするプロセスであり、ライブラリーごとに並列化できる。このためArray jobとしてつくる。また、このコマンドは\を使った継続行として続いているのでその分まで読み込んで書き込む。.sub(/-j48/,"-j8")で48CPU使用指定になっている所は8CPUに引き下げる。

      if line =~/^(abyss-joindist|Overlap|SimpleGraph|MergePaths|cat|PathOverlap|\s+\|PathConsensus)/
        pejob.puts line.sub(/-j48/,"-j8")
        do_pe = true
      end

いよいよ最終フェーズ、paired-end assemblyのフェーズです。いくつかのコマンドが有りますが、マッチしたものを書き出します。

  pejob.puts "echo ending"
  pejob.puts "date"
  end;
  alignjob.puts "echo ending"
  alignjob.puts "date"
  end;
  postsejob.puts "echo ending"
  postsejob.puts "date"
  end;
#  pejob.close
#  alignjob.close
#  postsejob.close

openで開きっぱなしでループを回していたファイルに終了前の時刻確認コマンドを加えて閉じます。open() do |file| .. end;のendで自動的にcloseされますので、明示的close呼び出しはありません。

  hold_jid = ''
  if do_se
    cmd = "qsub sek#{k}job"
    log.puts cmd
    sejobid=`#{cmd}`.split[2]
    log.puts "sejobid = #{sejobid}"
    hold_jid = "-hold_jid #{sejobid}"
  end

qsubコマンドによりSGEにjobを渡し、jobidを受け取ります。受け取ったjobidは引き続いて実行されるジョブにhod_jidオプションにわたし、実行順序を保証させます。

  if do_postse
    cmd = "qsub #{hold_jid} postsek#{k}job"
    log.puts cmd
    postsejobid = `#{cmd}`.split[2]
    log.puts "postsejobid = #{postsejobid}"
    hold_jid = "-hold_jid #{postsejobid}"
  end

  if do_ali
    cmd = "qsub #{hold_jid} -t 1-#{alignjobarraycount} alignk#{k}job"
    log.puts cmd
    alignjobid = `#{cmd}`.split[2].sub(/\..*/,"")
    log.puts "alignjobid = #{alignjobid}"
    hold_jid = "-hold_jid #{alignjobid}"
  end

アラインメントジョブはSGEに投入した時に返ってくる文字列の構造がやや違うので少し違う対処をします。

  if do_pe
    cmd = "qsub #{hold_jid} pek#{k}job"
    log.puts cmd
    pejobid = `#{cmd}`.split[2]
    log.puts "pejobid = #{pejobid}"
  end

最後のpaired-end assembly phaseですが、これは普通に投入するだけです。この先は有りませんので、hold_jidとして記録する必要はありません。

ABySSの実行 (trans ABySS用)

ディレクトリーの作成