Mate-pair mapping with bwa

quiverでエラーを補正する前処理としてpbalign でリードデータを参照ゲノムにアラインメントする。

問題点

の2点です。(落とし穴が条件依存的なので、lpm scriptにしていません)

では、blasrに修正が必要となっておりましたが、今ではその必要はなさそうです。

前提

基本方針

手順詳細

smrtanalysisのインストール

$ lpm install smrtanalysis-centos

pbcoreのインストール

$ git clone https://github.com/PacificBiosciences/pbcore
$ cd pbcore
$ python setup.py install

smranalysis-2.0.1に付属のpbcoreではバージョンが古すぎるので、新しいものをインストールする

hdf5_toolsのインストール

BLASRのコンパイルにhdf5_toolsは必要であるのでBLASRより先に入れる。WarningやRemarkはあるが、深刻なエラーは発生しない。

BLASRのインストール

smrtanalysis-2.0.1に付属のblasrは版が古く、対応していないoptionをつけて呼び出されるエラーが発生する。

$ git clone https://github.com/PacificBiosciences/blasr
$ cd blasr
$ export HDF5INCLUDEDIR=$HOME/lcl/include
$ export HDF5LIBDIR=$HOME/lcl/lib
$ make CXX=/usr/bin/g++ -j 14
$ make install

C++のコンパイラーとしてはOS付属のものを使い、新しいもの(g++-4.8.1)を使うのを避ける。それでも、PrintMSA.cpp:18などからincludeされるhash_mapにdepricated扱いのメッセージがでる。一応コンパイルされるからよしとする。

pbh5toolsのインストール

$ git clone https://github.com/PacificBiosciences/pbh5tools
$ cd pbh5tools
$ python setup.py install

これにより、/home/tomoaki/lcl/opt/smrtanalysis/smrtanalysis-2.0.1/redist/python2.7/lib/python2.7/site-packages/ に pbtools.pbh5tools-0.75.0-py2.7-linux-x86_64.egg がインストールされるが、/lcl/opt/smrtanalysis/smrtanalysis-2.0.1/analysis/lib/python2.7/ にもより古い版のものが残っている点に注意する必要がある。

旧版ライブラリーの削除

$ rm -fr  $HOME/lcl/opt/smrtanalysis/smrtanalysis-2.0.1/analysis/lib/python2.7/pbtools.pbh5tools-0.75.0-py2.7-linux-x86_64.egg
$ rm -fr  $HOME/lcl/opt/smrtanalysis/smrtanalysis-2.0.1/analysis/lib/python2.7/pbcore-0.6.0-py2.7.egg

このライブラリーが残っていると優先的に利用されて、エラーが発生する。一方、$HOME/lcl/opt/smrtanalysis/smrtanalysis-2.0.1/analysis/lib/python2.7 をPYTHONPATHからはずすのはGenomicConsensusが見つからなくなってquiverでエラーが発生する事になる。よって削除する。

pbalignのインストール

$ git clone https://github.com/PacificBiosciences/pbalign.git
$ cd pbalign
$ python setup.py install

動作確認

PacificBiosciences社が公開している大腸菌データで動作確認を行う。

$ wget http://files.pacb.com/datasets/secondary-analysis/ecoli-k12-P4C2-20KSS/ecoliK12.tar.gz
$ tar xzvf ecoliK12.tar.gz
$ cd ecoliK12
$ cp ../U00096.3.fasta .
$ time pbalign.py --nproc 14 --forQuiver Analysis_Results/m130404_014004_sidney_c100506902550000001823076808221337_s1_p0.1.bax.h5 U00096.3.fasta all_against_celera.cmp.h5

real	3m19.342s
user	29m34.283s
sys	0m13.868s
time quiver -j 14 all_against_celera.cmp.h5 -r U00096.3.fasta -o consensus.fa -o consensus.fq -o variants.gff
real	7m28.705s
user	101m56.746s
sys	0m10.853s

U00096.3.fastaはDDBJなどから入手する。MG1665のreference genomeである。

上記実行例の--nproc 14は、qlogin した空いているnodeで行う事を想定している。jobで投入する時は --nproc $NSLOTSとすべきである。